00001
00006 #ifndef matrix_array_MATRIX_CPP
00007 #define matrix_array_MATRIX_CPP
00008
00009 #include "Matrix.h"
00010
00011 using namespace std;
00012
00013 #if defined(_WIN32) || defined(_WIN64)
00014 #define strncasecmp _strnicmp
00015 #define strcasecmp _stricmp
00016 #endif
00017
00018
00019
00021 template<class T>inline T* matrix(long Rows)
00022 {
00023 T *m=new T[Rows];
00024
00025 if (m == NULL) {
00026 throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", Rows*sizeof(T));
00027 }
00028 return m;
00029 }
00030
00032 template<class T>inline T** matrix(long Rows, long Columns)
00033 {
00034
00035 T **m=new T*[Rows];
00036
00037 if (m == NULL) {
00038 throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", Rows * sizeof(T));
00039 }
00040
00041 m[0] = new T[Rows * Columns];
00042
00043 if (m[0] == NULL) {
00044 throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", Rows * Columns * sizeof(T));
00045 }
00046
00047 for(long i=1; i<Rows; i++) m[i] = m[i-1] + Columns;
00048 return m;
00049 }
00050
00052 template<class T>inline T*** matrix(int size_x, int size_y, int size_z)
00053 {
00054
00055 T ***m=new T**[size_x];
00056
00057 if (m == NULL) {
00058 throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", size_x * sizeof(T));
00059 }
00060 for (int x = 0; x < size_x; x++) {
00061
00062 m[x] = new T*[size_y];
00063
00064 if (m[x] == NULL) {
00065 throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", size_y * sizeof(T));
00066 }
00067 }
00068
00069 m[0][0] = new T[size_x * size_y * size_z];
00070
00071 if (m[0][0] == NULL) {
00072 throw error_msg("MEMORY_ERRROR", "Memory can't be initialized: %s size", size_x * size_y * size_z * sizeof(T));
00073 }
00074 for (int x = 0; x < size_x; x++) {
00075 for (int y = 0; y < size_y; y++) {
00076
00077 m[x][y] = m[0][0] + (x*size_y + y)*size_z;
00078 }
00079 }
00080
00081 return m;
00082 }
00083
00085 template<class T>inline void free_matrix(T* m) {
00086 delete m;
00087 }
00088
00090 template<class T>inline void free_matrix(T** m) {
00091 delete[](m[0]);
00092 delete[](m);
00093 }
00094
00096 template<class T>inline void free_matrix(T*** m, int size_x, int size_y) {
00097 delete[](m[0][0]);
00098 for (int x = 0; x < size_x; x++) {
00099 delete[](m[x]);
00100 }
00101 delete[](m);
00102 }
00103
00105
00106
00107
00109
00119 template<class T>
00120 Matrix1D<T>::Matrix1D( int x_size , string name) {
00121 initialized = false;
00122 this->name = name;
00123
00124 AllocateMemory(x_size);
00125 }
00126
00135 template<class T>
00136 Matrix1D<T>::Matrix1D( int x_size ) {
00137 initialized = false;
00138
00139 AllocateMemory(x_size);
00140 }
00141
00150 template<class T>
00151 Matrix1D<T>::Matrix1D( const Matrix1D<T> &M ) {
00152 initialized = false;
00153 this->operator = (M);
00154 }
00155
00159 template<class T>
00160 Matrix1D<T>::~Matrix1D() {
00161 if (initialized) free_matrix<T>(matrix_array);
00162 }
00163
00169 template<class T>
00170 void Matrix1D<T>::AllocateMemory( int x_size ) {
00171 this->size_x = x_size;
00172
00173 matrix_array = matrix<T>(x_size);
00174 initialized = true;
00175
00176 }
00177
00178
00185 template<class T>
00186 T& Matrix1D<T>::operator[](int i) {
00187 #ifdef DEBUG_MODE
00188 if (!initialized) {
00189 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00190 }
00191 #endif
00192
00193 return matrix_array[i];
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00213 template<class T>
00214 T& Matrix1D<T>::operator()(int x) {
00215 #ifdef DEBUG_MODE
00216 if (!initialized) {
00217 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00218 }
00219 #endif
00220
00221 return matrix_array[x];
00222 }
00223
00230 template<class T>
00231 Matrix1D<T>& Matrix1D<T>::operator= (const Matrix1D<T> &M) {
00232
00233 if (this!=&M) {
00234
00235 if (M.initialized) {
00236
00237 if (initialized) free_matrix<T>(matrix_array);
00238
00239 this->size_x = M.size_x;
00240 this->name = M.name;
00241
00242 this->AllocateMemory(this->size_x);
00243
00244 memcpy( matrix_array, M.matrix_array, this->size_x * sizeof( T ) );
00245 } else {
00246 this->initialized = false;
00247 }
00248 }
00249 return *this;
00250 }
00251
00258 template<class T>
00259 Matrix1D<T>& Matrix1D<T>::operator= (T Val) {
00260 for (int x = 0; x < this->size_x; x++)
00261 matrix_array[x] = Val;
00262 return *this;
00263 }
00264
00271 template<class T>
00272 Matrix1D<T> Matrix1D<T>::operator/ (T Val) {
00273 int x;
00274 Matrix1D<T> Tmp(*this);
00275 for (x = 0; x < this->size_x; x++)
00276 Tmp[x] = matrix_array[x] / Val;
00277 return Tmp;
00278 }
00279
00286 template<class T>
00287 Matrix1D<T> Matrix1D<T>::operator* (T Val) {
00288 Matrix1D<T> Tmp(*this);
00289 for (int x = 0; x < this->size_x; x++)
00290 Tmp[x] = matrix_array[x] * Val;
00291 return Tmp;
00292 }
00293
00294
00300 template<class T>
00301 Matrix1D<T> Matrix1D<T>::operator/ (Matrix1D<T> &M) {
00302 int x;
00303 Matrix1D<T> Tmp(*this);
00304 for (x = 0; x < this->size_x; x++)
00305 Tmp[x] = matrix_array[x] / M[x];
00306 return Tmp;
00307 }
00308
00314 template<class T>
00315 Matrix1D<T> Matrix1D<T>::operator* (Matrix1D<T> &M) {
00316 int x;
00317 Matrix1D<T> Tmp(*this);
00318 for (x = 0; x < this->size_x; x++)
00319 Tmp[x] = matrix_array[x] * M[x];
00320 return Tmp;
00321 }
00322
00326 template<class T>
00327 void Matrix1D<T>::writeToFile(string filename) {
00328 int x;
00329 ofstream output(filename.c_str());
00330 if (output==NULL && (filename.find("Debug") == string::npos)) {
00331 throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
00332 }
00333 output << "VARIABLES = \"" << ((this->name!="")?this->name:"function") << "\" "<< endl;
00334 output << "ZONE T=\"" << filename << "\", I=" << size_x << endl;
00335 output.setf(ios_base::scientific, ios_base::floatfield);
00336 for (x = 0; x < size_x; x++) {
00337 output << matrix_array[x] << endl;
00338 }
00339 output.close();
00340 }
00341
00345 template<class T>
00346 void Matrix1D<T>::writeToFile(string filename, Matrix1D<T> &grid_x) {
00347 int x;
00348 ofstream output(filename.c_str());
00349 if (output==NULL && (filename.find("Debug") == string::npos)) {
00350 throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
00351 }
00352 output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((this->name!="")?this->name:"function") << "\" "<< endl;
00353 output << "ZONE T=\"" << filename << "\", I=" << size_x << endl;
00354 output.setf(ios_base::scientific, ios_base::floatfield);
00355 for (x = 0; x < size_x; x++) {
00356 output << "\t" << grid_x[x] << "\t" << "\t" << matrix_array[x] << endl;
00357 }
00358 output.close();
00359 }
00360
00361
00365 template<class T>
00366 void Matrix1D<T>::readFromFile(string filename) {
00367 int x;
00368 string inBuf;
00369 if (!initialized) {
00370 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00371 } else {
00372 ifstream input(filename.c_str());
00373 if (input != NULL && !input.eof()) {
00374
00375 input >> inBuf;
00376
00377 while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
00378 input >> inBuf;
00379 }
00380
00381 input.ignore(9999, '\n');
00382 for (x = 0; x < size_x; x++) {
00383 input >> matrix_array[x];
00384 }
00385 } else {
00386 throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
00387
00388
00389
00390 }
00391 input.close();
00392 }
00393 }
00394
00398 template<class T>
00399 void Matrix1D<T>::readFromFile(string filename, Matrix1D<T> &grid_x) {
00400 int x;
00401 string inBuf;
00402 double loaded_x;
00403 double err = 1e-8;
00404 if (!initialized) {
00405 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00406 } else {
00407 ifstream input(filename.c_str());
00408 if (input != NULL && !input.eof()) {
00409
00411 input >> inBuf;
00412
00413 while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
00414 input >> inBuf;
00415 }
00416
00417 input.ignore(9999, '\n');
00418 for (x = 0; x < size_x; x++) {
00419 input >> loaded_x;
00420
00421 if (fabs(loaded_x - grid_x[x]) > err) {
00422 throw error_msg("MATRIX_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: %e\nGrid: %e\n", filename.c_str(), loaded_x, grid_x[x]);
00423 } else {
00424 input >> matrix_array[x];
00425 }
00426 }
00427 } else {
00428 throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
00429
00430
00431
00432 }
00433 input.close();
00434 }
00435 }
00436
00437
00439
00440
00441
00443
00451 template<class T>
00452 Matrix2D<T>::Matrix2D( int x_size, int y_size ) {
00453 initialized = false;
00454
00455 AllocateMemory(x_size, y_size);
00456 }
00457
00464 template<class T>
00465 Matrix2D<T>::Matrix2D( const Matrix2D<T> &M ) {
00466 initialized = false;
00467 this->operator = (M);
00468 }
00469
00473 template<class T>
00474 Matrix2D<T>::~Matrix2D() {
00475 if (initialized) free_matrix<T>(matrix_array);
00476 }
00477
00484 template<class T>
00485 void Matrix2D<T>::AllocateMemory( int x_size, int y_size ) {
00486 this->size_x = x_size;
00487 this->size_y = y_size;
00488
00489 matrix_array = matrix<T>(x_size, y_size);
00490 initialized = true;
00491
00492
00493
00494 }
00495
00499 template<class T>
00500 int Matrix2D<T>::index1d(int x, int y) {
00501 return x*size_y + y;
00502 }
00503
00504
00510 template<class T>
00511 Matrix2D<T>& Matrix2D<T>::operator= (const Matrix2D<T> &M) {
00512
00513 if (this!=&M) {
00514
00515 if (M.initialized) {
00516
00517 if (initialized) free_matrix<T>(matrix_array);
00518 this->size_x = M.size_x;
00519 this->size_y = M.size_y;
00520 this->name = M.name;
00521
00522 this->AllocateMemory(this->size_x, this->size_y);
00523
00524
00525
00526 memcpy(matrix_array[0], M.matrix_array[0], this->size_x * this->size_y * sizeof(T) );
00527 } else {
00528 this->initialized = false;
00529 }
00530 }
00531 return *this;
00532 }
00533
00539 template<class T>
00540 Matrix2D<T>& Matrix2D<T>::operator= (T val) {
00541 int x, y;
00542 if (initialized) {
00543 for (x = 0; x < this->size_x; x++)
00544 for (y = 0; y < this->size_y; y++)
00545 matrix_array[x][y] = val;
00546 } else {
00547 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00548 }
00549
00550 return *this;
00551 }
00552
00558 template<class T>
00559 Matrix2D<T> Matrix2D<T>::operator/ (T Val) {
00560 int x, y;
00561 Matrix2D<T> Tmp(*this);
00562 for (x = 0; x < this->size_x; x++)
00563 for (y = 0; y < this->size_y; y++)
00564 Tmp[x][y] = matrix_array[x][y] / Val;
00565 return Tmp;
00566 }
00567
00573 template<class T>
00574 Matrix2D<T> Matrix2D<T>::operator* (T Val) {
00575 int x, y;
00576 Matrix2D<T> Tmp(*this);
00577 for (x = 0; x < this->size_x; x++)
00578 for (y = 0; y < this->size_y; y++)
00579 Tmp[x][y] = matrix_array[x][y] * Val;
00580 return Tmp;
00581 }
00582
00588 template<class T>
00589 Matrix2D<T> Matrix2D<T>::operator/ (Matrix2D<T> &M) {
00590 int x, y;
00591 Matrix2D<T> Tmp(*this);
00592 for (x = 0; x < this->size_x; x++)
00593 for (y = 0; y < this->size_y; y++)
00594 Tmp[x][y] = matrix_array[x][y] / M[x][y];
00595 return Tmp;
00596 }
00597
00603 template<class T>
00604 Matrix2D<T> Matrix2D<T>::operator* (Matrix2D<T> &M) {
00605 int x, y;
00606 Matrix2D<T> Tmp(*this);
00607 for (x = 0; x < this->size_x; x++)
00608 for (y = 0; y < this->size_y; y++)
00609 Tmp[x][y] = matrix_array[x][y] * M[x][y];
00610 return Tmp;
00611 }
00612
00618 template<class T>
00619 Matrix2D<T> Matrix2D<T>::max(T val) {
00620 int x, y;
00621 Matrix2D<T> Tmp(*this);
00622 for (x = 0; x < this->size_x; x++)
00623 for (y = 0; y < this->size_y; y++)
00624 Tmp[x][y] = (matrix_array[x][y]>val)?matrix_array[x][y]:val;
00625 return Tmp;
00626 }
00627
00636 template<class T>
00637 void Matrix2D<T>::writeToFile(string filename) {
00638 int x, y;
00639 ofstream output(filename.c_str());
00640 if (output==NULL && (filename.find("Debug") == string::npos)) {
00641 throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
00642 }
00643 output << "VARIABLES = \""<< ((this->name!="")?this->name:"f") << "\" "<< endl;
00644 output << "ZONE T=\"" << filename << "\", I=" << size_y << ", J= " << size_x << endl;
00645 for (x = 0; x < size_x; x++) {
00646 for (y = 0; y < size_y; y++) {
00647 output << matrix_array[x][y] << endl;
00648 }
00649 }
00650 output.close();
00651 }
00652
00653
00658 template<class T>
00659 void Matrix2D<T>::writeToFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y) {
00660 int x, y;
00661 ofstream output(filename.c_str());
00662 output << "VARIABLES = \"" << ((grid_x.name!="")?grid_x.name:"x") << "\", \"" << ((grid_y.name!="")?grid_y.name:"y") << "\", \"" << ((this->name!="")?this->name:"f") << "\" "<< endl;
00663 output << "ZONE T=\"" << filename << "\", I=" << size_y << ", J=" << size_x << endl;
00664 for (x = 0; x < size_x; x++) {
00665 for (y = 0; y < size_y; y++) {
00666 output << "\t" << grid_x[x][y] << "\t" << grid_y[x][y] << "\t" << matrix_array[x][y] << endl;
00667 }
00668 }
00669 output.close();
00670 }
00671
00672
00676 template<class T>
00677 void Matrix2D<T>::readFromFile(string filename) {
00678 int x, y;
00679 string inBuf;
00680 if (!initialized) {
00681 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00682 } else {
00683 ifstream input(filename.c_str());
00684 if (input != NULL && !input.eof()) {
00685
00686 input >> inBuf;
00687
00688 while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
00689 input >> inBuf;
00690 }
00691
00692 input.ignore(9999, '\n');
00693 for (x = 0; x < size_x; x++) {
00694 for (y = 0; y < size_y; y++) {
00695 input >> matrix_array[x][y];
00696 }
00697 }
00698 } else {
00699 throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
00700
00701
00702
00703 }
00704 input.close();
00705 }
00706 }
00707
00711 template<class T>
00712 void Matrix2D<T>::readFromFile(string filename, Matrix2D<T> &grid_x, Matrix2D<T> &grid_y) {
00713 int x, y;
00714 string inBuf;
00715 double err = 1e-6;
00716 double loaded_x, loaded_y;
00717 if (!initialized) {
00718 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00719 } else {
00720 ifstream input(filename.c_str());
00721 if (input != NULL && !input.eof()) {
00722
00723 input >> inBuf;
00724
00725 while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
00726 input >> inBuf;
00727 }
00728
00729 input.ignore(9999, '\n');
00730 for (x = 0; x < size_x; x++) {
00731 for (y = 0; y < size_y; y++) {
00732 input >> loaded_x >> loaded_y;
00733
00734 if (fabs(loaded_x - grid_x[x][y]) > err || fabs(loaded_y - grid_y[x][y]) > err) {
00735 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]);
00736 } else {
00737 input >> matrix_array[x][y];
00738 }
00739 }
00740 }
00741 } else {
00742 throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
00743
00744
00745
00746 }
00747 input.close();
00748 }
00749 }
00750
00751
00753
00754
00755
00757
00762 template<class T>
00763 Matrix3D<T>::Matrix3D( int x_size, int y_size, int z_size ) {
00764 initialized = false;
00765
00766 AllocateMemory(x_size, y_size, z_size);
00767 }
00768
00775 template<class T>
00776 Matrix3D<T>::Matrix3D( const Matrix3D<T> &M ) {
00777 initialized = false;
00778 this->operator = (M);
00779 }
00780
00784 template<class T>
00785 Matrix3D<T>::~Matrix3D() {
00786 if (initialized) free_matrix<T>(matrix_array, size_x, size_y);
00787
00788 }
00789
00795 template<class T>
00796 void Matrix3D<T>::AllocateMemory( int x_size, int y_size, int z_size) {
00797 this->size_x = x_size;
00798 this->size_y = y_size;
00799 this->size_z = z_size;
00800 matrix_array = matrix<T>(x_size, y_size, z_size);
00801 plane_array = matrix_array[0][0];
00802 initialized = true;
00803
00804 for (int i = 0; i < size_x; i++)
00805 for (int j = 0; j < size_y; j++)
00806 for (int k = 0; k < size_z; k++)
00807 matrix_array[i][j][k] = 0;
00808 }
00809
00813 template<class T>
00814 int Matrix3D<T>::index1d(int x, int y, int z) {
00815 return (x*size_y + y)*size_z + z;
00816 }
00817
00824 template<class T>
00825 T** Matrix3D<T>::operator[](int i) {
00826 #ifdef DEBUG_MODE
00827 if (!initialized) {
00828 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00829 }
00830 #endif
00831
00832 return matrix_array[i];
00833 }
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00850 template<class T>
00851 T& Matrix3D<T>::operator()(int x, int y, int z) {
00852 #ifdef DEBUG_MODE
00853 if (!initialized) {
00854 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
00855 }
00856 #endif
00857
00858 return plane_array[(x*size_y + y)*size_z + z];
00859 }
00860
00866 template<class T>
00867 Matrix3D<T>& Matrix3D<T>::operator= (const Matrix3D<T> &M) {
00868
00869 if (this!=&M) {
00870
00871 if (M.initialized) {
00872
00873 if (initialized) free_matrix<T>(matrix_array, this->size_x, this->size_y);
00874 this->size_x = M.size_x;
00875 this->size_y = M.size_y;
00876 this->size_z = M.size_z;
00877 this->name = M.name;
00878
00879 this->AllocateMemory(this->size_x, this->size_y, this->size_z);
00880
00881
00882
00883
00884 memcpy( matrix_array[0][0], M.matrix_array[0][0], this->size_x * this->size_y * this->size_z * sizeof( T ) );
00885 } else {
00886 this->initialized = false;
00887 }
00888 }
00889 return *this;
00890 }
00891
00896 template<class T>
00897 Matrix3D<T>& Matrix3D<T>::operator= (Matrix2D<T> &M) {
00898 if (M.initialized) {
00899 if (initialized) free_matrix<T>(matrix_array, this->size_x, this->size_y);
00900 this->size_x = 1;
00901 this->size_y = M.size_x;
00902 this->size_z = M.size_y;
00903
00904 this->AllocateMemory(this->size_x, this->size_y, this->size_z);
00905 for (int y = 0; y < this->size_y; y++)
00906
00907 for (int z = 0; z < this->size_z; z++)
00908 this->matrix_array[0][y][z] = M[y][z];
00909 } else {
00910 this->initialized = false;
00911 }
00912 return *this;
00913 }
00914
00918 template<class T>
00919 Matrix3D<T>& Matrix3D<T>::operator= (T Val) {
00920 for (int x = 0; x < this->size_x; x++)
00921 for (int y = 0; y < this->size_y; y++)
00922 for (int z = 0; z < this->size_z; z++)
00923 this->matrix_array[x][y][z] = Val;
00924 return *this;
00925 }
00926
00930 template<class T>
00931 Matrix3D<T> Matrix3D<T>::operator/ (T Val) {
00932 int x, y, z;
00933 Matrix3D<T> Tmp(*this);
00934 for (x = 0; x < this->size_x; x++)
00935 for (y = 0; y < this->size_y; y++)
00936 for (z = 0; z < this->size_z; z++)
00937 Tmp[x][y][z] = matrix_array[x][y][z] / Val;
00938 return Tmp;
00939 }
00940
00944 template<class T>
00945 Matrix3D<T> Matrix3D<T>::operator* (T Val) {
00946 int x, y, z;
00947 Matrix3D<T> Tmp(*this);
00948 for (x = 0; x < this->size_x; x++)
00949 for (y = 0; y < this->size_y; y++)
00950 for (z = 0; z < this->size_z; z++)
00951 Tmp[x][y][z] = matrix_array[x][y][z] * Val;
00952 return Tmp;
00953 }
00954
00958 template<class T>
00959 Matrix3D<T> Matrix3D<T>::operator/ (Matrix3D<T> &M) {
00960 int x, y, z;
00961 Matrix3D<T> Tmp(*this);
00962 for (x = 0; x < this->size_x; x++)
00963 for (y = 0; y < this->size_y; y++)
00964 for (z = 0; z < this->size_z; z++)
00965 Tmp[x][y][z] = matrix_array[x][y][z] / M[x][y][z];
00966 return Tmp;
00967 }
00968
00972 template<class T>
00973 Matrix3D<T> Matrix3D<T>::operator* (Matrix3D<T> &M) {
00974 int x, y, z;
00975 Matrix3D<T> Tmp(*this);
00976 for (x = 0; x < this->size_x; x++)
00977 for (y = 0; y < this->size_y; y++)
00978 for (z = 0; z < this->size_z; z++)
00979 Tmp[x][y][z] = matrix_array[x][y][z] * M[x][y][z];
00980 return Tmp;
00981 }
00982
00986 template<class T>
00987 Matrix3D<T> Matrix3D<T>::operator+ (Matrix3D<T> &M) {
00988 int x, y, z;
00989 Matrix3D<T> Tmp(*this);
00990 for (x = 0; x < this->size_x; x++)
00991 for (y = 0; y < this->size_y; y++)
00992 for (z = 0; z < this->size_z; z++)
00993 Tmp[x][y][z] = matrix_array[x][y][z] + M[x][y][z];
00994 return Tmp;
00995 }
00996
01000 template<class T>
01001 Matrix3D<T> Matrix3D<T>::operator- (Matrix3D<T> &M) {
01002 int x, y, z;
01003 Matrix3D<T> Tmp(*this);
01004 for (x = 0; x < this->size_x; x++)
01005 for (y = 0; y < this->size_y; y++)
01006 for (z = 0; z < this->size_z; z++)
01007 Tmp[x][y][z] = matrix_array[x][y][z] - M[x][y][z];
01008 return Tmp;
01009 }
01010
01015 template<class T>
01016 void Matrix3D<T>::writeToFile(string filename) {
01017 int x, y, z;
01018 ofstream output(filename.c_str());
01019 output << "VARIABLES = \""<< ((this->name!="")?this->name:"f") <<"\" "<< endl;
01020 output << "ZONE T=\"" << filename << "\", I=" << size_z << ", J=" << size_y << ", K=" << size_x << endl;
01021 for (x = 0; x < size_x; x++) {
01022 for (y = 0; y < size_y; y++) {
01023 for (z = 0; z < size_z; z++) {
01024 output << matrix_array[x][y][z] << endl;
01025 }
01026 }
01027 }
01028 output.close();
01029 }
01030
01035 template<class T>
01036 void Matrix3D<T>::writeToFile(string filename, Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z) {
01037 int x, y, z;
01038 ofstream output(filename.c_str());
01039 if (output==NULL && (filename.find("Debug") == string::npos)) {
01040 throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
01041 }
01042 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;
01043 output << "ZONE T=\"" << filename << "\", I=" << size_z << ", J=" << size_y << ", K=" << size_x << endl;
01044 for (x = 0; x < size_x; x++) {
01045 for (y = 0; y < size_y; y++) {
01046 for (z = 0; z < size_z; z++) {
01047 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;
01048 }
01049 }
01050 }
01051 output.close();
01052 }
01053
01057 template<class T>
01058 void Matrix3D<T>::readFromFile(string filename) {
01059 int x, y, z;
01060 string inBuf;
01061 if (!initialized) {
01062 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
01063 } else {
01064 ifstream input(filename.c_str());
01065 if (input != NULL && !input.eof()) {
01066
01067 input >> inBuf;
01068
01069 while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
01070
01071 input >> inBuf;
01072 }
01073
01074 input.ignore(9999, '\n');
01075 for (x = 0; x < size_x; x++) {
01076 for (y = 0; y < size_y; y++) {
01077 for (z = 0; z < size_z; z++) {
01078 input >> matrix_array[x][y][z];
01079 }
01080 }
01081 }
01082 } else {
01083 throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
01084
01085
01086
01087 }
01088 input.close();
01089 }
01090 }
01091
01095 template<class T>
01096 void Matrix3D<T>::readFromFile(string filename, Matrix3D<T> &grid_x, Matrix3D<T> &grid_y, Matrix3D<T> &grid_z) {
01097 int x, y, z;
01098 string inBuf;
01099 if (!initialized) {
01100 throw error_msg("MATRIX_ERROR", "Using unitialized matrix");
01101 } else {
01102 ifstream input(filename.c_str());
01103 if (input != NULL && !input.eof()) {
01104
01105 input >> inBuf;
01106
01107
01108 while (strcmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) {
01109 input >> inBuf;
01110 }
01111
01112 input.ignore(9999, '\n');
01113 for (x = 0; x < size_x; x++) {
01114 for (y = 0; y < size_y; y++) {
01115 for (z = 0; z < size_z; z++) {
01116 input >> grid_x[x][y][z] >> grid_y[x][y][z] >> grid_z[x][y][z] >> matrix_array[x][y][z];
01117
01118 input.ignore(9999, '\n');
01119
01120
01121
01122
01123
01124
01125
01126
01127 }
01128 }
01129 }
01130 } else {
01131 throw error_msg("MATRIX_LOAD_ERROR", "Error reading file %s.\n", filename.c_str());
01132
01133
01134
01135 }
01136 input.close();
01137 }
01138 }
01139
01140
01144 template<class T>
01145 T Matrix3D<T>::max() {
01146 T tmp = 0;
01147 int x, y, z;
01148 for (x = 0; x < size_x; x++) {
01149 for (y = 0; y < size_y; y++) {
01150 for (z = 0; z < size_z; z++) {
01151 tmp = (tmp>matrix_array[x][y][z])?tmp:matrix_array[x][y][z];
01152 }
01153 }
01154 }
01155 return tmp;
01156 }
01157
01161 template<class T>
01162 T Matrix3D<T>::maxabs() {
01163 T tmp = 0;
01164 int x, y, z;
01165 for (x = 0; x < size_x; x++) {
01166 for (y = 0; y < size_y; y++) {
01167 for (z = 0; z < size_z; z++) {
01168 tmp = (tmp>fabs((double)matrix_array[x][y][z]))?tmp:fabs((double)matrix_array[x][y][z]);
01169 }
01170 }
01171 }
01172 return tmp;
01173 }
01174
01178 template<class T>
01179 Matrix3D<T> Matrix3D<T>::abs() {
01180 Matrix3D<T> tmp(this->size_x, this->size_y, this->size_z);
01181 int x, y, z;
01182 for (x = 0; x < size_x; x++) {
01183 for (y = 0; y < size_y; y++) {
01184 for (z = 0; z < size_z; z++) {
01185 tmp[x][y][z] = (matrix_array[x][y][z]>0)?matrix_array[x][y][z]:-matrix_array[x][y][z];
01186 }
01187 }
01188 }
01189 return tmp;
01190 }
01191
01192
01196 template<class T>
01197 Matrix2D<T> Matrix3D<T>::xSlice(int p_x) {
01198 int y, z;
01199 Matrix2D<T> tmp(this->size_y, this->size_z);
01200 tmp.name = this->name + "_slice";
01201 for (y = 0; y < size_y; y++) {
01202 for (z = 0; z < size_z; z++) {
01203 tmp[y][z] = matrix_array[p_x][y][z];
01204 }
01205 }
01206 return tmp;
01207 }
01208
01212 template<class T>
01213 Matrix2D<T> Matrix3D<T>::ySlice(int p_y) {
01214 int x, z;
01215 Matrix2D<T> tmp(this->size_x, this->size_z);
01216 tmp.name = this->name + "_slice";
01217 for (x = 0; x < size_x; x++) {
01218 for (z = 0; z < size_z; z++) {
01219 tmp[x][z] = matrix_array[x][p_y][z];
01220 }
01221 }
01222 return tmp;
01223 }
01224
01228 template<class T>
01229 Matrix2D<T> Matrix3D<T>::zSlice(int p_z) {
01230 int x, y;
01231 Matrix2D<T> tmp(this->size_x, this->size_y);
01232 tmp.name = this->name + "_slice";
01233 for (x = 0; x < size_x; x++) {
01234 for (y = 0; y < size_y; y++) {
01235 tmp[x][y] = matrix_array[x][y][p_z];
01236 }
01237 }
01238 return tmp;
01239 }
01240
01241
01243
01244
01245
01247
01252 template<class T>
01253 void Matrix1D<T>::Interpolate(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
01254
01255
01256
01257 if (this->size_x != new_grid.size_x) {
01258 throw error_msg("INTERPOLATION_ERROR", "New grid and function size mismatch!\n");
01259
01260
01261 }
01262
01263
01264
01265
01266 Maths::Interpolation::Linear LinInterp(old_grid.size_x, old_grid.matrix_array, old_function.matrix_array);
01267
01268 int i;
01269 for (i = 0; i < new_grid.size_x; i++)
01270
01271 matrix_array[i] = LinInterp.getValue(new_grid.matrix_array[i]);
01272
01273
01274
01275 }
01276
01293 template<class T>
01294 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) {
01295
01296
01297
01298 if (this->size_x != new_grid.size_x) {
01299 throw error_msg("SPLINE_INTERPOLATION", "New grid and function size mismatch!");
01300 }
01301
01302 if (this->size_x == 1) {
01303 this->matrix_array[0] = old_function.matrix_array[0];
01304 return;
01305 }
01306
01307
01308 int i=0, i_min=0, i_max=old_grid.size_x-1, interp_size=0;
01309
01310
01311
01312
01313
01314
01315
01316 i_min = 1; i_max = old_grid.size_x-2;
01317 interp_size = i_max-i_min+1;
01318
01319 double yp0, ypn;
01320
01321
01322
01323
01324
01325 yp0 = 1e99;
01326
01327 ypn = 1e99;
01328
01329 Matrix1D<T> y2(old_grid.size_x);
01330
01331
01332
01333 spline(&old_grid.matrix_array[i_min], &old_function.matrix_array[i_min], interp_size, yp0, ypn, &y2.matrix_array[i_min]);
01334
01335
01336 for (i = 0; i < old_grid.size_x; i++) {
01337
01338 if (fabs(y2[i]) >= max_second_der) {
01339
01340 if (y2[i] > 0) {
01341
01342
01343 y2[i] = max_second_der + (y2[i] - max_second_der)*(1.0 - lin_spline_coef);
01344 } else {
01345 y2[i] = - max_second_der + (y2[i] + max_second_der)*(1.0 - lin_spline_coef);
01346 }
01347 }
01348 }
01349
01350
01351 for (i = 0; i < new_grid.size_x; i++) {
01352 if (i == 0) {
01353 this->matrix_array[i] = lb;
01354 } else if (i == new_grid.size_x-1) {
01355 this->matrix_array[i] = ub;
01356 } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[0] ) {
01357 this->matrix_array[i] = lb;
01358
01359 } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-1]) {
01360 this->matrix_array[i] = ub;
01361
01362 } else if (new_grid.matrix_array[i] <= old_grid.matrix_array[1] ) {
01363
01364 double a = (new_grid.matrix_array[i] - old_grid.matrix_array[0])/(old_grid.matrix_array[1] - old_grid.matrix_array[0]);
01365 this->matrix_array[i] = old_function.matrix_array[0] + a*(old_function.matrix_array[1] - old_function.matrix_array[0]);
01366 } else if (new_grid.matrix_array[i] >= old_grid.matrix_array[old_grid.size_x-2]) {
01367
01368
01369
01370 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]);
01371 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]);
01372 } else if (i <= 1) {
01373 int k = 1;
01374
01375
01376 while(old_grid.matrix_array[k] < new_grid.matrix_array[i]) {k++;}
01377 double a = (new_grid.matrix_array[i] - old_grid.matrix_array[k-1])/(old_grid.matrix_array[k] - old_grid.matrix_array[k-1]);
01378 this->matrix_array[i] = old_function.matrix_array[k-1] + a*(old_function.matrix_array[k] - old_function.matrix_array[k-1]);
01379 } else {
01380
01381
01382
01384 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]);
01385 }
01386 }
01387
01388 }
01389
01394 template<class T>
01395 void Matrix1D<T>::Polilinear(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid, double lb, double ub) {
01396
01397
01398 if (this->size_x != new_grid.size_x) {
01399 throw error_msg("POLINOMIAL_INTERPOLATION", "New grid and function size mismatch!");
01400 }
01401
01402 int i;
01403 for (i = 0; i < new_grid.size_x; i++)
01404 matrix_array[i] = polilinear( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], lb, ub);
01405
01406 }
01407
01412 template<class T>
01413 void Matrix1D<T>::Polint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
01414
01415
01416 if (this->size_x != new_grid.size_x) {
01417 throw error_msg("POLINT_INTERPOLATION", "New grid and function size mismatch!");
01418 }
01419
01420 int i;
01421 double err;
01422 for (i = 0; i < new_grid.size_x; i++)
01423 polint( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i], &err);
01424
01425 }
01426
01427
01428 template<class T>
01429 void Matrix1D<T>::Ratint(Matrix1D<T> &old_function, Matrix1D<T> &old_grid, Matrix1D<T> &new_grid) {
01430
01431
01432 if (this->size_x != new_grid.size_x) {
01433 throw error_msg("RATIONAL_INTERPOLATION", "New grid and function size mismatch!");
01434 }
01435
01436 int i;
01437 double err;
01438 for (i = 0; i < new_grid.size_x; i++)
01439 ratint( old_grid.matrix_array, old_function.matrix_array, old_grid.size_x, new_grid.matrix_array[i], &this->matrix_array[i], &err);
01440
01441 }
01442
01447 template<class T>
01448 double Linear2D(double x, double y, T &old_grid_x, T &old_grid_y, T &f) {
01449
01450 int i = 0, j = 0;
01451
01452 while (x >= old_grid_x[++i][0] && i < old_grid_y.size_y-1);
01453 while (y >= old_grid_y[0][++j] && j < old_grid_x.size_x-1);
01454
01455 double t = (x - old_grid_x[i-1][0])/(old_grid_x[i][0] - old_grid_x[i-1][0]);
01456 double u = (y - old_grid_y[0][j-1])/(old_grid_y[0][j] - old_grid_y[0][j-1]);
01457
01458 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];
01459 }
01460
01461
01467 template<class T>
01468 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) {
01469
01470 Matrix2D<T> res(new_grid_x.size_x, new_grid_y.size_y);
01471
01472 if (old_function.size_x != old_grid_x.size_x || old_function.size_y != old_grid_y.size_y) {
01473 throw error_msg("LINEAR_INTERPOLATION", "Old grid and function size mismatch!");
01474 }
01475
01476 if (this->size_x != new_grid_x.size_x || this->size_y != new_grid_y.size_y) {
01477 throw error_msg("LINEAR_INTERPOLATION", "New grid and function size mismatch!");
01478 }
01479
01480 int i, j;
01481 for (i = 0; i < new_grid_x.size_x; i++) {
01482 for (j = 0; j < new_grid_y.size_y; j++) {
01483 this->matrix_array[i][j] = Linear2D(new_grid_x[i][j], new_grid_y[i][j], old_grid_x, old_grid_y, old_function);
01484 }
01485 }
01486
01487 return;
01488
01489 }
01490
01491
01499 CalculationMatrix::CalculationMatrix(int x_size, int y_size, int z_size, int n_of_diags) {
01500 this->initialized = false;
01501 Initialize(x_size, y_size, z_size, n_of_diags);
01502 }
01503
01507 void CalculationMatrix::Initialize(int x_size, int y_size, int z_size, int n_of_diags) {
01508 this->initialized = false;
01509 this->x_size = x_size;
01510 this->y_size = y_size;
01511 this->total_size = x_size;
01512 if (y_size > 0) this->total_size = this->total_size * y_size;
01513 if (z_size > 0) this->total_size = this->total_size * z_size;
01514
01515
01516
01517 int x, y, z, id;
01518
01519 if (z_size > 1) {
01520 for (x = -n_of_diags; x <= n_of_diags; x++) {
01521 for (y = -n_of_diags; y <= n_of_diags; y++) {
01522 for (z = -n_of_diags; z <= n_of_diags; z++) {
01523
01524 id = this->index1d(x, y, z);
01525
01526 (*this)[id] = Matrix1D<double>(this->total_size);
01527 (*this)[id] = 0;
01528 }
01529 }
01530 }
01531
01532 } else if (y_size > 1) {
01533 for (x = -n_of_diags; x <= n_of_diags; x++) {
01534 for (y = -n_of_diags; y <= n_of_diags; y++) {
01535
01536 id = this->index1d(x, y);
01537
01538 (*this)[id] = Matrix1D<double>(this->total_size);
01539 (*this)[id] = 0;
01540 }
01541 }
01542
01543 } else if (x_size > 1) {
01544 for (x = -n_of_diags; x <= n_of_diags; x++) {
01545
01546 id = this->index1d(x);
01547
01548 (*this)[id] = Matrix1D<double>(this->total_size);
01549 (*this)[id] = 0;
01550 }
01551 }
01552
01553
01554 if (y_size == 0) this->x_size = 0;
01555
01556 if (z_size == 0) this->y_size = 0;
01557
01558 this->initialized = true;
01559 }
01560
01564 int CalculationMatrix::index1d(int x, int y, int z) {
01565 return (z*y_size + y)*x_size + x;
01566 }
01567
01571 void CalculationMatrix::writeToFile(string filename) {
01572 int in;
01573
01574 ofstream output(filename.c_str());
01575 if (output==NULL && (filename.find("Debug") == string::npos)) {
01576 throw error_msg("FILE", "Unable to output file: %s", filename.c_str());
01577 }
01578 output << "VARIABLES = \"";
01579 for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
01580 output << "\"" << it->first << "\", ";
01581 }
01582 output << endl;
01583 output << "ZONE T=\"" << "\", I=" << total_size << endl;
01584 for (in = 0; in < this->total_size; in++) {
01585 for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
01586 output << "\t" << it->second[in];
01587 }
01588 output << endl;
01589 }
01590 output.close();
01591 }
01592
01593
01595
01597
01598 template class Matrix1D<double>;
01599
01600
01601 template class Matrix2D<double>;
01602
01603 template class Matrix3D<double>;
01604 template class Matrix3D<char>;
01605
01606
01607 #endif