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
22 template<
class T>
inline T*
matrix(
long Rows)
27 throw error_msg(
"MEMORY_ERRROR",
"Memory can't be initialized: %s size", Rows*
sizeof(
T));
34 template<
class T>
inline T**
matrix(
long Rows,
long Columns)
40 throw error_msg(
"MEMORY_ERRROR",
"Memory can't be initialized: %s size", Rows *
sizeof(
T));
43 m[0] =
new T[Rows * Columns];
46 throw error_msg(
"MEMORY_ERRROR",
"Memory can't be initialized: %s size", Rows * Columns *
sizeof(
T));
49 for(
long i=1; i<Rows; i++) m[i] = m[i-1] + Columns;
55 template<
class T>
inline T***
matrix(
int size_x,
int size_y,
int size_z)
58 T ***
m=
new T**[size_x];
61 throw error_msg(
"MEMORY_ERRROR",
"Memory can't be initialized: %s size", size_x *
sizeof(
T));
63 for (
int x = 0; x < size_x; x++) {
65 m[x] =
new T*[size_y];
68 throw error_msg(
"MEMORY_ERRROR",
"Memory can't be initialized: %s size", size_y *
sizeof(
T));
72 m[0][0] =
new T[size_x * size_y * size_z];
74 if (m[0][0] == NULL) {
75 throw error_msg(
"MEMORY_ERRROR",
"Memory can't be initialized: %s size", size_x * size_y * size_z *
sizeof(
T));
77 for (
int x = 0; x < size_x; x++) {
78 for (
int y = 0; y < size_y; y++) {
80 m[x][y] = m[0][0] + (x*size_y + y)*size_z;
104 for (
int x = 0; x < size_x; x++) {
130 AllocateMemory(x_size);
145 AllocateMemory(x_size);
159 this->operator = (M);
167 if (initialized) free_matrix<T>(matrix_array);
177 this->size_x = x_size;
179 matrix_array = matrix<T>(x_size);
195 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
199 return matrix_array[i];
212 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
216 return matrix_array[i];
230 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
234 return matrix_array[x];
250 if (initialized) free_matrix<T>(matrix_array);
255 this->AllocateMemory(this->size_x);
257 memcpy( matrix_array, M.
matrix_array, this->size_x *
sizeof(
T ) );
259 this->initialized =
false;
273 for (
int x = 0; x < this->size_x; x++)
274 matrix_array[x] = Val;
287 for (
int x = 0; x < this->size_x; x++)
288 Tmp[x] = matrix_array[x] * Val;
302 for (x = 0; x < this->size_x; x++)
303 Tmp[x] = matrix_array[x] / Val;
316 for (x = 0; x < this->size_x; x++)
330 for (x = 0; x < this->size_x; x++)
342 for (x = 0; x < this->size_x; x++) {
343 res += matrix_array[x] * matrix_array[x];
353 if (this->size_x != W.
size_x)
354 throw error_msg(
"DOT_PRODUCT",
"Size is different");
358 for (x = 0; x < this->size_x; x++) {
359 res += matrix_array[x] * W[x];
371 ofstream output(filename.c_str());
372 if (output==NULL && (filename.find(
"Debug") == string::npos)) {
373 throw error_msg(
"FILE",
"Unable to output file: %s", filename.c_str());
375 output <<
"VARIABLES = \"" << ((this->name!=
"")?this->name:
"function") <<
"\" "<< endl;
376 output <<
"ZONE T=\"" << filename <<
"\", I=" << size_x << endl;
377 output.setf(ios_base::scientific, ios_base::floatfield);
378 for (x = 0; x < size_x; x++) {
379 output << matrix_array[x] << endl;
390 ofstream output(filename.c_str());
391 if (output==NULL && (filename.find(
"Debug") == string::npos)) {
392 throw error_msg(
"FILE",
"Unable to output file: %s", filename.c_str());
394 output <<
"VARIABLES = \"" << ((grid_x.
name!=
"")?grid_x.
name:
"x") <<
"\", \"" << ((this->name!=
"")?this->name:
"function") <<
"\" "<< endl;
395 output <<
"ZONE T=\"" << filename <<
"\", I=" << size_x << endl;
396 output.setf(ios_base::scientific, ios_base::floatfield);
397 for (x = 0; x < size_x; x++) {
398 output <<
"\t" << grid_x[x] <<
"\t" <<
"\t" << matrix_array[x] << endl;
412 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
414 ifstream input(filename.c_str());
415 if (input != NULL && !input.eof()) {
419 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
423 input.ignore(9999,
'\n');
424 for (x = 0; x < size_x; x++) {
425 input >> matrix_array[x];
428 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
447 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
449 ifstream input(filename.c_str());
450 if (input != NULL && !input.eof()) {
455 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
459 input.ignore(9999,
'\n');
460 for (x = 0; x < size_x; x++) {
463 if (fabs(loaded_x - grid_x[x]) > err) {
464 throw error_msg(
"MATRIX_LOAD_GRID_ERR",
"Loading %s: grid mismatch.\nLoaded: %e\nGrid: %e\n", filename.c_str(), loaded_x, grid_x[x]);
466 input >> matrix_array[x];
470 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
497 AllocateMemory(x_size, y_size);
509 this->operator = (M);
517 if (initialized) free_matrix<T>(matrix_array);
528 this->size_x = x_size;
529 this->size_y = y_size;
531 matrix_array = matrix<T>(x_size, y_size);
550 if (initialized && (size_x != M.
size_x || size_y != M.
size_y)) {
551 free_matrix<T>(matrix_array);
558 this->AllocateMemory(this->size_x, this->size_y);
564 memcpy(matrix_array[0], M.
matrix_array[0], this->size_x * this->size_y *
sizeof(
T) );
566 this->initialized =
false;
581 for (x = 0; x < this->size_x; x++)
582 for (y = 0; y < this->size_y; y++)
583 matrix_array[x][y] = val;
585 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
600 for (x = 0; x < this->size_x; x++)
601 for (y = 0; y < this->size_y; y++)
602 Tmp[x][y] = matrix_array[x][y] * Val;
615 for (x = 0; x < this->size_x; x++)
616 for (y = 0; y < this->size_y; y++)
617 Tmp[x][y] = matrix_array[x][y] / Val;
630 for (x = 0; x < this->size_x; x++)
631 for (y = 0; y < this->size_y; y++)
632 Tmp[x][y] = matrix_array[x][y] / M[x][y];
645 for (x = 0; x < this->size_x; x++)
646 for (y = 0; y < this->size_y; y++)
647 Tmp[x][y] = matrix_array[x][y] * M[x][y];
660 for (x = 0; x < this->size_x; x++)
661 for (y = 0; y < this->size_y; y++)
662 Tmp[x][y] = (matrix_array[x][y]>val)?matrix_array[x][y]:val;
687 ofstream output(filename.c_str());
688 if (output==NULL && (filename.find(
"Debug") == string::npos)) {
689 throw error_msg(
"FILE",
"Unable to output file: %s", filename.c_str());
691 output <<
"VARIABLES = \""<< ((this->name!=
"")?this->name:
"f") <<
"\" "<< endl;
692 output <<
"ZONE T=\"" << filename <<
"\", I=" << size_y <<
", J= " << size_x << endl;
693 output.setf(ios_base::scientific, ios_base::floatfield);
694 for (x = 0; x < size_x; x++) {
695 for (y = 0; y < size_y; y++) {
696 output << matrix_array[x][y] << endl;
710 ofstream output(filename.c_str());
711 output <<
"VARIABLES = \"" << ((grid_x.
name!=
"")?grid_x.
name:
"x") <<
"\", \"" << ((grid_y.
name!=
"")?grid_y.
name:
"y") <<
"\", \"" << ((this->name!=
"")?this->name:
"f") <<
"\" "<< endl;
712 output <<
"ZONE T=\"" << filename <<
"\", I=" << size_y <<
", J=" << size_x << endl;
713 output.setf(ios_base::scientific, ios_base::floatfield);
714 for (x = 0; x < size_x; x++) {
715 for (y = 0; y < size_y; y++) {
716 output <<
"\t" << grid_x[x][y] <<
"\t" << grid_y[x][y] <<
"\t" << matrix_array[x][y] << endl;
731 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
733 ifstream input(filename.c_str());
734 if (input != NULL && !input.eof()) {
738 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
742 input.ignore(9999,
'\n');
743 for (x = 0; x < size_x; x++) {
744 for (y = 0; y < size_y; y++) {
745 input >> matrix_array[x][y];
749 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
766 double loaded_x, loaded_y;
768 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
770 ifstream input(filename.c_str());
771 if (input != NULL && !input.eof()) {
775 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
779 input.ignore(9999,
'\n');
780 for (x = 0; x < size_x; x++) {
781 for (y = 0; y < size_y; y++) {
782 input >> loaded_x >> loaded_y;
784 if (fabs(loaded_x - grid_x[x][y]) > err || fabs(loaded_y - grid_y[x][y]) > err) {
785 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]);
787 input >> matrix_array[x][y];
792 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
816 AllocateMemory(x_size, y_size, z_size);
828 this->operator = (M);
836 if (initialized) free_matrix<T>(matrix_array, size_x, size_y);
845 this->size_x = x_size;
846 this->size_y = y_size;
847 this->size_z = z_size;
848 matrix_array = matrix<T>(x_size, y_size, z_size);
849 plane_array = matrix_array[0][0];
853 for (
int i = 0; i < size_x; i++)
854 for (
int j = 0; j < size_y; j++)
855 for (
int k = 0; k < size_z; k++)
856 matrix_array[i][j][k] = 0;
870 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
874 return matrix_array[i];
886 throw error_msg(
"MATRIX_ERROR",
"Using not initialized matrix");
888 if ((x < 0 || x > size_x-1) || (y < 0 || y > size_y-1) || (z < 0 || z > size_z-1)) {
889 throw error_msg(
"MATRIX_ERROR",
"Index is out of bound");
893 return plane_array[(x*size_y + y)*size_z + z];
908 if (initialized && (size_x != M.
size_x || size_y != M.
size_y || size_z != M.
size_z)) {
909 free_matrix<T>(matrix_array, this->size_x, this->size_y);
917 this->AllocateMemory(this->size_x, this->size_y, this->size_z);
924 memcpy( matrix_array[0][0], M.
matrix_array[0][0], this->size_x * this->size_y * this->size_z *
sizeof(
T ) );
926 this->initialized =
false;
963 for (
int x = 0; x < this->size_x; x++)
964 for (
int y = 0; y < this->size_y; y++)
965 for (
int z = 0; z < this->size_z; z++)
966 this->matrix_array[x][y][z] = Val;
977 for (x = 0; x < this->size_x; x++)
978 for (y = 0; y < this->size_y; y++)
979 for (z = 0; z < this->size_z; z++)
990 for (x = 0; x < this->size_x; x++)
991 for (y = 0; y < this->size_y; y++)
992 for (z = 0; z < this->size_z; z++)
1003 for (x = 0; x < this->size_x; x++)
1004 for (y = 0; y < this->size_y; y++)
1005 for (z = 0; z < this->size_z; z++)
1006 matrix_array[x][y][z] *= Val;
1016 for (x = 0; x < this->size_x; x++)
1017 for (y = 0; y < this->size_y; y++)
1018 for (z = 0; z < this->size_z; z++)
1019 matrix_array[x][y][z] /= Val;
1029 for (x = 0; x < this->size_x; x++)
1030 for (y = 0; y < this->size_y; y++)
1031 for (z = 0; z < this->size_z; z++)
1032 matrix_array[x][y][z] += Val;
1042 for (x = 0; x < this->size_x; x++)
1043 for (y = 0; y < this->size_y; y++)
1044 for (z = 0; z < this->size_z; z++)
1045 matrix_array[x][y][z] -= Val;
1055 for (x = 0; x < this->size_x; x++)
1056 for (y = 0; y < this->size_y; y++)
1057 for (z = 0; z < this->size_z; z++)
1068 for (x = 0; x < this->size_x; x++)
1069 for (y = 0; y < this->size_y; y++)
1070 for (z = 0; z < this->size_z; z++)
1082 for (x = 0; x < this->size_x; x++)
1083 for (y = 0; y < this->size_y; y++)
1084 for (z = 0; z < this->size_z; z++)
1085 Tmp[x][y][z] = matrix_array[x][y][z] + M.
matrix_array[x][y][z];
1096 for (x = 0; x < this->size_x; x++)
1097 for (y = 0; y < this->size_y; y++)
1098 for (z = 0; z < this->size_z; z++)
1099 Tmp[x][y][z] = matrix_array[x][y][z] - M.
matrix_array[x][y][z];
1111 for (x = 0; x < this->size_x; x++)
1112 for (y = 0; y < this->size_y; y++)
1113 for (z = 0; z < this->size_z; z++)
1114 Tmp[x][y][z] = matrix_array[x][y][z] * Val;
1126 for (x = 0; x < this->size_x; x++)
1127 for (y = 0; y < this->size_y; y++)
1128 for (z = 0; z < this->size_z; z++)
1129 Tmp[x][y][z] = matrix_array[x][y][z] / Val;
1140 for (x = 0; x < this->size_x; x++)
1141 for (y = 0; y < this->size_y; y++)
1142 for (z = 0; z < this->size_z; z++)
1143 Tmp[x][y][z] = matrix_array[x][y][z] * M.
matrix_array[x][y][z];
1154 for (x = 0; x < this->size_x; x++)
1155 for (y = 0; y < this->size_y; y++)
1156 for (z = 0; z < this->size_z; z++)
1157 Tmp[x][y][z] = matrix_array[x][y][z] / M.
matrix_array[x][y][z];
1169 ofstream output(filename.c_str());
1170 output <<
"VARIABLES = \""<< ((this->name!=
"")?this->name:
"f") <<
"\" "<< endl;
1171 output <<
"ZONE T=\"" << filename <<
"\", I=" << size_z <<
", J=" << size_y <<
", K=" << size_x << endl;
1172 output.setf(ios_base::scientific, ios_base::floatfield);
1173 for (x = 0; x < size_x; x++) {
1174 for (y = 0; y < size_y; y++) {
1175 for (z = 0; z < size_z; z++) {
1176 output << matrix_array[x][y][z] << endl;
1190 ofstream output(filename.c_str());
1191 if (output==NULL && (filename.find(
"Debug") == string::npos)) {
1192 throw error_msg(
"FILE",
"Unable to output file: %s", filename.c_str());
1194 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;
1195 output <<
"ZONE T=\"" << filename <<
"\", I=" << size_z <<
", J=" << size_y <<
", K=" << size_x << endl;
1196 output.setf(ios_base::scientific, ios_base::floatfield);
1197 for (x = 0; x < size_x; x++) {
1198 for (y = 0; y < size_y; y++) {
1199 for (z = 0; z < size_z; z++) {
1200 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;
1215 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
1217 ifstream input(filename.c_str());
1218 if (input != NULL && !input.eof()) {
1222 while (strcasecmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
1227 input.ignore(9999,
'\n');
1228 for (x = 0; x < size_x; x++) {
1229 for (y = 0; y < size_y; y++) {
1230 for (z = 0; z < size_z; z++) {
1231 input >> matrix_array[x][y][z];
1236 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
1253 throw error_msg(
"MATRIX_ERROR",
"Using unitialized matrix");
1255 ifstream input(filename.c_str());
1256 if (input != NULL && !input.eof()) {
1261 while (strcmp(inBuf.c_str(),
"ZONE") != 0 && !input.eof() ) {
1265 input.ignore(9999,
'\n');
1266 for (x = 0; x < size_x; x++) {
1267 for (y = 0; y < size_y; y++) {
1268 for (z = 0; z < size_z; z++) {
1269 input >> grid_x[x][y][z] >> grid_y[x][y][z] >> grid_z[x][y][z] >> matrix_array[x][y][z];
1271 input.ignore(9999,
'\n');
1284 throw error_msg(
"MATRIX_LOAD_ERROR",
"Error reading file %s.\n", filename.c_str());
1300 return (x*size_y + y)*size_z + z;
1311 for (x = 0; x < size_x; x++) {
1312 for (y = 0; y < size_y; y++) {
1313 for (z = 0; z < size_z; z++) {
1314 tmp = (tmp>matrix_array[x][y][z])?tmp:matrix_array[x][y][z];
1328 for (x = 0; x < size_x; x++) {
1329 for (y = 0; y < size_y; y++) {
1330 for (z = 0; z < size_z; z++) {
1331 tmp = (tmp>fabs((
double)matrix_array[x][y][z]))?tmp:fabs((
double)matrix_array[x][y][z]);
1343 Matrix3D<T> tmp(this->size_x, this->size_y, this->size_z);
1345 for (x = 0; x < size_x; x++) {
1346 for (y = 0; y < size_y; y++) {
1347 for (z = 0; z < size_z; z++) {
1348 tmp[x][y][z] = (matrix_array[x][y][z]>0)?matrix_array[x][y][z]:-matrix_array[x][y][z];
1363 tmp.
name = this->name +
"_slice";
1364 for (y = 0; y < size_y; y++) {
1365 for (z = 0; z < size_z; z++) {
1366 tmp[y][z] = matrix_array[p_x][y][z];
1379 tmp.
name = this->name +
"_slice";
1380 for (x = 0; x < size_x; x++) {
1381 for (z = 0; z < size_z; z++) {
1382 tmp[x][z] = matrix_array[x][p_y][z];
1395 tmp.
name = this->name +
"_slice";
1396 for (x = 0; x < size_x; x++) {
1397 for (y = 0; y < size_y; y++) {
1398 tmp[x][y] = matrix_array[x][y][p_z];
1421 if (this->size_x != new_grid.
size_x) {
1422 throw error_msg(
"INTERPOLATION_ERROR",
"New grid and function size mismatch!\n");
1433 for (i = 0; i < new_grid.
size_x; i++)
1462 if (this->size_x != new_grid.
size_x) {
1463 throw error_msg(
"SPLINE_INTERPOLATION",
"New grid and function size mismatch!");
1466 if (this->size_x == 1) {
1472 int i=0, i_min=0, i_max=old_grid.
size_x-1, interp_size=0;
1480 i_min = 1; i_max = old_grid.
size_x-2;
1481 interp_size = i_max-i_min+1;
1500 for (i = 0; i < old_grid.
size_x; i++) {
1502 if (fabs(y2[i]) >= max_second_der) {
1506 y2[i] = max_second_der + (y2[i] - max_second_der)*(1.0 - lin_spline_coef);
1508 y2[i] = - max_second_der + (y2[i] + max_second_der)*(1.0 - lin_spline_coef);
1515 for (i = 0; i < new_grid.
size_x; i++) {
1517 this->matrix_array[i] = lb;
1518 }
else if (i == new_grid.
size_x-1) {
1519 this->matrix_array[i] = ub;
1521 this->matrix_array[i] = lb;
1524 this->matrix_array[i] = ub;
1628 if (this->size_x != new_grid.
size_x) {
1629 throw error_msg(
"SPLINE_INTERPOLATION",
"New grid and function size mismatch!");
1632 if (this->size_x == 1) {
1638 int i=0, i_min=0, i_max=old_grid.
size_x-1, interp_size=0;
1646 i_min = 0; i_max = old_grid.
size_x-1;
1647 interp_size = i_max-i_min+1;
1664 for (i = 0; i < new_grid.
size_x; i++) {
1679 if (this->size_x != new_grid.
size_x) {
1680 throw error_msg(
"POLINOMIAL_INTERPOLATION",
"New grid and function size mismatch!");
1684 for (i = 0; i < new_grid.
size_x; i++)
1697 if (this->size_x != new_grid.
size_x) {
1698 throw error_msg(
"POLINT_INTERPOLATION",
"New grid and function size mismatch!");
1703 for (i = 0; i < new_grid.
size_x; i++)
1713 if (this->size_x != new_grid.
size_x) {
1714 throw error_msg(
"RATIONAL_INTERPOLATION",
"New grid and function size mismatch!");
1719 for (i = 0; i < new_grid.
size_x; i++)
1862 matrix_array[0] = lb;
1863 matrix_array[size_x-1] = ub;
1880 double Linear2D(
double x,
double y,
T &old_grid_x,
T &old_grid_y,
T &f) {
1884 while (x >= old_grid_x[++i][0] && i < old_grid_y.size_y-1);
1885 while (y >= old_grid_y[0][++j] && j < old_grid_x.size_x-1);
1887 double t = (x - old_grid_x[i-1][0])/(old_grid_x[i][0] - old_grid_x[i-1][0]);
1888 double u = (y - old_grid_y[0][j-1])/(old_grid_y[0][j] - old_grid_y[0][j-1]);
1890 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];
1905 throw error_msg(
"LINEAR_INTERPOLATION",
"Old grid and function size mismatch!");
1908 if (this->size_x != new_grid_x.
size_x || this->size_y != new_grid_y.
size_y) {
1909 throw error_msg(
"LINEAR_INTERPOLATION",
"New grid and function size mismatch!");
1913 for (i = 0; i < new_grid_x.
size_x; i++) {
1914 for (j = 0; j < new_grid_y.
size_y; j++) {
1915 this->matrix_array[i][j] =
Linear2D(new_grid_x[i][j], new_grid_y[i][j], old_grid_x, old_grid_y, old_function);
1932 this->initialized =
false;
1933 Initialize(x_size, y_size, z_size, n_of_diags);
1940 this->initialized =
false;
1941 this->x_size = x_size;
1942 this->y_size = y_size;
1943 this->total_size = x_size;
1944 if (y_size > 0) this->total_size = this->total_size * y_size;
1945 if (z_size > 0) this->total_size = this->total_size * z_size;
1952 for (x = -n_of_diags; x <= n_of_diags; x++) {
1953 for (y = -n_of_diags; y <= n_of_diags; y++) {
1954 for (z = -n_of_diags; z <= n_of_diags; z++) {
1956 id = this->index1d(x, y, z);
1964 }
else if (y_size > 1) {
1965 for (x = -n_of_diags; x <= n_of_diags; x++) {
1966 for (y = -n_of_diags; y <= n_of_diags; y++) {
1968 id = this->index1d(x, y);
1975 }
else if (x_size > 1) {
1976 for (x = -n_of_diags; x <= n_of_diags; x++) {
1978 id = this->index1d(x);
1986 if (y_size == 0) this->x_size = 0;
1988 if (z_size == 0) this->y_size = 0;
1990 this->initialized =
true;
1997 return (z*y_size + y)*x_size + x;
2006 ofstream output(filename.c_str());
2007 if (output==NULL && (filename.find(
"Debug") == string::npos)) {
2008 throw error_msg(
"FILE",
"Unable to output file: %s", filename.c_str());
2010 output <<
"VARIABLES = ";
2011 for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2012 output <<
"\"" << it->first <<
"\", ";
2015 output <<
"ZONE T=\"" <<
"\", I=" << total_size << endl;
2016 output.setf(ios_base::scientific, ios_base::floatfield);
2017 for (in = 0; in < this->total_size; in++) {
2018 for (DiagMatrix::iterator it = (*this).begin(); it != (*this).end(); it++) {
2019 output <<
"\t" << it->second[in];