6 #ifndef matrix_array_MATRIX_CPP
7 #define matrix_array_MATRIX_CPP
13 #if defined(_WIN32) || defined(_WIN64)
14 #define strncasecmp _strnicmp
15 #define strcasecmp _stricmp
21 template<
class T>
inline T*
matrix(
long Rows)
26 throw error_msg(
"MEMORY_ERROR",
"Memory can't be initialized: %s size", Rows*
sizeof(
T));
32 template<
class T>
inline T**
matrix(
long Rows,
long Columns)
38 throw error_msg(
"MEMORY_ERROR",
"Memory can't be initialized: %s size", Rows *
sizeof(
T));
41 m[0] =
new T[Rows * Columns];
44 throw error_msg(
"MEMORY_ERROR",
"Memory can't be initialized: %s size", Rows * Columns *
sizeof(
T));
47 for(
long i=1; i<Rows; i++) m[i] = m[i-1] + Columns;
52 template<
class T>
inline T***
matrix(
int size_x,
int size_y,
int size_z)
55 T ***
m=
new T**[size_x];
58 throw error_msg(
"MEMORY_ERROR",
"Memory can't be initialized: %s size", size_x *
sizeof(
T));
60 for (
int x = 0; x < size_x; x++) {
62 m[x] =
new T*[size_y];
65 throw error_msg(
"MEMORY_ERROR",
"Memory can't be initialized: %s size", size_y *
sizeof(
T));
69 m[0][0] =
new T[size_x * size_y * size_z];
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));
74 for (
int x = 0; x < size_x; x++) {
75 for (
int y = 0; y < size_y; y++) {
77 m[x][y] = m[0][0] + (x*size_y + y)*size_z;
99 template<
class T>
inline void free_matrix(
T***
m,
int size_x,
int size_y) {
101 for (
int x = 0; x < size_x; x++) {
127 AllocateMemory(x_size);
142 AllocateMemory(x_size);
156 this->operator = (M);
164 if (initialized) free_matrix<T>(matrix_array);
174 this->size_x = x_size;
176 matrix_array = matrix<T>(x_size);
192 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
196 return matrix_array[i];
209 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
213 return matrix_array[i];
227 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
231 return matrix_array[x];
247 if (initialized) free_matrix<T>(matrix_array);
252 this->AllocateMemory(this->size_x);
254 memcpy( matrix_array, M.
matrix_array, this->size_x *
sizeof(
T ) );
256 this->initialized =
false;
270 for (
int x = 0; x < this->size_x; x++)
271 matrix_array[x] = Val;
284 for (
int x = 0; x < this->size_x; x++)
285 Tmp[x] = matrix_array[x] * Val;
299 for (x = 0; x < this->size_x; x++)
300 Tmp[x] = matrix_array[x] / Val;
313 for (x = 0; x < this->size_x; x++)
327 for (x = 0; x < this->size_x; x++)
339 for (x = 0; x < this->size_x; x++) {
340 res += matrix_array[x] * matrix_array[x];
350 if (this->size_x != W.
size_x)
351 throw error_msg(
"DOT_PRODUCT",
"Size is different");
355 for (x = 0; x < this->size_x; x++) {
356 res += matrix_array[x] * W[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());
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;
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());
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;
409 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
411 ifstream input(filename.c_str());
412 if (input != NULL && !input.eof()) {
416 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
420 input.ignore(9999,
'\n');
421 for (x = 0; x < size_x; x++) {
422 input >> matrix_array[x];
425 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
444 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
446 ifstream input(filename.c_str());
447 if (input != NULL && !input.eof()) {
452 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
456 input.ignore(9999,
'\n');
457 for (x = 0; x < size_x; x++) {
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]);
463 input >> matrix_array[x];
467 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
494 AllocateMemory(x_size, y_size);
506 this->operator = (M);
514 if (initialized) free_matrix<T>(matrix_array);
525 this->size_x = x_size;
526 this->size_y = y_size;
528 matrix_array = matrix<T>(x_size, y_size);
547 if (initialized && (size_x != M.
size_x || size_y != M.
size_y)) {
548 free_matrix<T>(matrix_array);
555 this->AllocateMemory(this->size_x, this->size_y);
561 memcpy(matrix_array[0], M.matrix_array[0], this->size_x * this->size_y *
sizeof(
T) );
563 this->initialized =
false;
578 for (x = 0; x < this->size_x; x++)
579 for (y = 0; y < this->size_y; y++)
580 matrix_array[x][y] = val;
582 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
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;
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;
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];
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];
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;
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());
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;
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;
728 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
730 ifstream input(filename.c_str());
731 if (input != NULL && !input.eof()) {
735 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
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];
746 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
763 double loaded_x, loaded_y;
765 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
767 ifstream input(filename.c_str());
768 if (input != NULL && !input.eof()) {
772 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
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;
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]);
784 input >> matrix_array[x][y];
789 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
813 AllocateMemory(x_size, y_size, z_size);
825 this->operator = (M);
833 if (initialized) free_matrix<T>(matrix_array, size_x, size_y);
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];
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;
867 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
871 return matrix_array[i];
883 throw error_msg(
"MATRIX_ERROR",
"Using not initialized matrix");
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");
890 return plane_array[(x*size_y + y)*size_z + z];
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);
914 this->AllocateMemory(this->size_x, this->size_y, this->size_z);
921 memcpy( matrix_array[0][0], M.matrix_array[0][0], this->size_x * this->size_y * this->size_z *
sizeof(
T ) );
923 this->initialized =
false;
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;
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];
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];
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;
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;
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;
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;
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];
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];
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];
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];
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;
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;
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];
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];
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;
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());
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;
1212 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
1214 ifstream input(filename.c_str());
1215 if (input != NULL && !input.eof()) {
1219 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
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];
1233 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
1250 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
1252 ifstream input(filename.c_str());
1253 if (input != NULL && !input.eof()) {
1258 while (strcmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
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];
1268 input.ignore(9999,
'\n');
1281 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
1297 return (x*size_y + y)*size_z + 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];
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]);
1340 Matrix3D<T> tmp(this->size_x, this->size_y, this->size_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];
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];
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];
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];
1418 if (this->size_x != new_grid.
size_x) {
1419 throw error_msg(
"INTERPOLATION_ERROR",
"New grid and function size mismatch!\n");
1430 for (i = 0; i < new_grid.
size_x; i++)
1459 if (this->size_x != new_grid.
size_x) {
1460 throw error_msg(
"SPLINE_INTERPOLATION",
"New grid and function size mismatch!");
1463 if (this->size_x == 1) {
1469 int i=0, i_min=0, i_max=old_grid.
size_x-1, interp_size=0;
1477 i_min = 1; i_max = old_grid.
size_x-2;
1478 interp_size = i_max-i_min+1;
1497 for (i = 0; i < old_grid.
size_x; i++) {
1499 if (fabs(y2[i]) >= max_second_der) {
1503 y2[i] = max_second_der + (y2[i] - max_second_der)*(1.0 - lin_spline_coef);
1505 y2[i] = - max_second_der + (y2[i] + max_second_der)*(1.0 - lin_spline_coef);
1512 for (i = 0; i < new_grid.
size_x; i++) {
1514 this->matrix_array[i] = lb;
1515 }
else if (i == new_grid.
size_x-1) {
1516 this->matrix_array[i] = ub;
1518 this->matrix_array[i] = lb;
1521 this->matrix_array[i] = ub;
1625 if (this->size_x != new_grid.
size_x) {
1626 throw error_msg(
"SPLINE_INTERPOLATION",
"New grid and function size mismatch!");
1629 if (this->size_x == 1) {
1635 int i=0, i_min=0, i_max=old_grid.
size_x-1, interp_size=0;
1643 i_min = 0; i_max = old_grid.
size_x-1;
1644 interp_size = i_max-i_min+1;
1661 for (i = 0; i < new_grid.
size_x; i++) {
1676 if (this->size_x != new_grid.
size_x) {
1677 throw error_msg(
"POLINOMIAL_INTERPOLATION",
"New grid and function size mismatch!");
1681 for (i = 0; i < new_grid.
size_x; i++)
1694 if (this->size_x != new_grid.
size_x) {
1695 throw error_msg(
"POLINT_INTERPOLATION",
"New grid and function size mismatch!");
1700 for (i = 0; i < new_grid.
size_x; i++)
1710 if (this->size_x != new_grid.
size_x) {
1711 throw error_msg(
"RATIONAL_INTERPOLATION",
"New grid and function size mismatch!");
1716 for (i = 0; i < new_grid.
size_x; i++)
1859 matrix_array[0] = lb;
1860 matrix_array[size_x-1] = ub;
1877 double Linear2D(
double x,
double y,
T &old_grid_x,
T &old_grid_y,
T &f) {
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);
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]);
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];
1902 throw error_msg(
"LINEAR_INTERPOLATION",
"Old grid and function size mismatch!");
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!");
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);
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);
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;
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++) {
1953 id = this->index1d(x, y, z);
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++) {
1965 id = this->index1d(x, y);
1972 }
else if (x_size > 1) {
1973 for (x = -n_of_diags; x <= n_of_diags; x++) {
1975 id = this->index1d(x);
1983 if (y_size == 0) this->x_size = 0;
1985 if (z_size == 0) this->y_size = 0;
1987 this->initialized =
true;
1994 return (z*y_size + y)*x_size + x;
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());
2007 output <<
"VARIABLES = ";
2008 for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2009 output <<
"\"" << it->first <<
"\", ";
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];
void readFromFile(string filename)
bool initialized
Flag, equal true if initialized.
void Initialize(int x_size, int y_size=1, int z_size=1, int n_of_diags=1)
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)
Matrix1D times(const Matrix1D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
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)
const Matrix3D & operator+() const
Return itself as positive version of values.
Matrix3D()
Default constructor. Do nothing.
const Matrix3D operator-() const
Return negative version of values.
Matrix3D operator*(const T Val) const
string name
name of the Matrix
Matrix2D< T > ySlice(int p_y) const
bool initialized
Flag, equal true if initialized.
T dot(const Matrix1D< T > &M) const
Dot product.
Matrix1D & operator=(const Matrix1D< T > &M)
void writeToFile(string filename)
Save matrix to a file.
Matrix1D operator*(const T Val) const
Matrix2D times(const Matrix2D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
void AllocateMemory(int size_x, int size_y)
Matrix3D operator/(const T Val) const
void Ratint(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
rational
Matrix3D & times_equal(const Matrix3D< T > &M)
Arraywise multiplication (A.*B), stores result in the matrix it's applied to.
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)
Linear interpolation of two arrays.
Matrix3D & divide_equal(const Matrix3D< T > &M)
Arraywise division (A./B), stores result in the matrix it's applied to.
Matrix3D divide(const Matrix3D< T > &M) const
Arraywise division (A./B), stores result in a new matrix.
void AllocateMemory(int x_size)
bool initialized
Flag, equal true if initialized.
int index1d(int x, int y=0, int z=0)
void splint(double *xa, double *ya, double *y2a, int n, double x, double *y)
void AllocateMemory(int size_x, int size_y, int size_z)
string name
name of the Matrix
int index1d(int x, int y) const
one dimensional matrix class
void Polint(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
void polint(double *xa, double *ya, int n, double x, double *y, double *dy)
Matrix3D times(const Matrix3D< T > &M) const
Arraywise multiplication (A.*B), stores result in a new matrix.
Matrix2D operator/(const T Val) const
T & operator()(int x, int y, int z)
Return the (x,y,z) value of matrix.
Matrix3D & operator+=(const Matrix3D< T > &M)
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]...
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)
int size_y
size x, size y, size z
T & operator[](int i)
Return the i-th value of matrix.
Matrix3D & operator=(const Matrix3D< T > &M)
void writeToFile(string filename)
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)
Matrix2D & operator=(const Matrix2D< T > &M)
Matrix1D divide(const Matrix1D< T > &M) const
Arraywise division (A./B), stores result in a new matrix.
T * matrix_array
Array to keep the values.
double T(double alpha)
Function for bounce time.
void Interpolate(Matrix1D< T > &old_function, Matrix1D< T > &old_grid, Matrix1D< T > &new_grid)
Matrix2D< T > xSlice(int p_x) const
double Linear2D(double x, double y, T &old_grid_x, T &old_grid_y, T &f)
three dimensional matrix class
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
void ratint(double *xa, double *ya, int n, double x, double *y, double *dy)
Matrix2D divide(const Matrix2D< T > &M) const
Arraywise division (A./B), stores result in a new matrix.
two dimensional matrix class
Error message - stack of single_errors.
double getValue(double x)
Matrix3D & operator-=(const Matrix3D< T > &M)
void readFromFile(string filename)
Load matrix from a file.
T * matrix(long Rows)
Allocating memory for 1D matrix.
void readFromFile(string filename)
string name
name of the Matrix
void writeToFile(string filename)
Matrix2D operator*(const T Val) const
void writeToFile(string filename)
int index1d(int x, int y, int z)
Returns index of the element (x,y,z) in 1d array.
double polilinear(double *xa, double *ya, int n, double x, double lb, double ub)
int size_x
size x, size y, size z
T & operator()(int x)
Return the x-th value of matrix.
Matrix2D< T > zSlice(int p_z) const
Matrix1D operator/(const T Val) const
Matrix3D & operator*=(const T Val)
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)
Matrix3D & operator/=(const T Val)
int size_z
size x, size y, size z