VERB_code_2.3
Matrix.cpp
Go to the documentation of this file.
1 
6 #ifndef matrix_array_MATRIX_CPP
7 #define matrix_array_MATRIX_CPP
8 
9 #include "Matrix.h"
10 
11 using namespace std;
12 
13 #if defined(_WIN32) || defined(_WIN64)
14 #define strncasecmp _strnicmp
15 #define strcasecmp _stricmp
16 #endif
17 
18 // Memory related functions
19 
21 template<class T>inline T* matrix(long Rows)
22 {
23  T *m=new T[Rows];
24  // assert(m!=NULL);
25  if (m == NULL) {
26  throw error_msg("MEMORY_ERROR", "Memory can't be initialized: %s size", Rows*sizeof(T));
27  }
28  return m;
29 }
30 
32 template<class T>inline T** matrix(long Rows, long Columns)
33 {
34  // allocating memory for array of pinters
35  T **m=new T*[Rows];
36  // assert(m!=NULL);
37  if (m == NULL) {
38  throw error_msg("MEMORY_ERROR", "Memory can't be initialized: %s size", Rows * sizeof(T));
39  }
40  // allocating memory for data array
41  m[0] = new T[Rows * Columns];
42  // assert(m[0]!=NULL);
43  if (m[0] == NULL) {
44  throw error_msg("MEMORY_ERROR", "Memory can't be initialized: %s size", Rows * Columns * sizeof(T));
45  }
46  // assign pointers to data ranges
47  for(long i=1; i<Rows; i++) m[i] = m[i-1] + Columns;
48  return m;
49 }
50 
52 template<class T>inline T*** matrix(int size_x, int size_y, int size_z)
53 {
54  // allocating memory for array of pointers to pointers
55  T ***m=new T**[size_x];
56  // assert(m!=NULL);
57  if (m == NULL) {
58  throw error_msg("MEMORY_ERROR", "Memory can't be initialized: %s size", size_x * sizeof(T));
59  }
60  for (int x = 0; x < size_x; x++) {
61  // for each pointer allocating memory for array of pointers
62  m[x] = new T*[size_y];
63  // assert(m[x]!=NULL);
64  if (m[x] == NULL) {
65  throw error_msg("MEMORY_ERROR", "Memory can't be initialized: %s size", size_y * sizeof(T));
66  }
67  }
68  //allocating memory for data array
69  m[0][0] = new T[size_x * size_y * size_z]; // m[0][0] is pointer to array of T
70  // assert(m[0][0]!=NULL);
71  if (m[0][0] == NULL) {
72  throw error_msg("MEMORY_ERROR", "Memory can't be initialized: %s size", size_x * size_y * size_z * sizeof(T));
73  }
74  for (int x = 0; x < size_x; x++) {
75  for (int y = 0; y < size_y; y++) {
76  // assign pointers to data ranges
77  m[x][y] = m[0][0] + (x*size_y + y)*size_z;
78  }
79  }
80 
81  return m;
82 }
83 
86 template<class T>inline void free_matrix(T* m) {
87  delete m;
88 }
89 
92 template<class T>inline void free_matrix(T** m) {
93  delete[](m[0]);
94  delete[](m);
95 }
96 
99 template<class T>inline void free_matrix(T*** m, int size_x, int size_y) {
100  delete[](m[0][0]);
101  for (int x = 0; x < size_x; x++) {
102  delete[](m[x]);
103  }
104  delete[](m);
105 }
106 
108 //
109 // Matrix 1D
110 //
112 
122 template<class T>
123 Matrix1D<T>::Matrix1D( int x_size , string name) {
124  initialized = false;
125  this->name = name;
126  // allocating memory
127  AllocateMemory(x_size);
128 }
129 
138 template<class T>
139 Matrix1D<T>::Matrix1D( int x_size ) {
140  initialized = false;
141  // allocating memory
142  AllocateMemory(x_size);
143 }
144 
153 template<class T>
155  initialized = false;
156  this->operator = (M);
157 }
158 
162 template<class T>
164  if (initialized) free_matrix<T>(matrix_array);
165 }
166 
172 template<class T>
173 void Matrix1D<T>::AllocateMemory( int x_size ) {
174  this->size_x = x_size;
175  // using inline template for memory allocation
176  matrix_array = matrix<T>(x_size);
177  initialized = true;
178  // for (int i = 0; i < size_x; i++) matrix_array[i] = 2;
179 }
180 
181 
188 template<class T>
189 inline T& Matrix1D<T>::operator[](int i) {
190 #ifdef DEBUG_MODE
191  if (!initialized) {
192  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
193  }
194 #endif
195  // return i-value
196  return matrix_array[i];
197 }
198 
205 template<class T>
206 inline T& Matrix1D<T>::operator[](int i) const {
207 #ifdef DEBUG_MODE
208  if (!initialized) {
209  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
210  }
211 #endif
212  // return i-value
213  return matrix_array[i];
214 }
215 
223 template<class T>
224 inline T& Matrix1D<T>::operator()(int x) {
225 #ifdef DEBUG_MODE
226  if (!initialized) {
227  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
228  }
229 #endif
230  // return x-value
231  return matrix_array[x];
232 }
233 
240 template<class T>
242  // Check, if both matrix aren't the same one
243  if (this!=&M) {
244  // check, if RHS matrix was initialized
245  if (M.initialized) {
246  // free LHS matrix, if it was initialized
247  if (initialized) free_matrix<T>(matrix_array);
248  // copy parameters
249  this->size_x = M.size_x;
250  this->name = M.name;
251  // allocating memory for LHS matrix
252  this->AllocateMemory(this->size_x);
253  // fast values copying as a memory range
254  memcpy( matrix_array, M.matrix_array, this->size_x * sizeof( T ) );
255  } else {
256  this->initialized = false;
257  }
258  }
259  return *this;
260 }
261 
268 template<class T>
270  for (int x = 0; x < this->size_x; x++)
271  matrix_array[x] = Val;
272  return *this;
273 }
274 
281 template<class T>
282 inline Matrix1D<T> Matrix1D<T>::operator* (const T Val) const {
283  Matrix1D<T> Tmp(*this);
284  for (int x = 0; x < this->size_x; x++)
285  Tmp[x] = matrix_array[x] * Val;
286  return Tmp;
287 }
288 
295 template<class T>
296 inline Matrix1D<T> Matrix1D<T>::operator/ (const T Val) const {
297  int x;
298  Matrix1D<T> Tmp(*this);
299  for (x = 0; x < this->size_x; x++)
300  Tmp[x] = matrix_array[x] / Val;
301  return Tmp;
302 }
303 
309 template<class T>
310 inline Matrix1D<T> Matrix1D<T>::times (const Matrix1D<T> &M) const {
311  int x;
312  Matrix1D<T> Tmp(*this);
313  for (x = 0; x < this->size_x; x++)
314  Tmp[x] = matrix_array[x] * M.matrix_array[x];
315  return Tmp;
316 }
317 
323 template<class T>
324 inline Matrix1D<T> Matrix1D<T>::divide (const Matrix1D<T> &M) const {
325  int x;
326  Matrix1D<T> Tmp(*this);
327  for (x = 0; x < this->size_x; x++)
328  Tmp[x] = matrix_array[x] / M.matrix_array[x];
329  return Tmp;
330 }
331 
335 template<class T>
336 inline T Matrix1D<T>::norm() const {
337  T res = 0;
338  int x;
339  for (x = 0; x < this->size_x; x++) {
340  res += matrix_array[x] * matrix_array[x];
341  }
342  return sqrt(res);
343 }
344 
348 template<class T>
349 inline T Matrix1D<T>::dot( const Matrix1D<T> &W ) const {
350  if (this->size_x != W.size_x)
351  throw error_msg("DOT_PRODUCT", "Size is different");
352 
353  T res = 0;
354  int x;
355  for (x = 0; x < this->size_x; x++) {
356  res += matrix_array[x] * W[x];
357  }
358  return res;
359 }
360 
361 
365 template<class T>
366 void Matrix1D<T>::writeToFile(string filename) {
367  int x;
368  ofstream output(filename.c_str());
369  if (output==NULL && (filename.find("Debug") == string::npos)) {
370  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
371  }
372  output << "VARIABLES = \"" << ((this->name!="")?this->name:"function") << "\" "<< endl;
373  output << "ZONE T=\"" << filename << "\", I=" << size_x << endl;
374  output.setf(ios_base::scientific, ios_base::floatfield);
375  for (x = 0; x < size_x; x++) {
376  output << matrix_array[x] << endl;
377  }
378  output.close();
379 }
380 
384 template<class T>
385 void Matrix1D<T>::writeToFile(string filename, Matrix1D<T> &grid_x) {
386  int x;
387  ofstream output(filename.c_str());
388  if (output==NULL && (filename.find("Debug") == string::npos)) {
389  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
390  }
391  output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((this->name!="")?this->name:"function") << "\" "<< endl;
392  output << "ZONE T=\"" << filename << "\", I=" << size_x << endl;
393  output.setf(ios_base::scientific, ios_base::floatfield);
394  for (x = 0; x < size_x; x++) {
395  output << "\t" << grid_x[x] << "\t" << "\t" << matrix_array[x] << endl;
396  }
397  output.close();
398 }
399 
400 
404 template<class T>
405 void Matrix1D<T>::readFromFile(string filename) {
406  int x;
407  string inBuf;
408  if (!initialized) {
409  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
410  } else {
411  ifstream input(filename.c_str());
412  if (input != NULL && !input.eof()) {
413  // Skipping first two lines.
414  input >> inBuf;
415  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
416  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
417  input >> inBuf;
418  }
419  // read to the end of the line with 'zone'
420  input.ignore(9999, '\n');
421  for (x = 0; x < size_x; x++) {
422  input >> matrix_array[x];
423  }
424  } else {
425  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
426  //cout << "Error reading initial conditions input file " << filename << endl;
427  //getchar();
428  //exit(-1);
429  }
430  input.close();
431  }
432 }
433 
437 template<class T>
438 void Matrix1D<T>::readFromFile(string filename, Matrix1D<T> &grid_x) {
439  int x;
440  string inBuf;
441  double loaded_x;
442  double err = 1e-8;
443  if (!initialized) {
444  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
445  } else {
446  ifstream input(filename.c_str());
447  if (input != NULL && !input.eof()) {
448  // Skipping first two lines. Should serch for 'ZONE' and read from following line better - done.
450  input >> inBuf;
451  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
452  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
453  input >> inBuf;
454  }
455  // read to the end of the line with 'zone'
456  input.ignore(9999, '\n');
457  for (x = 0; x < size_x; x++) {
458  input >> loaded_x;
459  // check if grid is the same
460  if (fabs(loaded_x - grid_x[x]) > err) {
461  throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e\nGrid: %e\n", filename.c_str(), loaded_x, grid_x[x]);
462  } else {
463  input >> matrix_array[x];
464  }
465  }
466  } else {
467  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
468  //cout << "Error reading initial conditions input file " << filename << endl;
469  //getchar();
470  //exit(-1);
471  }
472  input.close();
473  }
474 }
475 
476 
478 //
479 // Matrix 2D
480 //
482 
490 template<class T>
491 Matrix2D<T>::Matrix2D( int x_size, int y_size ) {
492  initialized = false;
493  // allocating memory
494  AllocateMemory(x_size, y_size);
495 }
496 
503 template<class T>
505  initialized = false;
506  this->operator = (M);
507 }
508 
512 template<class T>
514  if (initialized) free_matrix<T>(matrix_array);
515 }
516 
523 template<class T>
524 void Matrix2D<T>::AllocateMemory( int x_size, int y_size ) {
525  this->size_x = x_size;
526  this->size_y = y_size;
527  // using matrix inline template to allocate memory
528  matrix_array = matrix<T>(x_size, y_size);
529  initialized = true;
530  // for (int i = 0; i < size_x; i++)
531  // for (int j = 0; j < size_y; j++)
532  // matrix_array[i][j] = 0;
533 }
534 
540 template<class T>
542  // check, if both matrix are the same one
543  if (this!=&M) {
544  // check, if RHS matrix in initialized
545  if (M.initialized) {
546  // free LHS matrix if it was initialized
547  if (initialized && (size_x != M.size_x || size_y != M.size_y)) {
548  free_matrix<T>(matrix_array);
549  initialized = false;
550  }
551  if (!initialized) {
552  this->size_x = M.size_x;
553  this->size_y = M.size_y;
554  // allocating memory for LHS matrix (and creating the correct pointers M[0] - M[N] to the matrix rows)
555  this->AllocateMemory(this->size_x, this->size_y);
556  }
557  this->name = M.name;
558  // for (int x = 0; x < this->size_x; x++)
559  // memcpy( matrix_array[x], M.matrix_array[x], this->size_y * sizeof( T ) );
560  // fast values copying as a memory range
561  memcpy(matrix_array[0], M.matrix_array[0], this->size_x * this->size_y * sizeof(T) );
562  } else {
563  this->initialized = false;
564  }
565  }
566  return *this;
567 }
568 
574 template<class T>
576  int x, y;
577  if (initialized) {
578  for (x = 0; x < this->size_x; x++)
579  for (y = 0; y < this->size_y; y++)
580  matrix_array[x][y] = val;
581  } else {
582  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
583  }
584 
585  return *this;
586 }
587 
593 template<class T>
594 inline Matrix2D<T> Matrix2D<T>::operator* (const T Val) const {
595  int x, y;
596  Matrix2D<T> Tmp(*this);
597  for (x = 0; x < this->size_x; x++)
598  for (y = 0; y < this->size_y; y++)
599  Tmp[x][y] = matrix_array[x][y] * Val;
600  return Tmp;
601 }
602 
608 template<class T>
609 inline Matrix2D<T> Matrix2D<T>::operator/(const T Val) const {
610  int x, y;
611  Matrix2D<T> Tmp(*this);
612  for (x = 0; x < this->size_x; x++)
613  for (y = 0; y < this->size_y; y++)
614  Tmp[x][y] = matrix_array[x][y] / Val;
615  return Tmp;
616 }
617 
623 template<class T>
624 inline Matrix2D<T> Matrix2D<T>::divide (const Matrix2D<T> &M) const {
625  int x, y;
626  Matrix2D<T> Tmp(*this);
627  for (x = 0; x < this->size_x; x++)
628  for (y = 0; y < this->size_y; y++)
629  Tmp[x][y] = matrix_array[x][y] / M[x][y];
630  return Tmp;
631 }
632 
638 template<class T>
639 inline Matrix2D<T> Matrix2D<T>::times (const Matrix2D<T> &M) const {
640  int x, y;
641  Matrix2D<T> Tmp(*this);
642  for (x = 0; x < this->size_x; x++)
643  for (y = 0; y < this->size_y; y++)
644  Tmp[x][y] = matrix_array[x][y] * M[x][y];
645  return Tmp;
646 }
647 
653 template<class T>
655  int x, y;
656  Matrix2D<T> Tmp(*this);
657  for (x = 0; x < this->size_x; x++)
658  for (y = 0; y < this->size_y; y++)
659  Tmp[x][y] = (matrix_array[x][y]>val)?matrix_array[x][y]:val;
660  return Tmp;
661 }
662 
666 template<class T>
667 inline int Matrix2D<T>::index1d(int x, int y) const {
668  return x*size_y + y;
669 }
670 
671 
672 
681 template<class T>
682 void Matrix2D<T>::writeToFile(string filename) {
683  int x, y;
684  ofstream output(filename.c_str());
685  if (output==NULL && (filename.find("Debug") == string::npos)) {
686  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
687  }
688  output << "VARIABLES = \""<< ((this->name!="")?this->name:"f") << "\" "<< endl;
689  output << "ZONE T=\"" << filename << "\", I=" << size_y << ", J= " << size_x << endl;
690  output.setf(ios_base::scientific, ios_base::floatfield);
691  for (x = 0; x < size_x; x++) {
692  for (y = 0; y < size_y; y++) {
693  output << matrix_array[x][y] << endl;
694  }
695  }
696  output.close();
697 }
698 
699 
704 template<class T>
705 void Matrix2D<T>::writeToFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y) {
706  int x, y;
707  ofstream output(filename.c_str());
708  output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((grid_y.name!="")?grid_y.name:"y") << "\", \"" << ((this->name!="")?this->name:"f") << "\" "<< endl;
709  output << "ZONE T=\"" << filename << "\", I=" << size_y << ", J=" << size_x << endl;
710  output.setf(ios_base::scientific, ios_base::floatfield);
711  for (x = 0; x < size_x; x++) {
712  for (y = 0; y < size_y; y++) {
713  output << "\t" << grid_x[x][y] << "\t" << grid_y[x][y] << "\t" << matrix_array[x][y] << endl;
714  }
715  }
716  output.close();
717 }
718 
719 
723 template<class T>
724 void Matrix2D<T>::readFromFile(string filename) {
725  int x, y;
726  string inBuf;
727  if (!initialized) {
728  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
729  } else {
730  ifstream input(filename.c_str());
731  if (input != NULL && !input.eof()) {
732  // Skipping first two lines.
733  input >> inBuf;
734  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
735  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
736  input >> inBuf;
737  }
738  // read to the end of the line with 'zone'
739  input.ignore(9999, '\n');
740  for (x = 0; x < size_x; x++) {
741  for (y = 0; y < size_y; y++) {
742  input >> matrix_array[x][y];
743  }
744  }
745  } else {
746  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
747  //cout << "Error reading initial conditions input file " << filename << endl;
748  //getchar();
749  //exit(-1);
750  }
751  input.close();
752  }
753 }
754 
758 template<class T>
759 void Matrix2D<T>::readFromFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y) {
760  int x, y;
761  string inBuf;
762  double err = 1e-6;
763  double loaded_x, loaded_y;
764  if (!initialized) {
765  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
766  } else {
767  ifstream input(filename.c_str());
768  if (input != NULL && !input.eof()) {
769  // Skipping first two lines.
770  input >> inBuf;
771  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
772  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
773  input >> inBuf;
774  }
775  // read to the end of the line with 'zone'
776  input.ignore(9999, '\n');
777  for (x = 0; x < size_x; x++) {
778  for (y = 0; y < size_y; y++) {
779  input >> loaded_x >> loaded_y;
780  // check if grid is the same
781  if (fabs(loaded_x - grid_x[x][y]) > err || fabs(loaded_y - grid_y[x][y]) > err) {
782  throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e, %e\nGrid: %e, %e\n", filename.c_str(), loaded_x, loaded_y, grid_x[x][y], grid_y[x][y]);
783  } else {
784  input >> matrix_array[x][y];
785  }
786  }
787  }
788  } else {
789  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
790  //cout << "Error reading initial conditions input file " << filename << endl;
791  //getchar();
792  //exit(-1);
793  }
794  input.close();
795  }
796 }
797 
798 
800 //
801 // Matrix 3D
802 //
804 
809 template<class T>
810 Matrix3D<T>::Matrix3D( int x_size, int y_size, int z_size ) {
811  initialized = false;
812  // allocating memory
813  AllocateMemory(x_size, y_size, z_size);
814 }
815 
822 template<class T>
824  initialized = false;
825  this->operator = (M);
826 }
827 
831 template<class T>
833  if (initialized) free_matrix<T>(matrix_array, size_x, size_y);
834  // delete[] plane_array; - it's get deleted with matrix_array[0][0][0] or something?
835 }
836 
840 template<class T>
841 void Matrix3D<T>::AllocateMemory( int x_size, int y_size, int z_size) {
842  this->size_x = x_size;
843  this->size_y = y_size;
844  this->size_z = z_size;
845  matrix_array = matrix<T>(x_size, y_size, z_size);
846  plane_array = matrix_array[0][0];
847  initialized = true;
848 #ifdef DEBUG_MODE
849  // should not initialize matrix with zeros, it can slow the code greatly in some cases
850  for (int i = 0; i < size_x; i++)
851  for (int j = 0; j < size_y; j++)
852  for (int k = 0; k < size_z; k++)
853  matrix_array[i][j][k] = 0;
854 #endif
855 }
856 
863 template<class T>
864 inline T** Matrix3D<T>::operator[] (int i) {
865 #ifdef DEBUG_MODE
866  if (!initialized) {
867  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
868  }
869 #endif
870  // return i-th pointer
871  return matrix_array[i];
872 }
873 
879 template<class T>
880 inline T& Matrix3D<T>::operator() (int x, int y, int z) {
881 #ifdef DEBUG_MODE
882  if (!initialized) {
883  throw error_msg("MATRIX_ERROR", "Using not initialized matrix");
884  }
885  if ((x < 0 || x > size_x-1) || (y < 0 || y > size_y-1) || (z < 0 || z > size_z-1)) {
886  throw error_msg("MATRIX_ERROR", "Index is out of bound");
887  }
888 #endif
889  // return (x,y,z) value
890  return plane_array[(x*size_y + y)*size_z + z];
891 }
892 
898 template<class T>
900  // check if not LHS and RHS matrix are the same
901  if (this!=&M) {
902  // check if RHS matrix is initialized
903  if (M.initialized) {
904  // free LHS matrix, if it is initialized
905  if (initialized && (size_x != M.size_x || size_y != M.size_y || size_z != M.size_z)) {
906  free_matrix<T>(matrix_array, this->size_x, this->size_y);
907  initialized = false;
908  }
909  if (!initialized) {
910  this->size_x = M.size_x;
911  this->size_y = M.size_y;
912  this->size_z = M.size_z;
913  // allocating memory for LHS matrix
914  this->AllocateMemory(this->size_x, this->size_y, this->size_z);
915  }
916  this->name = M.name;
917  // for (int x = 0; x < this->size_x; x++)
918  // for (int y = 0; y < this->size_y; y++)
919  // memcpy( matrix_array[x][y], M.matrix_array[x][y], this->size_z * sizeof( T ) );
920  // fast values copying as memory range
921  memcpy( matrix_array[0][0], M.matrix_array[0][0], this->size_x * this->size_y * this->size_z * sizeof( T ) );
922  } else {
923  this->initialized = false;
924  }
925  }
926  return *this;
927 }
928 
933 /* template<class T>
934 Matrix3D<T>& Matrix3D<T>::operator= (const Matrix2D<T> &M) {
935 double val;
936 if (M.initialized) {
937 if (initialized) free_matrix<T>(matrix_array, this->size_x, this->size_y);
938 this->size_x = 1;
939 this->size_y = M.size_x;
940 this->size_z = M.size_y;
941 //this->name = M.name;
942 this->AllocateMemory(this->size_x, this->size_y, this->size_z);
943 for (int y = 0; y < this->size_y; y++)
944 //memcpy( matrix_array[0][y], M.matrix_array[y], this->size_z * sizeof( T ) );
945 for (int z = 0; z < this->size_z; z++) {
946 this->matrix_array[0][y][z] = M[y][z];
947 }
948 } else {
949 this->initialized = false;
950 }
951 return *this;
952 }*/
953 
957 template<class T>
959  // std::fill_n(matrix_array[0][0][0], size_x*size*y*size*z, Val); // might be faster
960  for (int x = 0; x < this->size_x; x++)
961  for (int y = 0; y < this->size_y; y++)
962  for (int z = 0; z < this->size_z; z++)
963  this->matrix_array[x][y][z] = Val;
964  return *this;
965 }
966 
967 
971 template<class T>
973  int x, y, z;
974  for (x = 0; x < this->size_x; x++)
975  for (y = 0; y < this->size_y; y++)
976  for (z = 0; z < this->size_z; z++)
977  matrix_array[x][y][z] += M.matrix_array[x][y][z];
978  return *this;
979 }
980 
984 template<class T>
986  int x, y, z;
987  for (x = 0; x < this->size_x; x++)
988  for (y = 0; y < this->size_y; y++)
989  for (z = 0; z < this->size_z; z++)
990  matrix_array[x][y][z] -= M.matrix_array[x][y][z];
991  return *this;
992 }
993 
997 template<class T>
999  int x, y, z;
1000  for (x = 0; x < this->size_x; x++)
1001  for (y = 0; y < this->size_y; y++)
1002  for (z = 0; z < this->size_z; z++)
1003  matrix_array[x][y][z] *= Val;
1004  return *this;
1005 }
1006 
1010 template<class T>
1012  int x, y, z;
1013  for (x = 0; x < this->size_x; x++)
1014  for (y = 0; y < this->size_y; y++)
1015  for (z = 0; z < this->size_z; z++)
1016  matrix_array[x][y][z] /= Val;
1017  return *this;
1018 }
1019 
1023 template<class T>
1025  int x, y, z;
1026  for (x = 0; x < this->size_x; x++)
1027  for (y = 0; y < this->size_y; y++)
1028  for (z = 0; z < this->size_z; z++)
1029  matrix_array[x][y][z] += Val;
1030  return *this;
1031 }
1032 
1036 template<class T>
1038  int x, y, z;
1039  for (x = 0; x < this->size_x; x++)
1040  for (y = 0; y < this->size_y; y++)
1041  for (z = 0; z < this->size_z; z++)
1042  matrix_array[x][y][z] -= Val;
1043  return *this;
1044 }
1045 
1049 template<class T>
1051  int x, y, z;
1052  for (x = 0; x < this->size_x; x++)
1053  for (y = 0; y < this->size_y; y++)
1054  for (z = 0; z < this->size_z; z++)
1055  matrix_array[x][y][z] *= M.matrix_array[x][y][z];
1056  return *this;
1057 }
1058 
1062 template<class T>
1064  int x, y, z;
1065  for (x = 0; x < this->size_x; x++)
1066  for (y = 0; y < this->size_y; y++)
1067  for (z = 0; z < this->size_z; z++)
1068  matrix_array[x][y][z] /= M.matrix_array[x][y][z];
1069  return *this;
1070 }
1071 
1075 template<class T>
1077  int x, y, z;
1078  Matrix3D<T> Tmp(*this);
1079  for (x = 0; x < this->size_x; x++)
1080  for (y = 0; y < this->size_y; y++)
1081  for (z = 0; z < this->size_z; z++)
1082  Tmp[x][y][z] = matrix_array[x][y][z] + M.matrix_array[x][y][z];
1083  return Tmp;
1084 }
1085 
1089 template<class T>
1091  int x, y, z;
1092  Matrix3D<T> Tmp(*this);
1093  for (x = 0; x < this->size_x; x++)
1094  for (y = 0; y < this->size_y; y++)
1095  for (z = 0; z < this->size_z; z++)
1096  Tmp[x][y][z] = matrix_array[x][y][z] - M.matrix_array[x][y][z];
1097  return Tmp;
1098 }
1099 
1100 
1104 template<class T>
1105 inline Matrix3D<T> Matrix3D<T>::operator* (const T Val) const {
1106  int x, y, z;
1107  Matrix3D<T> Tmp(*this);
1108  for (x = 0; x < this->size_x; x++)
1109  for (y = 0; y < this->size_y; y++)
1110  for (z = 0; z < this->size_z; z++)
1111  Tmp[x][y][z] = matrix_array[x][y][z] * Val;
1112  return Tmp;
1113 }
1114 
1115 
1119 template<class T>
1120 inline Matrix3D<T> Matrix3D<T>::operator/ (const T Val) const {
1121  int x, y, z;
1122  Matrix3D<T> Tmp(*this);
1123  for (x = 0; x < this->size_x; x++)
1124  for (y = 0; y < this->size_y; y++)
1125  for (z = 0; z < this->size_z; z++)
1126  Tmp[x][y][z] = matrix_array[x][y][z] / Val;
1127  return Tmp;
1128 }
1129 
1133 template<class T>
1134 inline Matrix3D<T> Matrix3D<T>::times (const Matrix3D<T> &M) const {
1135  int x, y, z;
1136  Matrix3D<T> Tmp(*this);
1137  for (x = 0; x < this->size_x; x++)
1138  for (y = 0; y < this->size_y; y++)
1139  for (z = 0; z < this->size_z; z++)
1140  Tmp[x][y][z] = matrix_array[x][y][z] * M.matrix_array[x][y][z];
1141  return Tmp;
1142 }
1143 
1147 template<class T>
1149  int x, y, z;
1150  Matrix3D<T> Tmp(*this);
1151  for (x = 0; x < this->size_x; x++)
1152  for (y = 0; y < this->size_y; y++)
1153  for (z = 0; z < this->size_z; z++)
1154  Tmp[x][y][z] = matrix_array[x][y][z] / M.matrix_array[x][y][z];
1155  return Tmp;
1156 }
1157 
1158 
1163 template<class T>
1164 void Matrix3D<T>::writeToFile(string filename) {
1165  int x, y, z;
1166  ofstream output(filename.c_str());
1167  output << "VARIABLES = \""<< ((this->name!="")?this->name:"f") <<"\" "<< endl;
1168  output << "ZONE T=\"" << filename << "\", I=" << size_z << ", J=" << size_y << ", K=" << size_x << endl;
1169  output.setf(ios_base::scientific, ios_base::floatfield);
1170  for (x = 0; x < size_x; x++) {
1171  for (y = 0; y < size_y; y++) {
1172  for (z = 0; z < size_z; z++) {
1173  output << matrix_array[x][y][z] << endl;
1174  }
1175  }
1176  }
1177  output.close();
1178 }
1179 
1184 template<class T>
1185 void Matrix3D<T>::writeToFile(string filename, Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z) {
1186  int x, y, z;
1187  ofstream output(filename.c_str());
1188  if (output==NULL && (filename.find("Debug") == string::npos)) {
1189  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
1190  }
1191  output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((grid_y.name!="")?grid_y.name:"y") << "\", \"" << ((grid_z.name!="")?grid_z.name:"z") << "\", \"" << ((this->name!="")?this->name:"f") << "\" "<< endl;
1192  output << "ZONE T=\"" << filename << "\", I=" << size_z << ", J=" << size_y << ", K=" << size_x << endl;
1193  output.setf(ios_base::scientific, ios_base::floatfield);
1194  for (x = 0; x < size_x; x++) {
1195  for (y = 0; y < size_y; y++) {
1196  for (z = 0; z < size_z; z++) {
1197  output << "\t" << grid_x[x][y][z] << "\t" << grid_y[x][y][z] << "\t" << grid_z[x][y][z] << "\t" << matrix_array[x][y][z] << endl;
1198  }
1199  }
1200  }
1201  output.close();
1202 }
1203 
1207 template<class T>
1208 void Matrix3D<T>::readFromFile(string filename) {
1209  int x, y, z;
1210  string inBuf;
1211  if (!initialized) {
1212  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
1213  } else {
1214  ifstream input(filename.c_str());
1215  if (input != NULL && !input.eof()) {
1216  // Skipping first two lines.
1217  input >> inBuf;
1218  // !!! while (strcmp(strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
1219  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
1220 
1221  input >> inBuf;
1222  }
1223  // read to the end of the line with 'zone'
1224  input.ignore(9999, '\n');
1225  for (x = 0; x < size_x; x++) {
1226  for (y = 0; y < size_y; y++) {
1227  for (z = 0; z < size_z; z++) {
1228  input >> matrix_array[x][y][z];
1229  }
1230  }
1231  }
1232  } else {
1233  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
1234  //cout << "Error reading initial conditions input file " << filename << endl;
1235  //getchar();
1236  //exit(-1);
1237  }
1238  input.close();
1239  }
1240 }
1241 
1245 template<class T>
1246 void Matrix3D<T>::readFromFile(string filename, Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z) {
1247  int x, y, z;
1248  string inBuf;
1249  if (!initialized) {
1250  throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
1251  } else {
1252  ifstream input(filename.c_str());
1253  if (input != NULL && !input.eof()) {
1254  // Skipping first two lines.
1255  input >> inBuf;
1256  // !!! while (strcmp(_strupr((char *) inBuf.c_str()), "ZONE") != 0 && !input.eof() ) {
1257  // !!! while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
1258  while (strcmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
1259  input >> inBuf;
1260  }
1261  // read to the end of the line with 'zone'
1262  input.ignore(9999, '\n');
1263  for (x = 0; x < size_x; x++) {
1264  for (y = 0; y < size_y; y++) {
1265  for (z = 0; z < size_z; z++) {
1266  input >> grid_x[x][y][z] >> grid_y[x][y][z] >> grid_z[x][y][z] >> matrix_array[x][y][z];
1267  // skip till the end of the line
1268  input.ignore(9999, '\n');
1269 
1270  // input >> loaded_x >> loaded_y >> loaded_z;
1271  // check if grid is the same
1272  // if (fabs(loaded_x - grid_x[x][y][z]) > err || fabs(loaded_y - grid_y[x][y][z]) > err || fabs(loaded_z - grid_z[x][y][z]) > err) {
1273  // throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e, %e, %e\nGrid: %e, %e, %e\n", filename.c_str(), loaded_x, loaded_y, loaded_z, grid_x[x][y][z], grid_y[x][y][z], grid_z[x][y][z]);
1274  //} else {
1275  // input >> matrix_array[x][y][z];
1276  //}
1277  }
1278  }
1279  }
1280  } else {
1281  throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
1282  //cout << "Error reading initial conditions input file " << filename << endl;
1283  //getchar();
1284  //exit(-1);
1285  }
1286  input.close();
1287  }
1288 }
1289 
1290 
1291 
1295 template<class T>
1296 inline int Matrix3D<T>::index1d(int x, int y, int z) {
1297  return (x*size_y + y)*size_z + z;
1298 }
1299 
1300 
1304 template<class T>
1306  T tmp = 0;
1307  int x, y, z;
1308  for (x = 0; x < size_x; x++) {
1309  for (y = 0; y < size_y; y++) {
1310  for (z = 0; z < size_z; z++) {
1311  tmp = (tmp>matrix_array[x][y][z])?tmp:matrix_array[x][y][z];
1312  }
1313  }
1314  }
1315  return tmp;
1316 }
1317 
1321 template<class T>
1323  T tmp = 0;
1324  int x, y, z;
1325  for (x = 0; x < size_x; x++) {
1326  for (y = 0; y < size_y; y++) {
1327  for (z = 0; z < size_z; z++) {
1328  tmp = (tmp>fabs((double)matrix_array[x][y][z]))?tmp:fabs((double)matrix_array[x][y][z]);
1329  }
1330  }
1331  }
1332  return tmp;
1333 }
1334 
1338 template<class T>
1340  Matrix3D<T> tmp(this->size_x, this->size_y, this->size_z);
1341  int x, y, z;
1342  for (x = 0; x < size_x; x++) {
1343  for (y = 0; y < size_y; y++) {
1344  for (z = 0; z < size_z; z++) {
1345  tmp[x][y][z] = (matrix_array[x][y][z]>0)?matrix_array[x][y][z]:-matrix_array[x][y][z];
1346  }
1347  }
1348  }
1349  return tmp;
1350 }
1351 
1352 
1356 template<class T>
1358  int y, z;
1359  Matrix2D<T> tmp(this->size_y, this->size_z);
1360  tmp.name = this->name + "_slice";
1361  for (y = 0; y < size_y; y++) {
1362  for (z = 0; z < size_z; z++) {
1363  tmp[y][z] = matrix_array[p_x][y][z];
1364  }
1365  }
1366  return tmp;
1367 }
1368 
1372 template<class T>
1374  int x, z;
1375  Matrix2D<T> tmp(this->size_x, this->size_z);
1376  tmp.name = this->name + "_slice";
1377  for (x = 0; x < size_x; x++) {
1378  for (z = 0; z < size_z; z++) {
1379  tmp[x][z] = matrix_array[x][p_y][z];
1380  }
1381  }
1382  return tmp;
1383 }
1384 
1388 template<class T>
1390  int x, y;
1391  Matrix2D<T> tmp(this->size_x, this->size_y);
1392  tmp.name = this->name + "_slice";
1393  for (x = 0; x < size_x; x++) {
1394  for (y = 0; y < size_y; y++) {
1395  tmp[x][y] = matrix_array[x][y][p_z];
1396  }
1397  }
1398  return tmp;
1399 }
1400 
1401 
1403 //
1404 // Interpolations (mess)
1405 //
1407 
1408 
1413 template<class T>
1414 void Matrix1D<T>::Interpolate(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1415 
1416 
1417  //if (size_x != old_grid.size_x) {
1418  if (this->size_x != new_grid.size_x) {
1419  throw error_msg("INTERPOLATION_ERROR", "New grid and function size mismatch!\n");
1420  //cout << "new grid and function size mismatch!" << endl;
1421  //return;
1422  }
1423 
1424  //Matrix1D<T> res(new_grid.size_x);
1425  //Maths::Interpolation::Linear LinInterp(old_grid.size_x, old_grid.matrix_array, matrix_array);
1426  // initialize some class. It was like that when I downloaded the function.
1427  Maths::Interpolation::Linear LinInterp(old_grid.size_x, old_grid.matrix_array, old_function.matrix_array);
1428 
1429  int i;
1430  for (i = 0; i < new_grid.size_x; i++)
1431  //res[i] = LinInterp.getValue(new_grid.matrix_array[i]);
1432  matrix_array[i] = LinInterp.getValue(new_grid.matrix_array[i]);
1433 
1434  //return res;
1435 
1436 }
1437 
1454 template<class T>
1455 void Matrix1D<T>::Spline(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub, double lin_spline_coef, double max_second_der) {
1456 
1457  //if (size_x != old_grid.size_x) {
1458  // check, if size of the grids are equal
1459  if (this->size_x != new_grid.size_x) {
1460  throw error_msg("SPLINE_INTERPOLATION", "New grid and function size mismatch!");
1461  }
1462  // if size == 1 - just copy the value
1463  if (this->size_x == 1) {
1464  this->matrix_array[0] = old_function.matrix_array[0];
1465  return;
1466  }
1467 
1468  // Idea was to find the intersection between grids, to interpolate insite the intersection and do not extrapolate outside. But is done later.
1469  int i=0, i_min=0, i_max=old_grid.size_x-1, interp_size=0;
1470  // for(i = 0; old_grid.matrix_array[i] <= new_grid.matrix_array[0] && i < old_grid.size_x-1; i++) {}
1471  // i_min = i;
1472  // for(i = old_grid.size_x-1; old_grid.matrix_array[i] >= new_grid.matrix_array[old_grid.size_x-1] && i > 0; i--) {}
1473  // i_max = i;
1474  //cout << i_min << " " << i_max << endl;
1475 
1476  // new range for spline interpolation - values close to the boundaries are interpolated by linear interpolation
1477  i_min = 1; i_max = old_grid.size_x-2;
1478  interp_size = i_max-i_min+1;
1479 
1480  double yp0, ypn;
1481  //yp0 = 2*(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]) - (old_function[2] - old_function[1])/(old_grid[2] - old_grid[1]);
1482  //ypn = 2*(old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]) - (old_function[size_x-2] - old_function[size_x-3])/(old_grid[size_x-2] - old_grid[size_x-3]);
1483  // ypn = (old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1484 
1485  // flag, means use smooth second derivative on the boundary
1486  yp0 = 1e99;
1487  //yp0 = -(old_function[2] - old_function[0])/(old_grid[2] - old_grid[0]);
1488  ypn = 1e99;//0;
1489 
1490  Matrix1D<T> y2(old_grid.size_x);
1491  //spline(old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, yp0, ypn, y2.matrix_array);
1492  //spline(&old_grid.matrix_array[1], &old_function.matrix_array[1], old_grid.size_x-2, yp0, ypn, &y2.matrix_array[1]);
1493  // calculating array of second derivatives
1494  spline(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], interp_size, yp0, ypn, &y2.matrix_array[i_min]);
1495 
1496  // make it more linear (lower second derivatives)
1497  for (i = 0; i < old_grid.size_x; i++) {
1498  // if second derivative value is too big
1499  if (fabs(y2[i]) >= max_second_der) {
1500  if (y2[i] > 0) {
1501  // decreasing by lin_spline_coef big second derivative values - makes interpolation more linear
1502  // see numerical recipes
1503  y2[i] = max_second_der + (y2[i] - max_second_der)*(1.0 - lin_spline_coef);
1504  } else {
1505  y2[i] = - max_second_der + (y2[i] + max_second_der)*(1.0 - lin_spline_coef);
1506  }
1507  }
1508  }
1509 
1510  // for all grid points - looking for new values
1511  double a; int k; // interpolation parameters
1512  for (i = 0; i < new_grid.size_x; i++) {
1513  if (i == 0) { // if we are on the boundary
1514  this->matrix_array[i] = lb;
1515  } else if (i == new_grid.size_x-1) { // if we are on the boundary
1516  this->matrix_array[i] = ub;
1517  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) { // if it is extrapolation
1518  this->matrix_array[i] = lb;
1519  //this->matrix_array[i] = old_function.matrix_array[0];
1520  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) { // if it is extrapolation
1521  this->matrix_array[i] = ub;
1522  //this->matrix_array[i] = old_function.matrix_array[old_function.size_x-1];
1523  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) { // if the new grid point is between two boundary points of the old grid
1524  // linear
1525  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
1526  this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
1527  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) { // if the new grid point is between two boundary points of the old grid
1528  // linear
1529  //double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-2])/(old_grid.matrix_array[old_grid.size_x-1] - old_grid.matrix_array[old_grid.size_x-2]);
1530  //this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-2] + a*(old_function.matrix_array[old_grid.size_x-1] - old_function.matrix_array[old_grid.size_x-2]);
1531  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-1])/(old_grid.matrix_array[old_grid.size_x-2] - old_grid.matrix_array[old_grid.size_x-1]);
1532  this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-1] + a*(old_function.matrix_array[old_grid.size_x-2] - old_function.matrix_array[old_grid.size_x-1]);
1533  /* } else if (i <= 1) { // !!!!!!!! always linear for the first two points ???
1534  int k = 1;
1535  while(old_grid.matrix_array[k] < new_grid.matrix_array[i]) {k++;}
1536  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k-1])/(old_grid.matrix_array[k] - old_grid.matrix_array[k-1]);
1537  this->matrix_array[i] = old_function.matrix_array[k-1] + a*(old_function.matrix_array[k] - old_function.matrix_array[k-1]);
1538  } else if (i >= new_grid.size_x-2) { // !!!!!!!! always linear for the last two points ???
1539  int k = new_grid.size_x-1;
1540  while(old_grid.matrix_array[k] > new_grid.matrix_array[i]) {k--;}
1541  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k])/(old_grid.matrix_array[k+1] - old_grid.matrix_array[k]);
1542  this->matrix_array[i] = old_function.matrix_array[k] + a*(old_function.matrix_array[k+1] - old_function.matrix_array[k]);*/
1543  } else {
1544  //splint(old_grid.matrix_array, old_function.matrix_array, y2.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i]);
1545  //splint(&old_grid.matrix_array[1], &old_function.matrix_array[1], &y2.matrix_array[1], old_grid.size_x-2, new_grid.matrix_array[i], &this->matrix_array[i]);
1546  // Shouldn't we pass whole matrix for the least argument????
1548  splint(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], &y2.matrix_array[i_min], interp_size, new_grid.matrix_array[i], &this->matrix_array[i]);
1549  }
1550  }
1551 
1552 }
1553 
1554 /*template<class T>
1555 void Matrix1D<T>::Spline(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double ub, double lb) {
1556 
1557 //if (size_x != old_grid.size_x) {
1558 if (this->size_x != new_grid.size_x) {
1559 cout << "new grid and function size mismatch!" << endl;
1560 return;
1561 }
1562 
1563 double yp0, ypn;
1564 //yp0 = 2*(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]) - (old_function[2] - old_function[1])/(old_grid[2] - old_grid[1]);
1565 //ypn = 2*(old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]) - (old_function[size_x-2] - old_function[size_x-3])/(old_grid[size_x-2] - old_grid[size_x-3]);
1566 // ypn = (old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1567 
1568 yp0 = 1e99;//-(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]);
1569 ypn = 1e99;//0;
1570 
1571 Matrix1D<double> y2(old_grid.size_x);
1572 //spline(old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, yp0, ypn, y2.matrix_array);
1573 spline(&old_grid.matrix_array[1], &old_function.matrix_array[1], old_grid.size_x-2, yp0, ypn, &y2.matrix_array[1]);
1574 
1575 int i;
1576 for (i = 0; i < new_grid.size_x; i++) {
1577 if (i == 0) {
1578 this->matrix_array[i] = lb;
1579 } else if (i == new_grid.size_x-1) {
1580 this->matrix_array[i] = ub;
1581 } else
1582 if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) {
1583 this->matrix_array[i] = lb;//old_function.matrix_array[0];
1584 } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) {
1585 // linear
1586 double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
1587 this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
1588 } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) {
1589 this->matrix_array[i] = ub;//1.e-21; //old_function.matrix_array[old_function.size_x-1];
1590 } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) {
1591 // linear
1592 //double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-2])/(old_grid.matrix_array[old_grid.size_x-1] - old_grid.matrix_array[old_grid.size_x-2]);
1593 //this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-2] + a*(old_function.matrix_array[old_grid.size_x-1] - old_function.matrix_array[old_grid.size_x-2]);
1594 double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-1])/(old_grid.matrix_array[old_grid.size_x-2] - old_grid.matrix_array[old_grid.size_x-1]);
1595 this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-1] + a*(old_function.matrix_array[old_grid.size_x-2] - old_function.matrix_array[old_grid.size_x-1]);
1596 } else {
1597 //splint(old_grid.matrix_array, old_function.matrix_array, y2.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i]);
1598 splint(&old_grid.matrix_array[1], &old_function.matrix_array[1], &y2.matrix_array[1], old_grid.size_x-2, new_grid.matrix_array[i], &this->matrix_array[i]);
1599 }
1600 }
1601 
1602 }*/
1603 
1620 template<class T>
1621 void Matrix1D<T>::Spline2(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub, double lin_spline_coef, double max_second_der) {
1622 
1623  //if (size_x != old_grid.size_x) {
1624  // check, if size of the grids are equal
1625  if (this->size_x != new_grid.size_x) {
1626  throw error_msg("SPLINE_INTERPOLATION", "New grid and function size mismatch!");
1627  }
1628  // if size == 1 - just copy the value
1629  if (this->size_x == 1) {
1630  this->matrix_array[0] = old_function.matrix_array[0];
1631  return;
1632  }
1633 
1634  // Idea was to find the intersection between grids, to interpolate inside the intersection and do not extrapolate outside. But is done later.
1635  int i=0, i_min=0, i_max=old_grid.size_x-1, interp_size=0;
1636  // for(i = 0; old_grid.matrix_array[i] <= new_grid.matrix_array[0] && i < old_grid.size_x-1; i++) {}
1637  // i_min = i;
1638  // for(i = old_grid.size_x-1; old_grid.matrix_array[i] >= new_grid.matrix_array[old_grid.size_x-1] && i > 0; i--) {}
1639  // i_max = i;
1640  //cout << i_min << " " << i_max << endl;
1641 
1642  // new range for spline interpolation - values close to the boundaries are interpolated by linear interpolation
1643  i_min = 0; i_max = old_grid.size_x-1;
1644  interp_size = i_max-i_min+1;
1645 
1646  double yp0, ypn;
1647  //yp0 = 2*(old_function[1] - old_function[0])/(old_grid[1] - old_grid[0]) - (old_function[2] - old_function[1])/(old_grid[2] - old_grid[1]);
1648  //ypn = 2*(old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]) - (old_function[size_x-2] - old_function[size_x-3])/(old_grid[size_x-2] - old_grid[size_x-3]);
1649  // ypn = (old_function[size_x-1] - old_function[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1650 
1651  // flag, means use smooth second derivative on the boundary
1652  yp0 = 1e99;
1653  //yp0 = -(old_function[2] - old_function[0])/(old_grid[2] - old_grid[0]);
1654  ypn = 1e99;//0;
1655 
1656  Matrix1D<T> y2(old_grid.size_x);
1657  // calculating array of second derivatives
1658  spline(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], interp_size, yp0, ypn, &y2.matrix_array[i_min]);
1659 
1660  // for all grid points - looking for new values
1661  for (i = 0; i < new_grid.size_x; i++) {
1662  splint(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], &y2.matrix_array[i_min], interp_size, new_grid.matrix_array[i], &this->matrix_array[i]);
1663  }
1664 
1665 }
1666 
1667 
1672 template<class T>
1673 void Matrix1D<T>::Polilinear(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub) {
1674 
1675  //if (size_x != old_grid.size_x) {
1676  if (this->size_x != new_grid.size_x) {
1677  throw error_msg("POLINOMIAL_INTERPOLATION", "New grid and function size mismatch!");
1678  }
1679 
1680  int i;
1681  for (i = 0; i < new_grid.size_x; i++)
1682  matrix_array[i] = polilinear( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], lb, ub);
1683 
1684 }
1685 
1690 template<class T>
1691 void Matrix1D<T>::Polint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1692 
1693  //if (size_x != old_grid.size_x) {
1694  if (this->size_x != new_grid.size_x) {
1695  throw error_msg("POLINT_INTERPOLATION", "New grid and function size mismatch!");
1696  }
1697 
1698  int i;
1699  double err;
1700  for (i = 0; i < new_grid.size_x; i++)
1701  polint( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i], &err);
1702 
1703 }
1704 
1706 template<class T>
1707 void Matrix1D<T>::Ratint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1708 
1709  //if (size_x != old_grid.size_x) {
1710  if (this->size_x != new_grid.size_x) {
1711  throw error_msg("RATIONAL_INTERPOLATION", "New grid and function size mismatch!");
1712  }
1713 
1714  int i;
1715  double err;
1716  for (i = 0; i < new_grid.size_x; i++)
1717  ratint( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i], &err);
1718 
1719 }
1720 
1721 //template<class T>
1722 //Matrix1D<T> Matrix1D<T>::Spline(Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
1723 //
1724 // if (this->size_x != old_grid.size_x) {
1725 // //if (this->size_x != new_grid.size_x) {
1726 // cout << "old grid and function size mismatch!" << endl;
1727 // return -1;
1728 // }
1729 //
1730 // double yp0, ypn;
1731 // yp0 = (matrix_array[0] - matrix_array[1])/(old_grid[0] - old_grid[1]);
1732 // ypn = (matrix_array[size_x-1] - matrix_array[size_x-2])/(old_grid[size_x-1] - old_grid[size_x-2]);
1733 // Matrix1D<T> y2(old_grid.size_x);
1734 // spline(old_grid.matrix_array, this->matrix_array, old_grid.size_x, yp0, ypn, y2.matrix_array)
1735 //
1736 // int i;
1737 // Matrix1D<T> res(new_grid.size_x);
1738 // for (i = 0; i < new_grid.size_x; i++)
1739 // spline(old_grid.matrix_array, this->matrix_array, y2.matrix_array, new_grid.matrix_array[i], &res[i]);
1740 //
1741 // return res;
1742 //
1743 //}
1744 
1745 
1746 // interpolations
1747 
1751 template<class T>
1752 void Matrix1D<T>::Akima_interpolation(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, int extrapolation_type, double lb, double ub) {
1753 
1754  /*
1755  for (i = 0; i < new_grid.size_x; i++) {
1756  if (i == 0) { // if we are on the boundary
1757  this->matrix_array[i] = lb;
1758  } else if (i == new_grid.size_x-1) { // if we are on the boundary
1759  this->matrix_array[i] = ub;
1760  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) { // if it is extrapolation
1761  this->matrix_array[i] = lb;
1762  //this->matrix_array[i] = old_function.matrix_array[0];
1763  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) { // if it is extrapolation
1764  this->matrix_array[i] = ub;
1765  //this->matrix_array[i] = old_function.matrix_array[old_function.size_x-1];
1766  } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) { // if the new grid point is between two boundaryes poins of the old grid
1767  // linear
1768  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
1769  this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
1770  } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) { // if the new grid point is between two boundaryes poins of the old grid
1771  // linear
1772  //double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-2])/(old_grid.matrix_array[old_grid.size_x-1] - old_grid.matrix_array[old_grid.size_x-2]);
1773  //this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-2] + a*(old_function.matrix_array[old_grid.size_x-1] - old_function.matrix_array[old_grid.size_x-2]);
1774  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[old_grid.size_x-1])/(old_grid.matrix_array[old_grid.size_x-2] - old_grid.matrix_array[old_grid.size_x-1]);
1775  this->matrix_array[i] = old_function.matrix_array[old_grid.size_x-1] + a*(old_function.matrix_array[old_grid.size_x-2] - old_function.matrix_array[old_grid.size_x-1]);
1776  } else if (i <= 1) { // always linear for the first point?
1777  int k = 1;
1778  // Wrong! We need to compare to grids, not indexes. Changed in ver 1.04.
1779  // while(k < i) {k++;}
1780  while(old_grid.matrix_array[k] < new_grid.matrix_array[i]) {k++;}
1781  double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k-1])/(old_grid.matrix_array[k] - old_grid.matrix_array[k-1]);
1782  this->matrix_array[i] = old_function.matrix_array[k-1] + a*(old_function.matrix_array[k] - old_function.matrix_array[k-1]);
1783  } else {
1784  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, ub, lb);
1785  }
1786  }*/
1787 
1788 
1789  // add point
1790  /*int i_first = 0, i_last = new_grid.size_x-1;
1791  while(old_grid.matrix_array[i_first] <= new_grid.matrix_array[0]) {i_first++;}
1792  while(old_grid.matrix_array[i_last] >= new_grid.matrix_array[new_grid.size_x-1]) {i_last--;}
1793  Matrix1D<double> mod_grid(i_last-i_first+2);
1794  mod_grid[0] = lb;
1795  mod_grid[mod_grid.x_size-1] = ub;
1796 
1797  for (int i = 1; i < mod_grid.size_x-1; i++) {
1798  mod_grid[i] = old_grid.matrix_array[i_first + i - 1];
1799  }
1800 
1801  akima(mod_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);*/
1802 
1803 
1804  /*int i_first = 0, i_last = new_grid.size_x-1;
1805  while(new_grid.matrix_array[i_first] <= old_grid.matrix_array[0]) {i_first++;}
1806  while(new_grid.matrix_array[i_last] >= old_grid.matrix_array[old_grid.size_x-1]) {i_last--;}
1807  int new_size_x = i_last - i_first + 1;
1808 
1809  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x,
1810  &new_grid.matrix_array[i_first], &this->matrix_array[i_first], new_size_x,
1811  2,
1812  ub, lb);
1813 
1814  // boundaries (and extrapolation)
1815  int i;
1816  for (i = 0; i < i_first; i++) this->matrix_array[i] = lb;
1817  for (i = this->size_x-1; i > i_last; i--) this->matrix_array[i] = ub;*/
1818 
1819  // new grid is exactly one point wider than old grid
1820  /*int old_first = 0, old_last = old_grid.size_x-1;
1821  while(old_grid.matrix_array[old_first] <= new_grid.matrix_array[0]) {old_first++;}
1822  while(old_grid.matrix_array[old_last] >= new_grid.matrix_array[new_grid.size_x-4]) {old_last--;}
1823  int old_size_x = old_last - old_first + 1;
1824  int new_first = 0, new_last = new_grid.size_x-1;
1825  while(new_grid.matrix_array[new_first+1] <= old_grid.matrix_array[0]) {new_first++;}
1826  while(new_grid.matrix_array[new_last-1] >= old_grid.matrix_array[old_grid.size_x-4]) {new_last--;}
1827  int new_size_x = new_last - new_first + 1;
1828 
1829  akima(&old_grid.matrix_array[old_first], &old_function.matrix_array[old_first], old_size_x,
1830  &new_grid.matrix_array[new_first], &this->matrix_array[new_first], new_size_x, 1, ub, lb);
1831 
1832  // boundaries (and extrapolation)
1833  int i;
1834  for (i = 0; i < new_first; i++) this->matrix_array[i] = lb;
1835  for (i = this->size_x-1; i > new_last; i--) this->matrix_array[i] = ub;*/
1836 
1838  /* int i_first = 0, i_last = new_grid.size_x-1;
1839  while(old_grid.matrix_array[i_first] <= new_grid.matrix_array[0]) {i_first++;}
1840  while(old_grid.matrix_array[i_last] >= new_grid.matrix_array[new_grid.size_x-1]) {i_last--;}
1841  int new_size_x = i_last - i_first + 1;
1842 
1843  //akima(&old_grid.matrix_array[i_first], &old_function.matrix_array[i_first], new_size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);
1844  // akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);
1845  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x,
1846  &new_grid.matrix_array[i_first], &this->matrix_array[i_first], new_size_x,
1847  2,
1848  ub, lb,
1849  max(this->matrix_array[0], old_grid.matrix_array[0]), min(this->matrix_array[this->x_size-1], old_grid.matrix_array[old_grid.x_size-1]));
1850 
1851  // boundaries (and extrapolation)
1852  int i;
1853  for (i = 0; i < i_first; i++) this->matrix_array[i] = lb;
1854  for (i = this->size_x-1; i > i_last; i--) this->matrix_array[i] = ub;*/
1855 
1856 
1857 
1858  akima(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, 1, ub, lb);
1859  matrix_array[0] = lb;
1860  matrix_array[size_x-1] = ub;
1861 }
1862 
1863 
1867 template<class T>
1868 void Matrix1D<T>::Akima_interpolation2(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, int extrapolation_type, double lb, double ub) {
1869  akima2(old_grid.matrix_array, old_function.matrix_array, old_function.size_x, new_grid.matrix_array, this->matrix_array, this->size_x, ub, lb);
1870 }
1871 
1876 template<class T>
1877 double Linear2D(double x, double y, T &old_grid_x, T &old_grid_y, T &f) {
1878 
1879  int i = 0, j = 0;
1880 
1881  while (x >= old_grid_x[++i][0] && i < old_grid_y.size_y-1);
1882  while (y >= old_grid_y[0][++j] && j < old_grid_x.size_x-1);
1883 
1884  double t = (x - old_grid_x[i-1][0])/(old_grid_x[i][0] - old_grid_x[i-1][0]);
1885  double u = (y - old_grid_y[0][j-1])/(old_grid_y[0][j] - old_grid_y[0][j-1]);
1886 
1887  return (1.0-t)*(1.0-u)*f[i-1][j-1] + t*(1.0-u)*f[i][j-1] + t*u*f[i][j] + (1.0-t)*u*f[i-1][j];
1888 }
1889 
1890 
1896 template<class T>
1897 void Matrix2D<T>::Interpolate(Matrix2D<T> &old_function, Matrix2D<T> &old_grid_x, Matrix2D<T> &old_grid_y, Matrix2D<T> &new_grid_x, Matrix2D<T> &new_grid_y) {
1898 
1899  Matrix2D<T> res(new_grid_x.size_x, new_grid_y.size_y);
1900 
1901  if (old_function.size_x != old_grid_x.size_x || old_function.size_y != old_grid_y.size_y) {
1902  throw error_msg("LINEAR_INTERPOLATION", "Old grid and function size mismatch!");
1903  }
1904 
1905  if (this->size_x != new_grid_x.size_x || this->size_y != new_grid_y.size_y) {
1906  throw error_msg("LINEAR_INTERPOLATION", "New grid and function size mismatch!");
1907  }
1908 
1909  int i, j;
1910  for (i = 0; i < new_grid_x.size_x; i++) {
1911  for (j = 0; j < new_grid_y.size_y; j++) {
1912  this->matrix_array[i][j] = Linear2D(new_grid_x[i][j], new_grid_y[i][j], old_grid_x, old_grid_y, old_function);
1913  }
1914  }
1915 
1916  return;
1917 
1918 }
1919 
1920 
1928 CalculationMatrix::CalculationMatrix(int x_size, int y_size, int z_size, int n_of_diags) {
1929  this->initialized = false;
1930  Initialize(x_size, y_size, z_size, n_of_diags);
1931 }
1932 
1936 void CalculationMatrix::Initialize(int x_size, int y_size, int z_size, int n_of_diags) {
1937  this->initialized = false;
1938  this->x_size = x_size;
1939  this->y_size = y_size;
1940  this->total_size = x_size;
1941  if (y_size > 0) this->total_size = this->total_size * y_size;
1942  if (z_size > 0) this->total_size = this->total_size * z_size;
1943 
1944 
1945  // initializing diagonals
1946  int x, y, z, id;
1947  // !!! if (z_size > 0) { // means 3d
1948  if (z_size > 1) { // means 3d
1949  for (x = -n_of_diags; x <= n_of_diags; x++) {
1950  for (y = -n_of_diags; y <= n_of_diags; y++) {
1951  for (z = -n_of_diags; z <= n_of_diags; z++) {
1952  // calculating a diagonal number (id)
1953  id = this->index1d(x, y, z);
1954  // allocating memory for each diagonal
1955  (*this)[id] = Matrix1D<double>(this->total_size);
1956  (*this)[id] = 0;
1957  }
1958  }
1959  }
1960  // !!! } else if (y_size > 0) {
1961  } else if (y_size > 1) {
1962  for (x = -n_of_diags; x <= n_of_diags; x++) {
1963  for (y = -n_of_diags; y <= n_of_diags; y++) {
1964  // calculating a diagonal number (id)
1965  id = this->index1d(x, y);
1966  // allocating memory for each diagonal
1967  (*this)[id] = Matrix1D<double>(this->total_size);
1968  (*this)[id] = 0;
1969  }
1970  }
1971  // !!! } else if (x_size > 0) {
1972  } else if (x_size > 1) {
1973  for (x = -n_of_diags; x <= n_of_diags; x++) {
1974  // calculating a diagonal number (id)
1975  id = this->index1d(x);
1976  // allocating memory for each diagonal
1977  (*this)[id] = Matrix1D<double>(this->total_size);
1978  (*this)[id] = 0;
1979  }
1980  }
1981 
1982  // If im_size is not set, then it's 1d and ia_size should be zero for future
1983  if (y_size == 0) this->x_size = 0;
1984  // If il_size is not set, then it's 2d and im_size should be zero for future
1985  if (z_size == 0) this->y_size = 0;
1986 
1987  this->initialized = true;
1988 }
1989 
1993 int CalculationMatrix::index1d(int x, int y, int z) {
1994  return (z*y_size + y)*x_size + x;
1995 }
1996 
2000 void CalculationMatrix::writeToFile(string filename) {
2001  int in;
2002 
2003  ofstream output(filename.c_str());
2004  if (output==NULL && (filename.find("Debug") == string::npos)) {
2005  throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
2006  }
2007  output << "VARIABLES = ";
2008  for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2009  output << "\"" << it->first << "\", ";
2010  }
2011  output << endl;
2012  output << "ZONE T=\"" << "\", I=" << total_size << endl;
2013  output.setf(ios_base::scientific, ios_base::floatfield);
2014  for (in = 0; in < this->total_size; in++) {
2015  for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2016  output << "\t" << it->second[in];
2017  }
2018  output << endl;
2019  }
2020  output.close();
2021 }
2022 
2023 
2025 // Implementations
2027 
2028 template class Matrix1D<double>;
2029 template class Matrix2D<double>;
2030 template class Matrix3D<double>;
2031 
2032 
2033 #endif
void readFromFile(string filename)
Definition: Matrix.cpp:724
bool initialized
Flag, equal true if initialized.
Definition: Matrix.h:138
void Initialize(int x_size, int y_size=1, int z_size=1, int n_of_diags=1)
Definition: Matrix.cpp:1936
int size_y
size x, size y
Definition: Matrix.h:141
void Spline(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, double lb, double ub, double lin_spline_coef=0, double max_second_der=0)
Definition: Matrix.cpp:1455
Matrix1D times(const Matrix1D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
Definition: Matrix.cpp:310
static const double m
mass of electron in grams
void Polilinear(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, double lb, double ub)
Definition: Matrix.cpp:1673
const Matrix3D & operator+() const
Return itself as positive version of values.
Definition: Matrix.h:252
Matrix3D()
Default constructor. Do nothing.
Definition: Matrix.h:231
const Matrix3D operator-() const
Return negative version of values.
Definition: Matrix.h:253
Matrix3D operator*(const T Val) const
Definition: Matrix.cpp:1105
string name
name of the Matrix
Definition: Matrix.h:47
Matrix2D< T > ySlice(int p_y) const
Definition: Matrix.cpp:1373
bool initialized
Flag, equal true if initialized.
Definition: Matrix.h:45
T dot(const Matrix1D< T > &M) const
Dot product.
Definition: Matrix.cpp:349
Matrix2D max_of(T val)
Definition: Matrix.cpp:654
Matrix1D & operator=(const Matrix1D< T > &M)
Definition: Matrix.cpp:241
void writeToFile(string filename)
Save matrix to a file.
Definition: Matrix.cpp:1164
Matrix1D operator*(const T Val) const
Definition: Matrix.cpp:282
Matrix2D times(const Matrix2D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
Definition: Matrix.cpp:639
void AllocateMemory(int size_x, int size_y)
Definition: Matrix.cpp:524
Matrix3D operator/(const T Val) const
Definition: Matrix.cpp:1120
void Ratint(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
rational
Definition: Matrix.cpp:1707
Matrix3D & times_equal(const Matrix3D< T > &M)
Arraywise multiplication (A.*B), stores result in the matrix it's applied to.
Definition: Matrix.cpp:1050
void Akima_interpolation2(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, int extrapolation_type=0, double lb=0, double ub=0)
Definition: Matrix.cpp:1868
int size_x
size x, size y
Definition: Matrix.h:141
Linear interpolation of two arrays.
Definition: linear.h:16
T norm() const
Norm.
Definition: Matrix.cpp:336
Matrix3D & divide_equal(const Matrix3D< T > &M)
Arraywise division (A./B), stores result in the matrix it's applied to.
Definition: Matrix.cpp:1063
Matrix3D divide(const Matrix3D< T > &M) const
Arraywise division (A./B), stores result in a new matrix.
Definition: Matrix.cpp:1148
void AllocateMemory(int x_size)
Definition: Matrix.cpp:173
bool initialized
Flag, equal true if initialized.
Definition: Matrix.h:222
int index1d(int x, int y=0, int z=0)
Definition: Matrix.cpp:1993
void splint(double *xa, double *ya, double *y2a, int n, double x, double *y)
Definition: spline.cpp:86
void AllocateMemory(int size_x, int size_y, int size_z)
Definition: Matrix.cpp:841
string name
name of the Matrix
Definition: Matrix.h:227
Matrix3D< T > abs()
Definition: Matrix.cpp:1339
int index1d(int x, int y) const
Definition: Matrix.cpp:667
void free_matrix(T *m)
Definition: Matrix.cpp:86
one dimensional matrix class
Definition: Matrix.h:41
void Polint(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
Definition: Matrix.cpp:1691
void polint(double *xa, double *ya, int n, double x, double *y, double *dy)
Definition: polint.cpp:15
Matrix3D times(const Matrix3D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
Definition: Matrix.cpp:1134
Matrix2D operator/(const T Val) const
Definition: Matrix.cpp:609
~Matrix1D()
Definition: Matrix.cpp:163
T & operator()(int x, int y, int z)
Return the (x,y,z) value of matrix.
Definition: Matrix.cpp:880
Matrix3D & operator+=(const Matrix3D< T > &M)
Definition: Matrix.cpp:972
T ** operator[](int i)
Return the i-th pointer to 2d-array. Next [j][k] can be applied, so we have regular [i][j][k]...
Definition: Matrix.cpp:864
void Spline2(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, double lb, double ub, double lin_spline_coef=0, double max_second_der=0)
Definition: Matrix.cpp:1621
int size_y
size x, size y, size z
Definition: Matrix.h:225
T & operator[](int i)
Return the i-th value of matrix.
Definition: Matrix.cpp:189
Matrix3D & operator=(const Matrix3D< T > &M)
Definition: Matrix.cpp:899
T maxabs()
Definition: Matrix.cpp:1322
void writeToFile(string filename)
Definition: Matrix.cpp:366
~Matrix2D()
Definition: Matrix.cpp:513
void Interpolate(Matrix2D< T > &old_function, Matrix2D< T > &old_grid_x, Matrix2D< T > &old_grid_y, Matrix2D< T > &new_grid_x, Matrix2D< T > &new_grid_y)
Definition: Matrix.cpp:1897
Matrix2D & operator=(const Matrix2D< T > &M)
Definition: Matrix.cpp:541
Matrix1D divide(const Matrix1D< T > &M) const
Arraywise division (A./B), stores result in a new matrix.
Definition: Matrix.cpp:324
T * matrix_array
Array to keep the values.
Definition: Matrix.h:43
double T(double alpha)
Function for bounce time.
void Interpolate(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
Definition: Matrix.cpp:1414
Matrix2D< T > xSlice(int p_x) const
Definition: Matrix.cpp:1357
double Linear2D(double x, double y, T &old_grid_x, T &old_grid_y, T &f)
Definition: Matrix.cpp:1877
three dimensional matrix class
Definition: Matrix.h:215
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
Definition: spline.cpp:30
void ratint(double *xa, double *ya, int n, double x, double *y, double *dy)
Definition: ratint.cpp:16
Matrix2D divide(const Matrix2D< T > &M) const
Arraywise division (A./B), stores result in a new matrix.
Definition: Matrix.cpp:624
two dimensional matrix class
Definition: Matrix.h:132
Error message - stack of single_errors.
Definition: error.h:54
double getValue(double x)
Definition: linear.h:47
Matrix3D & operator-=(const Matrix3D< T > &M)
Definition: Matrix.cpp:985
void readFromFile(string filename)
Load matrix from a file.
Definition: Matrix.cpp:1208
T * matrix(long Rows)
Allocating memory for 1D matrix.
Definition: Matrix.cpp:21
void readFromFile(string filename)
Definition: Matrix.cpp:405
T max()
Definition: Matrix.cpp:1305
string name
name of the Matrix
Definition: Matrix.h:143
void writeToFile(string filename)
Definition: Matrix.cpp:682
~Matrix3D()
Definition: Matrix.cpp:832
Matrix2D operator*(const T Val) const
Definition: Matrix.cpp:594
void writeToFile(string filename)
Definition: Matrix.cpp:2000
int index1d(int x, int y, int z)
Returns index of the element (x,y,z) in 1d array.
Definition: Matrix.cpp:1296
double polilinear(double *xa, double *ya, int n, double x, double lb, double ub)
Definition: polilinear.cpp:13
int size_x
size x, size y, size z
Definition: Matrix.h:225
int size_x
size x
Definition: Matrix.h:46
T & operator()(int x)
Return the x-th value of matrix.
Definition: Matrix.cpp:224
Matrix2D< T > zSlice(int p_z) const
Definition: Matrix.cpp:1389
Matrix1D operator/(const T Val) const
Definition: Matrix.cpp:296
Matrix3D & operator*=(const T Val)
Definition: Matrix.cpp:998
void Akima_interpolation(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid, int extrapolation_type=0, double lb=0, double ub=0)
Definition: Matrix.cpp:1752
Matrix3D & operator/=(const T Val)
Definition: Matrix.cpp:1011
int size_z
size x, size y, size z
Definition: Matrix.h:225