Allink  v0.1
MatematicaAlgebra.cpp
1 #include "../include/Matematica.h"
2 #include <iostream>
3 #include <fstream>
4 void Matrice::Shout(const char * s, ...){
5 #ifdef MAT_DEBUG
6  va_list args;
7  va_start(args, s);
8  fprintf(stderr, "Matrice] ");
9  vfprintf(stderr, s, args);
10  fprintf(stderr, "\n");
11  va_end(args);
12 #else
13  return;
14 #endif
15 }
16 // Matrice::Matrice(int newNSize, int newNRow) { // the only public ctor
17 // if (newNSize <= 0) newNSize = 5;
18 // NSize = newNSize;
19 // if ((newNRow <= newNSize)&&(newNRow>0))
20 // NRow = newNRow;
21 // else
22 // NRow = newNSize;
23 // // since allocate() will first call delete[] on data:
24 // //data = new double [NSize*NSize];
25 // data = (double *)calloc(NSize*NSize,sizeof(double));
26 // NCol = NRow = NSize;
27 // // allocate();
28 // }
29 //----------------Constructors------------------------
30 Matrice::Matrice(int ExtNRow, int ExtNCol){
31  Shout("Matrice(row,col)");
32  NRow = ExtNRow;
33  NCol = ExtNCol;
34  NZed = 0;
35  NSize = NCol*NRow;
36  NDim = 2;
37  data = (double *)calloc(NSize,sizeof(double));
38 }
39 Matrice::Matrice(int ExtNRow, int ExtNCol, double *pointer){
40  Shout("Matrice(row,col)");
41  NRow = ExtNRow;
42  NCol = ExtNCol;
43  NSize = NCol*NRow;
44  NDim = 2;
45  data = (double *)malloc(sizeof(double));
46  data = pointer;
47 }
48 Matrice::Matrice(int ExtNRow){ // the only public data = (double *)calloc(NSize*NSize,sizeof(double));
49  Shout("Matrice(row=col)");
50  NRow = ExtNRow;
51  NCol = 0;
52  NZed = 0;
53  NSize = NRow;
54  NDim = 1;
55  data = (double *)calloc(NSize,sizeof(double));
56 }
57 Matrice::Matrice(int ExtNRow,int ExtNCol,int ExtNZed){
58  Shout("Matrice(row,col)");
59  NRow = ExtNRow;
60  NCol = ExtNCol;
61  NZed = ExtNZed;
62  NSize = NCol*NRow*NZed;
63  NDim = 3;
64  data = (double *)calloc(NSize,sizeof(double));
65 }
67  Shout("Matrice(SPLINE)");
68  NRow = NCol = 3;
69  NZed = 0;
70  NSize = NRow * NCol;
71  data = (double *)calloc(NSize,sizeof(double));
72  double MenoDue = - .75*Wg.a3 + 1.5*Wg.a4;
73  //MenoDue += .125*Wg.a1 - .125*Wg.a2; //O(h^4)
74  double MenoUno = -.5*Wg.a1 + Wg.a2 + 1.5*Wg.a3 - 6.*Wg.a4;
75  //MenoUno += -.125*Wg.a1 - .5*Wg.a2;
76  double Zero = Wg.a0 - 2.*Wg.a2 + 9.*Wg.a4;
77  //Zero += 0.75*Wg.a2;
78  double PiuUno = .5*Wg.a1 + Wg.a2 - 1.5*Wg.a3 - 6.*Wg.a4;
79  //PiuUno += .25*Wg.a1 + .5*Wg.a2;
80  double PiuDue = .75*Wg.a3 + 1.5*Wg.a4;
81  //PiuDue += -.125*Wg.a1 -.125*Wg.a2;
82  double Norm=0.;
83  for(int r=0;r<NCol;r++){
84  data[NRow*0+r] = MenoUno;
85  data[NRow*1+r] = Zero;
86  data[NRow*2+r] = PiuUno;
87  }
88  for(int r=0;r<NRow;r++)
89  for(int c=0;c<NCol;c++)
90  //Norm += data[NSize*r+c] > 0. ? data[NSize*r+c] : 0;
91  Norm += ASS(data[NRow*c+r]);
92  for(int r=0;r<NRow;r++)
93  for(int c=0;c<NCol;c++)
94  data[NRow*c+r] /= Norm*.5;
95  NCol = NRow;
96  NSize = NCol * NRow;
97  NDim = 2;
98 }
99 Matrice::Matrice(double *Axis,double Angle,int ExtNRow){
100  double c = cos(Angle);
101  double s = sin(Angle);
102  double t = 1. - c;
103  double Norm = 0.;
104  for(int d=0;d<3;d++){
105  Norm += SQR(Axis[d]);
106  }
107  Norm = Norm > 0. ? sqrt(Norm) : 1.;
108  double x = Axis[0]/Norm;
109  double y = Axis[1]/Norm;
110  double z = Axis[2]/Norm;
111  NDim = 2;
112  NZed = 0;
113  if(ExtNRow == 4){
114  NCol = NRow = 4;
115  NSize = NCol * NRow;
116  data = (double *)calloc(NSize,sizeof(double));
117  data[NRow*0+0] = t*x*x + c;
118  data[NRow*0+1] = t*x*y + z*s;
119  data[NRow*0+2] = t*x*z - y*s;
120  data[NRow*0+3] = 0.;
121 
122  data[NRow*1+0] = t*x*y - z*s;
123  data[NRow*1+1] = t*y*y + c;
124  data[NRow*1+2] = t*y*z + x*s;
125  data[NRow*1+3] = 0.;
126 
127  data[NRow*2+0] = t*x*z + y*s;
128  data[NRow*2+1] = t*y*z - x*s;
129  data[NRow*2+2] = t*z*z + c;
130  data[NRow*2+3] = 0.;
131 
132  data[NRow*3+0] = 0.;
133  data[NRow*3+1] = 0.;
134  data[NRow*3+2] = 0.;
135  data[NRow*3+3] = 1.;
136  }
137  else{
138  NCol = NRow = 3;
139  NSize = NCol * NRow;
140  data = (double *)calloc(NSize,sizeof(double));
141  data[NRow*0+0] = t*x*x + c;
142  data[NRow*0+1] = t*x*y + z*s;
143  data[NRow*0+2] = t*x*z - y*s;
144 
145  data[NRow*1+0] = t*x*y - z*s;
146  data[NRow*1+1] = t*y*y + c;
147  data[NRow*1+2] = t*y*z + x*s;
148 
149  data[NRow*2+0] = t*x*z + y*s;
150  data[NRow*2+1] = t*y*z - x*s;
151  data[NRow*2+2] = t*z*z + c;
152  }
153 }
155  //FIXME: the determinant is not zero!
156  NDim = 2;
157  NZed = 0;
158  if(dim == 4){
159  NCol = NRow = 4;
160  NSize = NCol * NRow;
161  data = (double *)calloc(NSize,sizeof(double));
162  data[NRow*0+0] = q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z;
163  data[NRow*0+1] = 2.*q.x*q.y + 2.*q.w*q.z;
164  data[NRow*0+2] = 2.*q.x*q.z - 2.*q.w*q.y;
165  data[NRow*0+3] = 0.;
166 
167  data[NRow*1+0] = 2.*q.x*q.y - 2.*q.w*q.z;
168  data[NRow*1+1] = q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z;
169  data[NRow*1+2] = 2.*q.y*q.z + 2.*q.w*q.x;
170  data[NRow*1+3] = 0.;
171 
172  data[NRow*2+0] = 2.*q.x*q.z + 2.*q.w*q.y;
173  data[NRow*2+1] = 2.*q.y*q.z - 2.*q.w*q.x;
174  data[NRow*2+2] = q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z;
175  data[NRow*2+3] = 0.;
176 
177  data[NRow*3+0] = 0.;
178  data[NRow*3+1] = 0.;
179  data[NRow*3+2] = 0.;
180  data[NRow*3+3] = q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z;
181  }
182  else{
183  NCol = NRow = 3;
184  NSize = NCol * NRow;
185  data = (double *)calloc(NSize,sizeof(double));
186  data[NRow*0+0] = 1. - 2.*SQR(q.y) - 2.*SQR(q.z);
187  data[NRow*0+1] = 2.*q.x*q.y + 2.*q.w*q.z;
188  data[NRow*0+2] = 2.*q.x*q.z - 2.*q.w*q.y;
189 
190  data[NRow*1+0] = 2.*q.x*q.y - 2.*q.w*q.z;
191  data[NRow*1+1] = 1. - 2.*SQR(q.x) - 2.*SQR(q.z);
192  data[NRow*1+2] = 2.*q.y*q.z + 2.*q.w*q.x;
193 
194  data[NRow*2+0] = 2.*q.x*q.z + 2.*q.w*q.y;
195  data[NRow*2+1] = 2.*q.y*q.z - 2.*q.w*q.x;
196  data[NRow*2+2] = 1. - 2.*SQR(q.x) - 2.*SQR(q.y);
197  }
198 }
199 Matrice::Matrice(double Roll,double Pitch,double Yaw,int ExtNRow){
200  //Defining Rrpy(r,p,y) := RzrRypRxy
201  double CR = cos(Roll);
202  double SR = sin(Roll);
203  double CP = cos(Pitch);
204  double SP = sin(Pitch);
205  double CY = cos(Yaw);
206  double SY = sin(Yaw);
207  NDim = 2;
208  NZed = 0;
209  if(ExtNRow == 4){
210  NCol = NRow = 4;
211  NSize = NCol * NRow;
212  data = (double *)calloc(NSize,sizeof(double));
213  data[NRow*0+0] = CY*CP;
214  data[NRow*0+1] = SY*CP;
215  data[NRow*0+2] = -SP;
216  data[NRow*0+3] = 0.;
217 
218  data[NRow*1+0] = CY*SP*SR - SY*CR;
219  data[NRow*1+1] = SY*SP*SR + CY*CR;
220  data[NRow*1+2] = CP*SR;
221  data[NRow*1+3] = 0.;
222 
223  data[NRow*2+0] = CY*SP*CR + SY*SR;
224  data[NRow*2+1] = SY*SP*CR - CY*SR;
225  data[NRow*2+2] = CP*CR;
226  data[NRow*2+3] = 0.;
227 
228  data[NRow*3+0] = 0.;
229  data[NRow*3+1] = 0.;
230  data[NRow*3+2] = 0.;
231  data[NRow*3+3] = 1.;
232  // Other definition of the angles (to check)
233  // data[NRow*0+0] = CR*CP*CY - SR*SY;
234  // data[NRow*0+1] = CR*CP*SY - SR*CY;
235  // data[NRow*0+2] = -CR*SP;
236  // data[NRow*0+3] = 0.;
237 
238  // data[NRow*1+0] = -SR*CP*CY - CR*SY;
239  // data[NRow*1+1] = -SR*CP*SY + CR*CY;
240  // data[NRow*1+2] = SR*SP;
241  // data[NRow*1+3] = 0.;
242 
243  // data[NRow*2+0] = SP*CY;
244  // data[NRow*2+1] = SP*SY;
245  // data[NRow*2+2] = CP;
246  // data[NRow*2+3] = 0.;
247 
248  // data[NRow*3+0] = 0.;
249  // data[NRow*3+1] = 0.;
250  // data[NRow*3+2] = 0.;
251  // data[NRow*3+3] = 1.;
252  }
253  if(ExtNRow == 3){
254  NCol = NRow = 3;
255  NSize = NCol * NRow;
256  data = (double *)calloc(NSize,sizeof(double));
257  data[NRow*0+0] = CR*CP;
258  data[NRow*0+1] = SR*CP;
259  data[NRow*0+2] = -SP;
260 
261  data[NRow*1+0] = CR*SP*SY - SR*CY;
262  data[NRow*1+1] = SR*SP*SY + CR*SY;
263  data[NRow*1+2] = CP*SY;
264 
265  data[NRow*2+0] = CR*SP*CY+SR*SY;
266  data[NRow*2+1] = SR*SP*CY-CR*SY;
267  data[NRow*2+2] = CR*CY;
268  }
269 }
270 Matrice::Matrice(double *M,int Row,int Col){
271  NRow = Row;
272  NCol = Col;
273  NZed = 0;
274  NSize = NRow*NCol;
275  //data = M;
276  NDim = 2;
277  data = (double *)calloc(NSize,sizeof(double));
278  for(int c=0;c<NCol;c++)
279  for(int r=0;r<NRow;r++)
280  data[NRow*c+r] = M[NRow*c+r];
281 }
283  Shout("Matrice(SPLINE,dim)");
284  NCol = NRow = 5;
285  NZed = 0;
286  NSize = NCol * NRow;
287  data = (double *)calloc(NSize,sizeof(double));
288  FillDiffOperator(Wg,Dim);
289 }
291  // delete[] data;
292  free(data);
293 }
294 void Matrice::RandomFill(double Max){
295  Shout("RandomFill");
296  for(int r=0;r<NRow;r++)
297  for(int c=0;c<NCol;c++)
298  data[NRow*c+r] = Max*drand48();//Casuale();
299 }
301  Shout("Matrice(SPLINE,dim)");
302  double MenoDue = - .75*Wg.a3 + 2.*Wg.a4;//1.5
303  //MenoDue += .125*Wg.a1 - .125*Wg.a2; //O(h^4)
304  MenoDue -= 1./12.*Wg.a2; //O(h^4)
305  double MenoUno = -.5*Wg.a1 + Wg.a2 + 1.5*Wg.a3 - 8.*Wg.a4;//6
306  //MenoUno += -.125*Wg.a1 - .5*Wg.a2;
307  MenoUno += 1./3.*Wg.a2; //O(h^4)
308  double Zero = Wg.a0 - 2.*Wg.a2 + 12.*Wg.a4;//9
309  //Zero += 0.75*Wg.a2;
310  Zero -= 0.5*Wg.a2; //O(h^4)
311  double PiuUno = .5*Wg.a1 + Wg.a2 - 1.5*Wg.a3 - 8.*Wg.a4;
312  //PiuUno += .25*Wg.a1 + .5*Wg.a2;
313  PiuUno += 1./3.*Wg.a2;
314  double PiuDue = .75*Wg.a3 + 2.*Wg.a4;
315  //PiuDue += -.125*Wg.a1 -.125*Wg.a2;
316  PiuDue -= 1./12.*Wg.a2; //O(h^4)
317  double Norm=0.;
318  if (Dim == 0){
319  for(int r=0;r<NCol;r++){
320  data[NRow*0+r] = MenoDue;
321  data[NRow*1+r] = MenoUno;
322  data[NRow*2+r] = Zero;
323  data[NRow*3+r] = PiuUno;
324  data[NRow*4+r] = PiuDue;
325  }
326  }
327  else if(Dim == 1){
328  int r=2;
329  data[NRow*0+r] = MenoDue;
330  data[NRow*1+r] = MenoUno;
331  data[NRow*2+r] = Zero;
332  data[NRow*3+r] = PiuUno;
333  data[NRow*4+r] = PiuDue;
334  }
335  else if (Dim == 2){
336  int r=2;
337  data[NRow*0+r] = MenoDue;
338  data[NRow*1+r] = MenoUno;
339  data[NRow*2+r] = Zero;
340  data[NRow*3+r] = PiuUno;
341  data[NRow*4+r] = PiuDue;
342  int c=2;
343  data[NRow*c+0] = MenoDue;
344  data[NRow*c+1] = MenoUno;
345  data[NRow*c+2] += Zero;
346  data[NRow*c+3] = PiuUno;
347  data[NRow*c+4] = PiuDue;
348  }
349  else if (Dim == 3){
350  for(int r=0;r<NCol;r++){
351  data[NRow*0+r] = MenoDue;
352  data[NRow*1+r] = MenoUno;
353  data[NRow*2+r] = Zero;
354  data[NRow*3+r] = PiuUno;
355  data[NRow*4+r] = PiuDue;
356  }
357  }
358 }
360  if(NRow != 5){
361  printf("I don't know the expresion for matrices of NRow != 5\n");
362  return ;
363  }
364  if(NCol != 5){
365  printf("I don't know the expresion for matrices of NRow != 5\n");
366  return ;
367  }
368  Set(0,0,2.); Set(0,1,4.); Set(0,2,5.); Set(0,3,4.); Set(0,4,2.);
369  Set(1,0,4.); Set(1,1,9.); Set(1,2,12.); Set(1,3,9.); Set(1,4,4.);
370  Set(2,0,5.); Set(2,1,12.); Set(2,2,15.); Set(2,3,12.); Set(2,4,5.);
371  Set(3,0,4.); Set(3,1,9.); Set(3,2,12.); Set(3,3,9.); Set(3,4,4.);
372  Set(4,0,2.); Set(4,1,4.); Set(4,2,5.); Set(4,3,4.); Set(4,4,2.);
373  Multiply(1./159.);
374 }
375 void Matrice::FillGaussian(double Sigma,double CutOff){
376  int Half = NRow/2;
377  double Norm = 0.;
378  if(NDim == 3){
379  for(int r=0;r<NRow;r++){
380  double x = CutOff*(r-Half)/(double)NRow;
381  for(int c=0;c<NCol;c++){
382  double y = CutOff*(c-Half)/(double)NCol;
383  for(int q=0;q<NCol;q++){
384  double z = CutOff*(q-Half)/(double)NCol;
385  double r2 = SQR(x) + SQR(y) + SQR(z);
386  double Gauss = exp(-r2*.5/SQR(Sigma));
387  data[(NRow*c+r)*NCol+q] = Gauss;
388  Norm += Gauss;
389  }
390  }
391  }
392  }
393  else if(NDim == 2){
394  for(int r=0;r<NRow;r++){
395  double x = CutOff*(r-Half)/(double)NRow;
396  for(int c=0;c<NCol;c++){
397  double y = CutOff*(c-Half)/(double)NCol;
398  double r2 = SQR(x) + SQR(y);
399  double Gauss = exp(-r2*.5/SQR(Sigma));
400  data[NRow*c+r] = Gauss;
401  Norm += Gauss;
402  }
403  }
404  }
405  else if(NDim == 1){
406  for(int r=0;r<NRow;r++){
407  double x = CutOff*(r-Half)/(double)NRow;
408  double r2 = SQR(x);
409  double Gauss = exp(-r2*.5/SQR(Sigma));
410  data[r] = Gauss;
411  Norm += Gauss;
412  }
413  }
414  Multiply(1./Norm);
415 }
417  if(NRow != 5 && NCol != 5){
418  printf("Invalid number of rows and cols %dx%d != 7x7\n",NRow,NCol);
419  return;
420  }
421  data[NRow*0+0] = 1;
422  data[NRow*0+1] = 4;
423  data[NRow*0+2] = 7;
424  data[NRow*0+3] = 4;
425  data[NRow*0+4] = 1;
426 
427  data[NRow*1+0] = 4;
428  data[NRow*1+1] = 16;
429  data[NRow*1+2] = 26;
430  data[NRow*1+3] = 16;
431  data[NRow*1+4] = 4;
432 
433  data[NRow*2+0] = 7;
434  data[NRow*2+1] = 26;
435  data[NRow*2+2] = 41;
436  data[NRow*2+3] = 26;
437  data[NRow*2+4] = 7;
438 
439  data[NRow*3+0] = 4;
440  data[NRow*3+1] = 16;
441  data[NRow*3+2] = 26;
442  data[NRow*3+3] = 16;
443  data[NRow*3+4] = 4;
444 
445  data[NRow*4+0] = 1;
446  data[NRow*4+1] = 4;
447  data[NRow*4+2] = 7;
448  data[NRow*4+3] = 4;
449  data[NRow*4+4] = 1;
450 
451  Multiply(1./273.);
452 }
453 
454 
455 //---------------Manage-entries-------------------
456 double Matrice ::Val(int row){
457 #ifdef MATR_DEBUG
458  if ( (row>=NRow) || (row<0) )
459  { printf("Values %d %d out of range %d %d\n",row,NRow);
460  return 1.; }
461 #endif
462  return data[row];
463 }
464 double Matrice ::Val(int row,int col){
465 #ifdef MATR_DEBUG
466  if ( (row>=NRow) || (col>=NCol)
467  || (row<0) || (col<0) )
468  { printf("Values %d %d out of range %d %d\n",row,col,NRow,NCol);
469  return 1.; }
470 #endif
471  return data[ NRow*col+row ];
472 }
473 double Matrice ::Val(int row,int col,int zed){
474  return data[ (NRow*col+row)*NCol+zed ];
475 }
476 bool Matrice ::Set(int row, int column, double newvalue) {
477 #ifdef MATR_DEBUG
478  if ( (row >= NRow) || (column >= NCol)
479  || (row<0) || (column<0) ) return false;
480 #endif
481  data[ NRow*column + row ] = newvalue;
482  return true;
483 }
484 bool Matrice ::Add(int row, int column, double Value) {
485 #ifdef MATR_DEBUG
486  if ( (row >= NRow) || (column >= NCol)
487  || (row<0) || (column<0) ) return false;
488 #endif
489  data[ NRow*column + row ] += Value;
490  return true;
491 }
493  if(NDim == 3){
494  for(int q=0;q<NZed;q++){
495  printf("%d)\n",q);
496  for(int r=0;r<NRow;r++){
497  printf("|");
498  for(int c=0;c<NCol;c++){
499  printf("%.2g ",data[(NRow*c+r)*NCol+q]);
500  }
501  printf("|\n");
502  }
503  printf("\n");
504  }
505  return;
506  }
507  else if(NDim == 1){
508  printf("|");
509  for(int r=0;r<NRow;r++){
510  printf("%.2g ",data[r]);
511  }
512  printf("|\n");
513  printf("\n");
514  return;
515  }
516  for(int r=0;r<NRow;r++){
517  printf("|");
518  for(int c=0;c<NCol;c++){
519  printf("%.2g ",data[NRow*c+r]);
520  }
521  printf("|\n");
522  }
523  printf("\n");
524 }
525 //---------------inversion-(-stolen,-to-be-checked-)---------
527  int worstdiagonal = 0;
528  double maxunitydeviation = 0.0;
529  double currentunitydeviation;
530  for ( int i = 0; i < NRow; i++ ) {
531  currentunitydeviation = data[i+i*NRow] - 1.;
532  if ( currentunitydeviation < 0.0) currentunitydeviation *= -1.;
533  if ( currentunitydeviation > maxunitydeviation ) {
534  maxunitydeviation = currentunitydeviation;
535  worstdiagonal = i;
536  }
537  }
538  int worstoffdiagonalrow = 0;
539  int worstoffdiagonalcolumn = 0;
540  double maxzerodeviation = 0.0;
541  double currentzerodeviation ;
542  for ( int i = 0; i < NRow; i++ ) {
543  for ( int j = 0; j < NRow; j++ ) {
544  if ( i == j ) continue; // we look only at non-diagonal terms
545  currentzerodeviation = data[i+j*NRow];
546  if ( currentzerodeviation < 0.0) currentzerodeviation *= -1.0;
547  if ( currentzerodeviation > maxzerodeviation ) {
548  maxzerodeviation = currentzerodeviation;
549  worstoffdiagonalrow = i;
550  worstoffdiagonalcolumn = j;
551  }
552 
553  }
554  }
555  cout << "Worst diagonal value deviation from unity: "
556  << maxunitydeviation << " at row/column " << worstdiagonal << endl;
557  cout << "Worst off-diagonal value deviation from zero: "
558  << maxzerodeviation << " at row = " << worstoffdiagonalrow
559  << ", column = " << worstoffdiagonalcolumn << endl;
560 }
562 // NRow = left.getNRow();
563 // if ( NSize < left.getNRow() ) {
564 // NSize = left.getNRow();
565 // allocate();
566 // }
567  for ( int i = 0; i < NRow; i++ )
568  for ( int j = 0; j < NRow; j++ ) {
569  double sum = 0.0;
570  double leftvalue, rightvalue;
571  bool success;
572  for (int c = 0; c < NRow; c++) {
573  left.getvalue(i,c,leftvalue,success);
574  right.getvalue(c,j,rightvalue,success);
575  sum += leftvalue * rightvalue;
576  }
577  Set(i,j,sum);
578  }
579 }
581 // NRow = source.getNRow();
582 // if ( NSize < source.getNRow() ) {
583 // NSize = source.getNRow();
584 // allocate();
585 // }
586  for ( int i = 0; i < NRow; i++ )
587  for ( int j = 0; j < NRow; j++ ) {
588  double value;
589  bool success;
590  source.getvalue(i,j,value,success);
591  data[i+j*NRow] = value;
592  }
593 }
594 void Matrice ::setNRow(int newNRow) {
595 // if ( newNRow > NSize )
596 // {
597 // NSize = newNRow ; // * 2; // wastes memory but saves
598 // // time otherwise required for
599 // // operation new[]
600 // allocate();
601 // }
602 // if (newNRow >= 0) NRow = newNRow;
603 }
604 int Matrice ::getNRow() { return NRow; }
605 void Matrice ::getvalue(int row, int column, double& returnvalue, bool& success) {
606 #ifdef MATR_DEBUG
607  if ( (row>=NRow) || (column>=NCol)
608  || (row<0) || (column<0) )
609  { success = false;
610  return; }
611 #endif
612  returnvalue = data[ row + column*NRow ];
613  success = true;
614 }
616 #ifdef MATR_DEBUG //CHECKIT
617  if (NRow <= 0) return; // sanity check
618  if (NRow == 1) return; // must be of dimension >= 2
619 #endif
620  for (int i=1; i < NRow; i++) data[i] /= data[0]; // normalize row 0
621  for (int i=1; i < NRow; i++) {
622  for (int j=i; j < NRow; j++) { // do a column of L
623  double sum = 0.0;
624  for (int k = 0; k < i; k++)
625  sum += data[j+k*NRow] * data[k+i*NRow];
626  data[j+i*NRow] -= sum;
627  }
628  if (i == NRow-1) continue;
629  for (int j=i+1; j < NRow; j++) { // do a row of U
630  double sum = 0.0;
631  for (int k = 0; k < i; k++)
632  sum += data[i+k*NRow]*data[k+j*NRow];
633  data[i+j*NRow] =
634  (data[i+j*NRow]-sum) / data[i+i*NRow];
635  }
636  }
637  for ( int i = 0; i < NRow; i++ ) // invert L
638  for ( int j = i; j < NRow; j++ ) {
639  double x = 1.0;
640  if ( i != j ) {
641  x = 0.0;
642  for ( int k = i; k < j; k++ )
643  x -= data[j+k*NRow]*data[k+i*NRow];
644  }
645  data[j+i*NRow] = x / data[j+j*NRow];
646  }
647  for ( int i = 0; i < NRow; i++ ) // invert U
648  for ( int j = i; j < NRow; j++ ) {
649  if ( i == j ) continue;
650  double sum = 0.0;
651  for ( int k = i; k < j; k++ )
652  sum += data[k+j*NRow]*( (i==k) ? 1.0 : data[i+k*NRow] );
653  data[i+j*NRow] = -sum;
654  }
655  for ( int i = 0; i < NRow; i++ ) // final inversion
656  for ( int j = 0; j < NRow; j++ ) {
657  double sum = 0.0;
658  for ( int k = ((i>j)?i:j); k < NRow; k++ )
659  sum += ((j==k)?1.0:data[j+k*NRow])*data[k+i*NRow];
660  data[j+i*NRow] = sum;
661  }
662 };
663 //-------------Other-operations-----------------
664 int Matrice ::Solve(double *Known,double *UnKnown){
665  Shout("Solve");
666 #ifdef MATR_DEBUG
667 #endif
668  Matrice Lower(NCol,NCol);
669  Matrice Upper(NCol,NCol);
670  double *AlmostKnown = (double *)calloc(NCol,sizeof(double));
671  if(!AlmostKnown) {printf("Not allocated\n");return 0;}
672  /* assegna i valori della diagonale di L */
673  for(int r = 0; r < NRow ; r ++){
674  Lower.Set(r,r,1.);
675  }
676  /* calcola gli elementi fuori della diagonale */
677  for(int j = 0; j < NCol ; j ++){
678  for(int i = 0; i <= j ; i++){
679  Upper.Set(i,j,Val(i,j));
680  for (int k = 0; k <= i - 1; k++) {
681  Upper.Add(i,j,-Lower.Val(i,k)*Upper.Val(k,j));
682  }
683  }
684  for (int i = j + 1; i < NCol ; i++) {
685  Lower.Set(i,j,Val(i,j));
686  for (int k = 0; k <= j - 1; k ++) {
687  Lower.Add(i,j,-Lower.Val(i,k)*Upper.Val(k,j));
688  }//FIXME
689  double Diag = POS(Upper.Val(j,j)) > 0. ? Lower.Val(i,j)/Upper.Val(j,j) : 1.;
690  Lower.Set(i,j,Diag);
691  }
692  }//Solve the Uz = y
693  //Upper = NULL;
694  for(int r=0;r<NRow;r++){
695  AlmostKnown[r] = Known[r];
696  for(int c=r-1;c>=0;c--){
697  AlmostKnown[r] -= AlmostKnown[c]*Lower.Val(r,c);
698  //printf("Coeff[%d][%d]=%lf %lf\n",r,c,UnKnown[r],Lower.Val(r,c));
699  }
700  if(POS(Lower.Val(r,r))> 0.)
701  AlmostKnown[r] /= POS(Lower.Val(r,r))>0. ? Lower.Val(r,r) : 1.;
702  else
703  AlmostKnown[r] = 0.;
704  //printf("Coeff[%d][%d]=%lf %lf\n",r,r,AlmostKnown[r],Lower.Val(r,r));
705  }
706  //Lower.Print();
707  for(int r=NRow-1;r>=0;r--){
708  UnKnown[r] = AlmostKnown[r];
709  for(int c=r+1;c<NCol;c++){
710  UnKnown[r] -= UnKnown[c]*Upper.Val(r,c);
711  // printf("Coeff[%d][%d]=%lf %lf\n",r,c,AlmostKnown[r],Upper.Val(r,c));
712  }
713  if(POS(Upper.Val(r,r))> 0.)
714  UnKnown[r] /= Upper.Val(r,r);
715  else
716  UnKnown[r] = 0.;
717  //printf("Coeff[%d][%d]=%lf %lf\n",r,r,UnKnown[r],Upper.Val(r,r));
718  }//finish solving Lx = z
719  //Upper.Print();
720  free(AlmostKnown);
721  return 1;
722 };
723 void Matrice::Apply(double *Known,double *UnKnown){
724  Shout("Apply");
725  for(int r=0;r<NRow;r++){
726  for(int c=0;c<NCol;c++){
727  UnKnown[r] += Known[c]*data[r+c*NRow];
728  }
729  }
730 };
731 double Matrice::Det(){
732  Shout("Det");
733  Matrice Lower(NCol);
734  Matrice Upper(NCol);
735  /* assegna i valori della diagonale di L */
736  for (int r = 0; r < NCol ; r ++) {
737  Lower.Set(r,r,1.);
738  }
739  /* calcola gli elementi fuori della diagonale */
740  for (int j = 0; j < NCol ; j ++) {
741  for (int i = 0; i <= j ; i++) {
742  Upper.Set(i,j,Val(i,j));
743  for (int k = 0; k <= i - 1; k++) {
744  Upper.Add(i,j,-Lower.Val(i,k)*Upper.Val(k,j));
745  }
746  }
747  for (int i = j + 1; i < NCol ; i++) {
748  Lower.Set(i,j,Val(i,j));
749  for (int k = 0; k <= j - 1; k ++) {
750  Lower.Add(i,j,-Lower.Val(i,k)*Upper.Val(k,j));
751  }
752  Lower.Set(i,j,Lower.Val(i,j)/Upper.Val(j,j));
753  }
754  }
755  double Det = 1.;
756  for(int r=0;r<NCol;r++)
757  Det *= Upper.Val(r,r);
758  return Det;
759 };
761  Shout("Transpose");
762  double *Temp = (double *)calloc(NCol*NRow,sizeof(double));
763  for(int r=0;r<NRow;r++){
764  for(int c=0;c<NCol;c++){
765  Temp[r+c*NRow] = data[r+c*NRow];
766  }
767  }
768  for(int r=0;r<NRow;r++){
769  for(int c=0;c<NCol;c++){
770  data[r+c*NRow] = Temp[c+r*NRow];
771  }
772  }
773  free(Temp);
774 }
776  Shout("Clear");
777  for(int r=0;r<NRow;r++){
778  for(int c=0;c<NCol;c++){
779  data[r+c*NRow] = 0.;
780  }
781  }
782 }
784  Shout("Normalize");
785  double Max=0.;
786  double Min=0.;
787  for(int r=0;r<NRow;r++){
788  for(int c=0;c<NCol;c++){
789  if(Max < data[r+c*NRow])
790  Max = data[r+c*NRow];
791  if(Min > data[r+c*NRow])
792  Min = data[r+c*NRow];
793  }
794  }
795  if(fabs(Min) > 0.0e-5){
796  for(int r=0;r<NRow;r++)
797  for(int c=0;c<NCol;c++)
798  data[r+c*NRow] = data[r+c*NRow]/(Max-Min)-Min;
799 
800  }
801  else {
802  for(int r=0;r<NRow;r++)
803  for(int c=0;c<NCol;c++)
804  data[r+c*NRow] = data[r+c*NRow]/Max;
805  }
806 }
807 void Matrice::Multiply(double Val){
808  for(int s=0;s<NSize;s++)
809  data[s] *= Val;
810  // for(int r=0;r<NRow;r++)
811  // for(int c=0;c<NCol;c++)
812  // data[r+c*NRow] = data[r+c*NRow]*Val;
813 }
815  for(int r=0;r<NRow;r++)
816  for(int c=0;c<NCol;c++)
817  B->Set(r,c,data[r+c*NRow]);
818 }
820  if(NCol != A.NCol || NRow != A.NCol){
821  printf("Incompatible matrices. Dim: %d/%d\n",NSize,A.NSize);
822  return 0;
823  }
824  Matrice Resp(NRow,NCol);
825  for(int r=0;r<NRow;r++)
826  for(int c=0;c<NCol;c++)
827  Resp.Set(r,c,Val(r,c)+A.Val(r,c));
828  return Resp;
829 }
831  if(NCol != A.NRow || NRow != A.NCol){
832  printf("Incompatible matrices. Dim: %d/%d\n",NSize,A.NSize);
833  return 0;
834  }
835  Matrice Resp(NRow,A.NCol);
836  for(int r=0;r<NRow;r++)
837  for(int c=0;c<A.NCol;c++){
838  double Temp=0.;
839  for(int i=0;i<NRow;i++)
840  //Temp += Val(r,i)*A.Val(i,c);
841  Temp += Resp.data[r+i*NRow]*A.data[i+c*NRow];
842  Resp.Set(r,c,Temp);
843  }
844  Resp.Print();
845  return Resp;
846 }
848  if(NRow != A.NRow || NCol != A.NCol){
849  printf("Incompatible matrices. Dim: %d/%d\n",NSize,A.NSize);
850  return 0;
851  }
852  for(int r=0;r<NRow;r++)
853  for(int c=0;c<NCol;c++)
854  Set(r,c,A.Val(r,c));
855  return *this;
856 }
857 Matrice Matrice::operator* (const double& Molt)const{
858  for(int r=0;r<NRow;r++)
859  for(int c=0;c<NCol;c++)
860  data[r+c*NRow] *= Molt;
861  return *this;
862 }
864  assert( NCol == A.NRow || NRow == A.NCol);
865  printf("Still to be written!!!\n");
866  Matrice B(NCol,NRow);
867  return B;
868 }
870  if(NCol != u.NDim){
871  printf("Incompatible matrices. Dim: %d/%d\n",NSize,u.NDim);
872  return 0;
873  }
874  Vettore Resp(NRow);
875  for(int r=0;r<NRow;r++){
876  double Temp=0.;
877  for(int c=0;c<NCol;c++)
878  Temp += data[r+c*NRow]*u.x[c];
879  Resp.x[r] = Temp;
880  }
881  return Resp;
882 }
884  if(B.NCol != A.NRow || B.NRow != A.NCol){
885  printf("Incompatible matrices. Dim: %d/%d\n",B.NSize,A.NSize);
886  return ;
887  }
888  for(int r=0;r<A.NRow;r++)
889  for(int c=0;c<B.NCol;c++){
890  double Temp=0.;
891  for(int i=0;i<A.NRow;i++)
892  Temp += A.Val(r,i)*B.Val(i,c);
893  Set(r,c,Temp);
894  }
895 }
897  if(NCol != A.NRow || NRow != A.NCol){
898  printf("Incompatible matrices. Dim: %d/%d\n",NSize,A.NSize);
899  return ;
900  }
901  for(int r=0;r<A.NRow;r++)
902  for(int c=0;c<NCol;c++){
903  double Temp=0.;
904  for(int i=0;i<A.NRow;i++)
905  Temp += A.Val(r,i)*Val(i,c);
906  Set(r,c,Temp);
907  }
908 }
910  Vettore Resp(NRow);
911  for(int r=0;r<NRow;r++){
912  double Temp=0.;
913  for(int c=0;c<NCol;c++)
914  Temp += A.data[r+c*NRow]*v.x[c];
915  Resp.x[r] = Temp;
916  }
917  return Resp;
918 }
920  Vettore Resp(NRow);
921  for(int r=0;r<NRow;r++){
922  double Temp=0.;
923  for(int c=0;c<NCol;c++)
924  Temp += data[r+c*NRow]*v.x[c];
925  Resp.x[r] = Temp;
926  }
927  return Resp;
928 }
930  for(int r=0;r<NRow;r++){
931  double Temp=0.;
932  for(int c=0;c<NCol;c++)
933  Temp += data[r+c*NRow]*v.x[c];
934  u.x[r] = Temp;
935  }
936 }
937 void Matrice::ConvoluteMatrix(double *Plot,int NGrid,int NDim,int IfMinImConv){
938  if(NDim == 1){
939  if(!IfMinImConv) ConvoluteMatrix1(Plot,NGrid);
940  else ConvoluteMatrix1MinImConv(Plot,NGrid);
941  }
942  else if(NDim == 2){
943  if(!IfMinImConv) ConvoluteMatrix2(Plot,NGrid);
944  else ConvoluteMatrix2MinImConv(Plot,NGrid);
945  }
946  else if(NDim == 3){
947  ConvoluteMatrix3(Plot,NGrid);
948  }
949 }
950 //without minimum image convention
951 void Matrice::ConvoluteMatrix1(double *Plot,int NGrid){
952  double *Plot2 = (double *)calloc(NGrid,sizeof(double));
953  int NMat2 = (int)pNRow()/2;
954  for(int gx=0;gx<NGrid;gx++){
955  if(gx <= NMat2 || gx >= NGrid - NMat2){
956  Plot2[gx] = Plot[gx];
957  continue;
958  }
959  for(int mx=0;mx<pNRow();mx++){
960  int g1x = gx + mx - NMat2;
961  if(g1x >= NGrid){
962  Plot2[gx] = Plot[gx];
963  break;
964  }
965  if(g1x < 0){
966  Plot2[gx] = Plot[gx];
967  break;
968  }
969  Plot2[gx] += Plot[g1x]*Val(mx);
970  }
971  }
972  for(int gx=0;gx<NGrid;gx++){
973  Plot[gx] = Plot2[gx];
974  }
975  free(Plot2);
976 }
977 //with minimum image convention
978 void Matrice::ConvoluteMatrix1MinImConv(double *Plot,int NGrid){
979  double *Plot2 = (double *)calloc(NGrid,sizeof(double));
980  int NMat2 = (int)pNRow()/2;
981  for(int gx=0;gx<NGrid;gx++){
982  if(gx <= NMat2 || gx >= NGrid - NMat2){
983  Plot2[gx] = Plot[gx];
984  continue;
985  }
986  for(int mx=0;mx<pNRow();mx++){
987  int g1x = gx + mx - NMat2;
988  if(g1x >= NGrid) g1x -= NGrid;
989  if(g1x < 0) g1x + NGrid;
990  Plot2[gx] += Plot[g1x]*Val(mx);
991  }
992  }
993  for(int gx=0;gx<NGrid;gx++){
994  Plot[gx] = Plot2[gx];
995  }
996  free(Plot2);
997 }
999 void Matrice::ConvoluteMatrix2(double *Plot,int NGrid){
1000  double *Plot2 = (double *)calloc(SQR(NGrid),sizeof(double));
1001  int NMat2 = (int)pNCol()/2;
1002  for(int gx=0;gx<NGrid;gx++){
1003  for(int gy=0;gy<NGrid;gy++){
1004  for(int mx=0;mx<pNRow();mx++){
1005  int g1x = gx + mx - NMat2;
1006  if(g1x >= NGrid){
1007  Plot2[gx*NGrid+gy] = Plot[gx*NGrid+gy];
1008  break;
1009  }
1010  if(g1x < 0) {
1011  Plot2[gx*NGrid+gy] = Plot[gx*NGrid+gy];
1012  break;
1013  }
1014  for(int my=0;my<pNCol();my++){
1015  int g1y = gy + my - NMat2;
1016  if(g1y >= NGrid){
1017  Plot2[gx*NGrid+gy] = Plot[gx*NGrid+gy];
1018  break;
1019  }
1020  if(g1y < 0){
1021  Plot2[gx*NGrid+gy] = Plot[gx*NGrid+gy];
1022  break;
1023  }
1024  Plot2[gx*NGrid+gy] += Plot[g1x*NGrid+g1y]*Val(mx,my);
1025  }
1026  }
1027  }
1028  }
1029  for(int gx=0;gx<NGrid;gx++){
1030  for(int gy=0;gy<NGrid;gy++){
1031  Plot[gx*NGrid+gy] = Plot2[gx*NGrid+gy];
1032  }
1033  }
1034  free(Plot2);
1035 }
1037 void Matrice::ConvoluteMatrix2MinImConv(double *Plot,int NGrid){
1038  double *Plot2 = (double *)calloc(SQR(NGrid),sizeof(double));
1039  int NMat2 = (int)pNCol()/2;
1040  for(int gx=0;gx<NGrid;gx++){
1041  for(int gy=0;gy<NGrid;gy++){
1042  for(int mx=0;mx<pNRow();mx++){
1043  int g1x = gx + mx - NMat2;
1044  if(g1x >= NGrid) continue;//g1x -= NGrid;
1045  if(g1x < 0) continue;//g1x + NGrid;
1046  for(int my=0;my<pNCol();my++){
1047  int g1y = gy + my - NMat2;
1048  if(g1y >= NGrid) continue;//g1y -= NGrid;
1049  if(g1y < 0) continue;//g1y + NGrid;
1050  Plot2[gx*NGrid+gy] += Plot[g1x*NGrid+g1y]*Val(mx,my);
1051  }
1052  }
1053  }
1054  }
1055  for(int gx=0;gx<NGrid;gx++){
1056  for(int gy=0;gy<NGrid;gy++){
1057  Plot[gx*NGrid+gy] = Plot2[gx*NGrid+gy];
1058  }
1059  }
1060  free(Plot2);
1061 }
1062 //with minimum image convention
1063 void Matrice::ConvoluteMatrix3(double *Plot,int NGrid){
1064  double *Plot2 = (double *)calloc(CUBE(NGrid),sizeof(double));
1065  int NMat2 = (int)pNCol()/2;
1066  printf("NMat2 %d\n",NMat2);
1067  for(int gx=0;gx<NGrid;gx++){
1068  for(int gy=0;gy<NGrid;gy++){
1069  for(int gz=0;gz<NGrid;gz++){
1070  for(int mx=0;mx<pNRow();mx++){
1071  int g1x = gx + mx - NMat2;
1072  if(g1x >= NGrid) g1x -= NGrid;
1073  if(g1x < 0) g1x += NGrid;
1074  for(int my=0;my<pNCol();my++){
1075  int g1y = gy + my - NMat2;
1076  if(g1y >= NGrid) g1y -= NGrid;
1077  if(g1y < 0) g1y += NGrid;
1078  for(int mz=0;mz<pNZed();mz++){
1079  int g1z = gz + mz - NMat2;
1080  if(g1z >= NGrid) g1z -= NGrid;
1081  if(g1z < 0) g1z += NGrid;
1082  Plot2[(gx*NGrid+gy)*NGrid+gz] += Plot[(g1x*NGrid+g1y)*NGrid+g1z]*Val(mx,my,mz);
1083  }
1084  }
1085  }
1086  }
1087  }
1088  }
1089  for(int gx=0;gx<NGrid;gx++){
1090  for(int gy=0;gy<NGrid;gy++){
1091  for(int gz=0;gz<NGrid;gz++){
1092  Plot[(gx*NGrid+gy)*NGrid+gz] = Plot2[(gx*NGrid+gy)*NGrid+gz];
1093  }
1094  }
1095  }
1096  free(Plot2);
1097 }
void FillGaussian5()
Fill the entries for the 5x5 Gauss blur.
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
int pNCol()
Number of columns.
void RandomFill(double Max)
Fill the entries randomly.
int Solve(double *Known, double *UnKnown)
Solve a system A|b = y.
int getNRow()
Return size.
void copymatrix(Matrice &source)
Copy matrix.
Matrice operator=(Matrice &)
Copy two matrices.
Geometrical operations on vectors.
Definition: MatematicaVect.h:9
void CopyOn(Matrice *B)
Copy on a matrix.
void FillDiffOperator(SPLINE Wg, int NDim)
Fill the entries of a differential operator.
void Mult(Matrice &A, Matrice &B)
Multiplication between two matrices.
Matrice operator*(Matrice &A)
Multiplies two matrices.
void ConvoluteMatrix2(double *Plot, int NGrid)
Convolute with a matrix.
double Det()
Computes the determinants.
double * x
Where the data are stored.
void ConvoluteMatrix(double *Plot, int NGrid, int NDim, int IfMinImConv)
Convolute with a matrix.
~Matrice()
Freeing.
double x
First basis component.
void Invert()
Invert.
void ConvoluteMatrix3(double *Plot, int NGrid)
Convolute with a matrix.
void Clear()
Set all the entries to zero.
void Normalize()
Normalize the matrix.
double a3
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double y
Second basis component.
void Transpose()
Transpose the matrix.
void ConvoluteMatrix1(double *Plot, int NGrid)
Convolute with a matrix.
void Multiply(double Val)
Multiply by a scalar.
double Val(int row)
Returns a value in 1d.
void ConvoluteMatrix2MinImConv(double *Plot, int NGrid)
Convolute with a matrix.
void setNRow(int newNRow)
Set new size.
void settoproduct(Matrice &left, Matrice &right)
Set to product.
bool Add(int row, int col, double Val)
Add the value of the coefficient to the previous one.
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double a1
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Matrice operator^(Matrice &) const
Tensor product?
Matrice computes the algebric operations on matrices.
bool Set(int row, int column, double Val)
Set a coefficient.
void Apply(double *Known, double *UnKnown)
Apply Ax = y.
int NDim
Dimension allocated.
double w
Forth basis component.
Quaternion class.
double z
Third basis component.
Matrice operator+(Matrice &)
Sums two matrices.
int pNZed()
Number of zed.
int pNRow()
Number of rows.
void getvalue(int row, int column, double &returnvalue, bool &success)
Return value.
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Coefficient of a spline.
void FillCanny()
Fill the entries for the Canny edge detector.
void Print()
Print the entries.
void comparetoidentity()
Compare to identity.
double * data
Stored entries.
void Shout(const char *s,...)
Name of the last function called. For debugging.
Matrice(int newNSize)
creates a square Matrice
void ConvoluteMatrix1MinImConv(double *Plot, int NGrid)
Convolute with a matrix.