00001
00019
00020
00021
00022
00023 #include "MatrixSolver.h"
00024 #include <math.h>
00025 #include <malloc.h>
00026 #include <string>
00027 #include <ctime>
00028
00029 using namespace std;
00030
00031
00032 #include "../VariousFunctions/variousFunctions.h"
00033 #include "../Exceptions/error.h"
00034
00035
00036 #include "../Lapack/Include/cpplapack.h"
00037
00038
00039 extern "C" {
00040 #include "../GMRES/itlin.h"
00041 }
00042
00043 void matvec(int n, double *x, double *b);
00044 void preconr(int n, double *x, double *b);
00045 extern "C" double gmres_norm2(int n, double *v);
00046
00047 extern "C" void gmres(int n, double *y, MATVEC *matvec,
00048 PRECON *preconr, PRECON *preconl, double *b,
00049 struct ITLIN_OPT *opt, struct ITLIN_INFO *info);
00050
00051 #include <iostream>
00052
00056 bool MakeMatrix(double *A,
00057 double *B1,
00058 double *C,
00059 double *R,
00060 double *x,
00061 double tau,
00062 double n,
00063 double f_lower,
00064 double f_upper,
00065 int nx,
00066 double dt,
00067 double *Dxx,
00068 double taulc,
00069 double alc,
00070 int g_flag,
00071 string lower_border_condition_type,
00072 string upper_border_condition_type,
00073 string approximationMethod)
00074 {
00075
00076
00077 double Coef1, Dm_05, Dp_05, tau_cur;
00078
00079
00080
00081 double *B = new double[nx];
00082
00083
00084 double *h = new double[nx];
00085 double *h_bar = new double[nx];
00086
00087
00088 int i, j;
00089 double coef_05m, coef_05p;
00090
00091
00092 for (j = 1; j <= nx-1; j++) {
00093 h[j] = x[j] - x[j-1];
00094 }
00095 for (j = 0; j <= nx-2; j++) {
00096 h_bar[j] = (h[j] + h[j+1])/2;
00097 }
00098
00099 for (j = 1; j <= nx-2; j++) {
00100
00101 if (approximationMethod == "AM_Split_LR") {
00102
00103 Coef1 = pow(x[j], n)/h_bar[j];
00104 if (g_flag == 1) {
00105 Coef1 = Coef1/VF::G(x[j]);
00106 coef_05m = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j-1]*VF::G(x[j-1])*pow(x[j-1], -n));
00107 coef_05p = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j+1]*VF::G(x[j+1])*pow(x[j+1], -n));
00108 } else {
00109 coef_05m = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j-1]*pow(x[j-1], -n));
00110 coef_05p = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j+1]*pow(x[j+1], -n));
00111 }
00112 A[j] = Coef1 * coef_05m / h[j];
00113 C[j] = Coef1 * coef_05p / h[j+1];
00114
00115 } else if (approximationMethod == "AM_Split_LR_oldway") {
00116
00117 Coef1 = pow(x[j], n);
00118 if (g_flag == 1) {
00119 Coef1 = Coef1/VF::G(x[j]);
00120 coef_05m = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j-1]*VF::G(x[j-1])*pow(x[j-1], -n));
00121 coef_05p = 0.5*(Dxx[j]*VF::G(x[j])*pow(x[j], -n) + Dxx[j+1]*VF::G(x[j+1])*pow(x[j+1], -n));
00122 } else {
00123 coef_05m = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j-1]*pow(x[j-1], -n));
00124 coef_05p = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j+1]*pow(x[j+1], -n));
00125 }
00126 A[j] = Coef1 * coef_05m / h[j] / h[j+1];
00127 C[j] = Coef1 * coef_05p / h[j+1] / h[j];
00128
00129 } else if (approximationMethod == "AM_Split_C") {
00130
00131 Coef1 = pow(x[j], n)/h_bar[j];
00132 if (g_flag == 1) {
00133 Coef1 = Coef1/VF::G(x[j]);
00134 }
00135 Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
00136 Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
00137 if (g_flag == 0) {
00138 A[j] = Coef1 * (Dm_05 * pow(x[j] - h[j]/2, -n) / h[j]);
00139 C[j] = Coef1 * (Dp_05 * pow(x[j] + h[j+1]/2, -n) / h[j+1]);
00140 } else {
00141 A[j] = Coef1 * (Dm_05 * VF::G(x[j] - h[j]/2) * pow(x[j] - h[j]/2, -n) / h[j]);
00142 C[j] = Coef1 * (Dp_05 * VF::G(x[j] + h[j+1]/2) * pow(x[j] + h[j+1]/2, -n) / h[j+1]);
00143 }
00144
00145 } else if (approximationMethod == "AM_Split_C_oldway") {
00146
00147 Coef1 = pow(x[j], n)/h_bar[j];
00148 if (g_flag == 1) {
00149 Coef1 = Coef1/VF::G(x[j]);
00150 }
00151 Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
00152 Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
00153 if (g_flag == 0) {
00154 A[j] = Coef1 * (Dm_05 * pow(x[j] - h[j]/2, -n) / h[j]);
00155 C[j] = Coef1 * (Dp_05 * pow(x[j] + h[j]/2, -n) / h[j+1]);
00156 } else {
00157 A[j] = Coef1 * (Dm_05 * VF::G(x[j] - h[j]/2) * pow(x[j] - h[j]/2, -n) / h[j]);
00158 C[j] = Coef1 * (Dp_05 * VF::G(x[j] + h[j]/2) * pow(x[j] + h[j]/2, -n) / h[j+1]);
00159 }
00160
00161 } else {
00162 throw error_msg("APPRERR", "Unknown approximation type: %d", approximationMethod.c_str());
00163 }
00164
00165 B[j] = -A[j] - C[j];
00166
00167
00168
00169 if ((x[j] < alc || x[j] > VC::pi - alc) && (g_flag == 1)) {
00170 tau_cur = taulc;
00171 } else {
00172 tau_cur = tau;
00173 }
00174
00175 B1[j] = B[j] - (1.0/dt + 1/tau_cur);
00176 }
00177
00178 for (i = 1; i <= nx-2; i++) {
00179 R[i] = - 1.0 / dt;
00180 }
00181
00182
00183
00184 i = 0;
00185 R[i] = f_lower;
00186 if (lower_border_condition_type == "BCT_CONSTANT_VALUE") {
00187 A[i] = 0; B1[i] = 1; C[i] = 0;
00188 } else if (lower_border_condition_type == "BCT_CONSTANT_DERIVATIVE") {
00189 A[i] = 0; B1[i] = -1; C[i] = 1;
00190 } else {
00191 throw error_msg("1D_DIFF_BOUNDARY", "1d_universal_solver: unknown lower boundary type: %d", lower_border_condition_type.c_str());
00192
00193
00194 }
00195
00196
00197 i = nx-1;
00198 R[i] = f_upper;
00199 if (upper_border_condition_type == "BCT_CONSTANT_VALUE") {
00200 A[i] = 0; B1[i] = 1; C[i] = 0;
00201 } else if (upper_border_condition_type == "BCT_CONSTANT_DERIVATIVE") {
00202 A[i] = -1; B1[i] = 1; C[i] = 0;
00203 } else {
00204 throw error_msg("1D_DIFF_BOUNDARY", "1d_universal_solver: unknown lower boundary type: %d", upper_border_condition_type.c_str());
00205
00206
00207 }
00208
00209
00210
00211
00212
00213
00214 delete B;
00215
00216
00217 delete h;
00218 delete h_bar;
00219
00220 return true;
00221 }
00222
00226 bool SolveMatrix(double *f,
00227 double *A,
00228 double *B1,
00229 double *C,
00230 double *R,
00231 double f_lower,
00232 double f_upper,
00233 int nx,
00234 double dt)
00235 {
00236
00237 double *R1 = new double[nx];
00238
00239
00240 int i;
00241
00242 for (i = 1; i <= nx-2; i++) {
00243 R1[i] = R[i] * f[i];
00244 }
00245
00246
00247
00248 i = 0;
00249 R1[i] = R[i];
00250
00251
00252 i = nx-1;
00253 R1[i] = R[i];
00254
00255
00256 tridag(A, B1, C, R1, f, nx);
00257
00258
00259 delete R1;
00260
00261 return true;
00262 }
00263
00264
00267 void AddBoundary(DiagMatrix &Matrix_A, string type, int in, int id1, double dh) {
00268 int id0 = 0;
00269 if (type == "BCT_CONSTANT_VALUE") {
00270 Matrix_A[id0][in] = 1;
00271 Matrix_A[id1][in] = 0;
00272 } else if (type == "BCT_CONSTANT_DERIVATIVE") {
00273 Matrix_A[id0][in] = -1/dh;
00274 Matrix_A[id1][in] = 1/dh;
00275 } else {
00276 throw error_msg("2D_DIFF_BOUNDARY", "1d_universal_solver: unknown boundary type: %d", type.c_str());
00277 }
00278 }
00279
00287 bool MakeModelMatrix_3D(
00288 CalculationMatrix &matr_A, CalculationMatrix &matr_B, CalculationMatrix &matr_C,
00289 Matrix3D<double> &L, Matrix3D<double> &pc, Matrix3D<double> &alpha,
00290 int L_size, int pc_size, int alpha_size,
00291 Matrix2D<double> &L_lowerBoundaryCondition,
00292 Matrix2D<double> &L_upperBoundaryCondition,
00293 Matrix2D<double> &pc_lowerBoundaryCondition,
00294 Matrix2D<double> &pc_upperBoundaryCondition,
00295 Matrix2D<double> &alpha_lowerBoundaryCondition,
00296 Matrix2D<double> &alpha_upperBoundaryCondition,
00297 string L_lowerBoundaryCondition_calculationType,
00298 string L_upperBoundaryCondition_calculationType,
00299 string pc_lowerBoundaryCondition_calculationType,
00300 string pc_upperBoundaryCondition_calculationType,
00301 string alpha_lowerBoundaryCondition_calculationType,
00302 string alpha_upperBoundaryCondition_calculationType,
00303 Matrix3D<double> &DLL,
00304 Matrix3D<double> &Dpcpc, Matrix3D<double> &DpcpcLpp,
00305 Matrix3D<double> &Daa, Matrix3D<double> &DaaLpp,
00306 Matrix3D<double> &Dpca, Matrix3D<double> &DpcaLpp,
00307 Matrix3D<double> &Jacobian, double dt, double Lpp,
00308 double tau, double tauLpp) {
00309
00310
00311
00312
00313
00314 int il, im, ia, id, in;
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 DiagMatrix::iterator it;
00336 for (it = matr_A.begin(); it != matr_A.end(); it++) it->second = 0;
00337 for (it = matr_B.begin(); it != matr_B.end(); it++) it->second = 0;
00338 for (it = matr_C.begin(); it != matr_C.end(); it++) it->second = 0;
00339
00340
00341
00342
00343
00344 double dh;
00345 for (il = 0; il < L_size; il++) {
00346 for (im = 0; im < pc_size; im++) {
00347 for (ia = 0; ia < alpha_size; ia++) {
00348
00349 in = matr_A.index1d(ia, im, il);
00350
00351
00352
00353
00354 if (il == 0 && L_size >= 3) {
00355
00356
00357
00358 matr_C[0][in] = L_lowerBoundaryCondition[im][ia];
00359 id = matr_A.index1d(ia, im, il+1) - in;
00360 dh = L[il+1][im][ia] - L[il][im][ia];
00361 AddBoundary(matr_A, L_lowerBoundaryCondition_calculationType, in, id, dh);
00362
00363 } else if (il == L_size-1 && L_size >= 3) {
00364
00365 matr_C[0][in] = L_upperBoundaryCondition[im][ia];
00366 id = matr_A.index1d(ia, im, il-1) - in;
00367 dh = L[il][im][ia] - L[il-1][im][ia];
00368 AddBoundary(matr_A, L_upperBoundaryCondition_calculationType, in, id, dh);
00369
00370 } else if (im == 0 && pc_size >= 3) {
00371
00372 matr_C[0][in] = pc_lowerBoundaryCondition[il][ia];
00373 id = matr_A.index1d(ia, im+1, il) - in;
00374 dh = pc[il][im+1][ia] - pc[il][im][ia];
00375 AddBoundary(matr_A, pc_lowerBoundaryCondition_calculationType, in, id, dh);
00376
00377 } else if (im == pc_size-1 && pc_size >= 3) {
00378
00379 matr_C[0][in] = pc_upperBoundaryCondition[il][ia];
00380 id = matr_A.index1d(ia, im-1, il) - in;
00381 dh = pc[il][im][ia] - pc[il][im-1][ia];
00382 AddBoundary(matr_A, pc_upperBoundaryCondition_calculationType, in, id, dh);
00383 } else if (ia == 0 && alpha_size >= 3) {
00384
00385 matr_C[0][in] = alpha_lowerBoundaryCondition[il][im];
00386 id = matr_A.index1d(ia+1, im, il) - in;
00387 dh = alpha[il][im][ia+1] - alpha[il][im][ia];
00388 AddBoundary(matr_A, alpha_lowerBoundaryCondition_calculationType, in, id, dh);
00389
00390 } else if (ia == alpha_size-1 && alpha_size >= 3) {
00391
00392 matr_C[0][in] = alpha_upperBoundaryCondition[il][im];
00393 id = matr_A.index1d(ia-1, im, il) - in;
00394 dh = alpha[il][im][ia] - alpha[il][im][ia-1];
00395 AddBoundary(matr_A, alpha_upperBoundaryCondition_calculationType, in, id, dh);
00396
00397 } else {
00398
00399
00400
00401
00402 matr_A[0][in] += 1.0/dt;
00403
00404 matr_B[0][in] += 1.0/dt;
00405
00406
00407 if (L_size >=3) {
00408
00409 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00410 "DT_L_left", "DT_L_right",
00411 L, pc, alpha, DLL, Jacobian, -0.5);
00412 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00413 "DT_L_right", "DT_L_left",
00414 L, pc, alpha, DLL, Jacobian, -0.5);
00415 }
00416
00417 if (pc_size >=3) {
00418
00419 if (L[il][im][ia] >= Lpp) {
00420 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00421 "DT_PC_left", "DT_PC_right",
00422 L, pc, alpha, Dpcpc, Jacobian, -0.5);
00423 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00424 "DT_PC_right", "DT_PC_left",
00425 L, pc, alpha, Dpcpc, Jacobian, -0.5);
00426 } else {
00427 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00428 "DT_PC_left", "DT_PC_right",
00429 L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
00430 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00431 "DT_PC_right", "DT_PC_left",
00432 L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
00433 }
00434 }
00435
00436 if (alpha_size >=3) {
00437
00438
00439
00440
00441 if (alpha[il][im][ia] <= VF::alc(L[il][im][0]) || alpha[il][im][ia] >= VC::pi - VF::alc(L[il][im][0])) {
00442 double taulc = VF::bounce_time(VF::Kfunc(pc[il][im][ia]), L[il][im][ia])/60./60./24.;
00443 matr_A[0][in] += 1.0/taulc;
00444 }
00445
00446 if (L[il][im][ia] >= Lpp) {
00447 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00448 "DT_ALPHA_left", "DT_ALPHA_right",
00449 L, pc, alpha, Daa, Jacobian, -0.5);
00450 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00451 "DT_ALPHA_right", "DT_ALPHA_left",
00452 L, pc, alpha, Daa, Jacobian, -0.5);
00453 } else {
00454 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00455 "DT_ALPHA_left", "DT_ALPHA_right",
00456 L, pc, alpha, DaaLpp, Jacobian, -0.5);
00457 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00458 "DT_ALPHA_right", "DT_ALPHA_left",
00459 L, pc, alpha, DaaLpp, Jacobian, -0.5);
00460 }
00461 }
00462
00463 if (alpha_size >= 3 && pc_size >= 3) {
00464
00465
00466 if (L[il][im][ia] >= Lpp) {
00467
00468
00469 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00470 "DT_ALPHA_left", "DT_PC_right",
00471 L, pc, alpha, Dpca, Jacobian, -0.25);
00472 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00473 "DT_ALPHA_right", "DT_PC_left",
00474 L, pc, alpha, Dpca, Jacobian, -0.25);
00475
00476
00477 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00478 "DT_PC_left", "DT_ALPHA_right",
00479 L, pc, alpha, Dpca, Jacobian, -0.25);
00480 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00481 "DT_PC_right", "DT_ALPHA_left",
00482 L, pc, alpha, Dpca, Jacobian, -0.25);
00483
00484
00485 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00486 "DT_ALPHA_left", "DT_PC_left",
00487 L, pc, alpha, Dpca, Jacobian, -0.25);
00488 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00489 "DT_ALPHA_right", "DT_PC_right",
00490 L, pc, alpha, Dpca, Jacobian, -0.25);
00491
00492
00493 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00494 "DT_PC_left", "DT_ALPHA_left",
00495 L, pc, alpha, Dpca, Jacobian, -0.25);
00496 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00497 "DT_PC_right", "DT_ALPHA_right",
00498 L, pc, alpha, Dpca, Jacobian, -0.25);
00499
00500 } else {
00501
00502
00503 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00504 "DT_ALPHA_left", "DT_PC_right",
00505 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00506 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00507 "DT_ALPHA_right", "DT_PC_left",
00508 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00509
00510
00511 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00512 "DT_PC_left", "DT_ALPHA_right",
00513 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00514 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00515 "DT_PC_right", "DT_ALPHA_left",
00516 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00517
00518
00519
00520 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00521 "DT_ALPHA_left", "DT_PC_left",
00522 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00523 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00524 "DT_ALPHA_right", "DT_PC_right",
00525 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00526
00527
00528 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00529 "DT_PC_left", "DT_ALPHA_left",
00530 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00531 SecondDerivativeApproximation_3D(matr_A, il, im, ia,
00532 "DT_PC_right", "DT_ALPHA_right",
00533 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
00534
00535 }
00536 }
00537 }
00538 }
00539 }
00540 }
00541
00542
00543
00544 matr_A.change_ind = clock();
00545 matr_B.change_ind = clock();
00546 matr_C.change_ind = clock();
00547
00548 return true;
00549 }
00550
00551
00552
00553
00554 static CalculationMatrix A_copy;
00555
00559 void matvec(int m_size, double *x, double *b)
00560 {
00561 DiagMatrix::iterator it;
00562 int in, ik;
00563 for (in = 0; in < m_size; in++) {
00564
00565 b[in] = 0;
00566
00567 for (it = A_copy.begin(); (it != A_copy.end()); it++) {
00568
00569
00570
00571
00572 ik = in + it->first;
00573 if (ik >= 0 && ik < m_size) {
00574
00575
00576
00577
00578
00579 b[in] += it->second[in] * x[ik];
00580
00581 }
00582 }
00583
00584 }
00585 }
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00601 void jacobi(int m_size, double *z, double *w)
00602 {
00603 int in;
00604 for (in = 0; in < m_size; in++) {
00605 if (A_copy[0][in] != 0) {
00606 w[in] = z[in] / A_copy[0][in];
00607 } else {
00608 w[in] = z[in];
00609 }
00610 }
00611 }
00612
00621 void SOR(int m_size, double *z, double *w)
00622 {
00623 int in;
00624 DiagMatrix::iterator it;
00625 double diag;
00626 double rel = 1;
00627 for (in = 0; in < m_size; in++) {
00628 diag = A_copy[0][in] / rel;
00629 w[in] = z[in] / diag;
00630
00631 for (it = A_copy.begin(); (it->first < 0 && it != A_copy.end()); it++) {
00632 if (in + it->first >= 0 && it->first < 0) {
00633
00634 w[in] -= w[in + it->first] * it->second[in] / diag;
00635 }
00636 }
00637 }
00638 }
00639
00643 void for_norm_A(int n, double *x, double *b)
00644 {
00645 int j;
00646 for (j=0;j<n;j++) b[j] = x[j];
00647 }
00648
00649
00656 void gmres_wrapout(CalculationMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X, int maxiter, int i_max, double max_err) {
00657 int m_size = A[0].size_x;
00658
00659
00660 int in;
00661
00662
00663 DiagMatrix::iterator it;
00664
00665
00666 std::clock_t start_time, current_time;
00667
00668 double err_max = 0.;
00669 double error;
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684 A_copy = A;
00685 Matrix1D<double> B_copy;
00686 B_copy = B;
00687
00688
00689
00690
00691
00692 double A0;
00693 for (in = 0; in < m_size; in++) {
00694 A0 = A_copy[0][in];
00695 B_copy[in] = B_copy[in] / A0;
00696 for (it = A_copy.begin(); it != A_copy.end(); it++) {
00697 it->second[in] = it->second[in]/A0;
00698 }
00699 }
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723 struct ITLIN_OPT *opt = new struct ITLIN_OPT;
00724 struct ITLIN_INFO *info = new struct ITLIN_INFO;
00725
00726 TERM_CHECK termcheck = CheckOnRestart;
00727
00728
00729 opt->tol = max_err;
00730 opt->i_max = i_max;
00731 opt->termcheck = termcheck;
00732 opt->maxiter = maxiter;
00733 opt->errorlevel = Verbose;
00734 opt->monitorlevel = Verbose;
00735 opt->datalevel = None;
00736 opt->errorfile = stdout;
00737 opt->monitorfile = stdout;
00738 opt->datafile = NULL;
00739
00740
00741 opt->iterfile = NULL;
00742 opt->resfile = NULL;
00743 opt->miscfile = NULL;
00744
00745 start_time = std::clock();
00746 Output::echo(5, "Inverting... \n");
00747
00748 gmres(m_size, &X[0], &matvec, NULL, NULL, &B_copy[0], opt, info);
00749
00750
00751
00752
00753 current_time = std::clock();
00754 Output::echo(5, "Done (%.2f sec)\n", (( current_time - start_time ) / (double)CLOCKS_PER_SEC ));
00755
00756
00757
00758
00759
00760 fprintf(stdout,"\n Return code: %6i\n",info->rcode);
00761 fprintf(stdout," Iterations: %6i\n",info->iter);
00762 fprintf(stdout," Matvec calls: %6i\n",info->nomatvec);
00763 fprintf(stdout," Right Precon calls: %6i\n",info->noprecr);
00764 fprintf(stdout," Left Precon calls: %6i\n",info->noprecl);
00765 fprintf(stdout," Estimated residual reduction: %e\n",info->precision);
00766
00767 error = 0; err_max = 0;
00768 for (in = 0; in < m_size; in++) {
00769
00770 error = B[in];
00771 for (it = A.begin(); it != A.end(); it++) {
00772 if (in + it->first >= 0 && in + it->first < m_size) {
00773 error -= it->second[in] * X[in + it->first];
00774 }
00775 }
00776 err_max = max(err_max, fabs(error));
00777 }
00778 Output::echo(0, "Error_norm = %e \n", err_max);
00779 err_max = 0;
00780
00781
00782 }
00783
00784
00785
00791 void Lapack(DiagMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X) {
00792
00793
00794 std::clock_t start_time, current_time;
00795
00796
00797 DiagMatrix::iterator it;
00798
00799
00800 int modelMatrixLineNumber;
00801
00802 long i;
00803 long m_size = A[0].size_x;
00804 it = A.begin();
00805 long kl = -it->first;
00806 it = A.end();
00807 it--;
00808 long ku = it->first;
00809
00810
00811 Output::echo(5, "Creating matrixes...\n");
00812 CPPL::dgbmatrix lap_AB(m_size, m_size, kl, ku);
00813 CPPL::dgbmatrix lap_AB_res(m_size, m_size, kl, ku);
00814
00815
00816 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
00817 for (int modelMatrixDiagNumber = -kl; modelMatrixDiagNumber <= ku; modelMatrixDiagNumber++) {
00818 if (modelMatrixDiagNumber + modelMatrixLineNumber >= 0 && modelMatrixDiagNumber + modelMatrixLineNumber < m_size) {
00819 lap_AB(modelMatrixLineNumber, modelMatrixDiagNumber + modelMatrixLineNumber) = 0;
00820 }
00821 }
00822 }
00823
00824
00825
00826 Output::echo(5, "Copying values...\n");
00827 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
00828 for (it = A.begin(); it != A.end(); it++) {
00829
00830
00831
00832
00833 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
00834
00835 lap_AB(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
00836 }
00837 }
00838 }
00839
00840
00841 Output::echo(5, "matr --> res\n");
00842 lap_AB_res = lap_AB;
00843
00844
00845 CPPL::dcovector lap_B(m_size);
00846 CPPL::dcovector lap_B_res(m_size);
00847 for(i = 0; i < lap_B.l; i++){
00848 lap_B(i) = B[i];
00849 }
00850 lap_B_res = lap_B;
00851
00852 start_time = std::clock();
00853 Output::echo(5, "inverting... ");
00854
00855
00856 lap_AB.dgbsv(lap_B);
00857
00858 current_time = std::clock();
00859 Output::echo(5, "done (%.2f sec)\n", (( current_time - start_time ) / (double)CLOCKS_PER_SEC ));
00860
00861 Output::echo(5, "Checking... ");
00862 lap_B_res = lap_AB_res*lap_B - lap_B_res;
00863
00864
00865 double max = 0;
00866 for (i = 0; i < lap_B_res.l; i++) {
00867 max = (max>lap_B_res(i))?max:lap_B_res(i);
00868 }
00869 Output::echo(5, " Max error: %e.\n", max);
00870
00871
00872 Output::echo(5, "Coping back... ");
00873 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
00874 X[modelMatrixLineNumber] = lap_B(modelMatrixLineNumber);
00875 }
00876 Output::echo(5, "done. \n");
00877
00878 }
00879
00880
00884 void over_relaxation_diag(DiagMatrix &A, Matrix1D<double> &B, Matrix1D<double> &X, int max_steps, double EE) {
00885 int i, j;
00886 double omega, sum;
00887 double err_max = 1e99;
00888 int n = A[0].size_x;
00889
00890 omega = 1.75;
00891
00892
00893 DiagMatrix::iterator it;
00894 int nsteps = 0;
00895
00896 while(err_max > EE && nsteps < max_steps) {
00897 nsteps++;
00898 for (i = 0; i < n; i++) {
00899 sum = 0;
00900 for (it = A.begin(); it != A.end(); it++) {
00901 j = i + it->first;
00902 if (j >= 0 && j < n && j != i) {
00903 sum = sum + it->second[i]/A[0][i] * X[j];
00904 }
00905 }
00906 X[i] = (1 - omega) * X[i] - omega * ( sum - B[i]/A[0][i]);
00907 }
00908
00909 err_max = 0.;
00910 DiagMatrix::iterator it;
00911 double error;
00912 for (j = 0; j < n; j++) {
00913
00914 error = B[j];
00915 for (it = A.begin(); it != A.end(); it++) {
00916 if (j + it->first >= 0 && j + it->first < n) {
00917 error -= it->second[j] * X[j + it->first];
00918 }
00919 }
00920 err_max = max(err_max, fabs(error));
00921 }
00922 }
00923 Output::echo(0, "[%d]Error_norm = %e \n", nsteps, err_max);
00924 }
00925
00926
00933 void GetDerivativeVector(string derivativeType, int &dL, int &dPc, int &dAlpha) {
00934 if (derivativeType == "DT_ALPHA_left") {
00935 dAlpha = -1;
00936 } else if (derivativeType == "DT_ALPHA_right") {
00937 dAlpha = 1;
00938 } else if (derivativeType == "DT_PC_left") {
00939 dPc = -1;
00940 } else if (derivativeType == "DT_PC_right") {
00941 dPc = 1;
00942 } else if (derivativeType == "DT_L_left") {
00943 dL = -1;
00944 } else if (derivativeType == "DT_L_right") {
00945 dL = 1;
00946 } else {
00947 throw error_msg("DERIVATIVE_TYPE", "Unknown derivative approximation: %s\n", derivativeType.c_str());
00948 }
00949 }
00950
00951
00960 void SecondDerivativeApproximation_3D(CalculationMatrix &matr_A,
00961 int il, int im, int ia,
00962 string FirstDerivative, string SecondDerivative,
00963 Matrix3D<double> &L, Matrix3D<double> &pc, Matrix3D<double> &alpha,
00964 Matrix3D<double> &Dxx, Matrix3D<double> &Jacobian,
00965 double multiplicator) {
00966
00967 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
00968 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
00969
00970
00971
00972
00973
00974
00975 GetDerivativeVector(FirstDerivative, dL1, dPc1, dAlpha1);
00976 GetDerivativeVector(SecondDerivative, dL2, dPc2, dAlpha2);
00977
00978
00979 int in = matr_A.index1d(ia, im, il);
00980
00981 double dh1, dh11, dh12;
00982
00983
00984
00985
00986 dh1 = (L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL1][im+dPc1][ia+dAlpha1] + pc[il+dL1][im+dPc1][ia+dAlpha1] + alpha[il+dL1][im+dPc1][ia+dAlpha1]);
00987 dh11 = (L[il][im][ia] + pc[il][im][ia] + alpha[il][im][ia]) - (L[il+dL2][im+dPc2][ia+dAlpha2] + pc[il+dL2][im+dPc2][ia+dAlpha2] + alpha[il+dL2][im+dPc2][ia+dAlpha2]);
00988 dh12 = (L[il+dL1][im+dPc1][ia+dAlpha1] + pc[il+dL1][im+dPc1][ia+dAlpha1] + alpha[il+dL1][im+dPc1][ia+dAlpha1]) - (L[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2] + pc[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2] + alpha[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2]);
00989
00990 int id;
00991
00992
00993 id = matr_A.index1d(ia, im, il) - in;
00994
00995 matr_A[id][in] += +multiplicator / Jacobian[il][im][ia] / dh1 / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
00996
00997 id = matr_A.index1d(ia+dAlpha2, im+dPc2, il+dL2) - in;
00998 matr_A[id][in] += -multiplicator / Jacobian[il][im][ia] / dh1 / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
00999
01000 id = matr_A.index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
01001 matr_A[id][in] += -multiplicator / Jacobian[il][im][ia] / dh1 / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
01002
01003 id = matr_A.index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
01004 matr_A[id][in] += +multiplicator / Jacobian[il][im][ia] / dh1 / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
01005
01006 }
01007
01008
01009
01010 #define EE 1.e-19
01011 #define I(l,k) (l)*n+(k)
01012
01014 double max2(double a, double b) {
01015 return (a>b)?a:b;
01016 }
01018 void mult_vect2(double *af, double *vf, double *wf, int nf) {
01019 int k,l,n;
01020 n = nf;
01021 for (k = 0; k < nf; k++) {
01022 wf[k] = 0.;
01023 for(l=0; l < nf; l++)
01024 wf[k]+=af[( (k)*n+(l) )]*vf[l];
01025 }
01026 }
01027
01028
01030 void gauss_solve(double *a, double *b, int n) {
01031 int j, k, l, l_max;
01032 double a_max;
01033 double *dum;
01034
01035 dum = (double*) malloc(n*sizeof(double));
01036
01037 for (k = 0; k < (n-1); k++) {
01038
01039
01040 a_max = fabs ( a[I(k,k)] );
01041 l_max = k;
01042
01043 for(l = k + 1; l < n; l++)
01044
01045 if( fabs( a[I(l,k)] ) > a_max ) {
01046 a_max = fabs( a[I(l,k)] );
01047 l_max = l;
01048 }
01049 if( a_max < EE ) {
01050 printf("GAUSS: zero element string=%d size=%d \n", k, n);
01051 printf(" a(k,k)=%e \n",a[I(k,k)]);
01052 return;
01053 }
01054 if( l_max != k) {
01055 memcpy( dum, &(a[I(k,k)]), (n-k)*sizeof(double) );
01056 memcpy( &(a[I(k,k)]), &(a[I(l_max,k)]), (n-k)*sizeof(double) );
01057 memcpy( &(a[I(l_max,k)]), dum, (n-k)*sizeof(double) );
01058 a_max = b[k];
01059 b[k] = b[l_max];
01060 b[l_max] = a_max;
01061 }
01062
01063
01064 b[k] = b[k] / a[I(k,k)];
01065 for (l = k + 1; l < n; l++)
01066 a[( (k)*n+(l) )] = a[( (k)*n+(l) )] / a[I(k,k)];
01067 a[I(k,k)] = 1.;
01068 for (l = k + 1; l < n; l++) {
01069 for(j = k + 1; j < n; j++)
01070 a[I(l,j)] = a[I(l,j)] - a[I(k,j)] * a[I(l,k)];
01071 b[l] = b[l] - b[k] * a[I(l,k)];
01072 a[I(l,k)] = 0.;
01073 }
01074 }
01075
01076 free(dum);
01077 }
01078
01080 bool tridag(double a[], double b[], double c[], double r[], double u[], long n) {
01081 long j;
01082 double bet, *gam;
01083 gam = new double[n];
01084
01085 if (b[0] == 0.0) {
01086 throw error_msg("TRIDAG_MATRIX_ERROR", "tridag: error, b[0] = 0");
01087
01088
01089 }
01090
01091 u[0]=r[0]/(bet=b[0]);
01092
01093 for (j=1;j<=n-1;j++) {
01094 gam[j]=c[j-1]/bet;
01095 bet=b[j]-a[j]*gam[j];
01096
01097 if (bet == 0.0) {
01098 throw error_msg("TRIDAG_MATRIX_ERROR", "tridag: error, bet = 0");
01099
01100
01101 }
01102
01103 u[j]=(r[j]-a[j]*u[j-1])/bet;
01104 }
01105
01106 for (j=(n-2);j>=0;j--)
01107 u[j] -= gam[j+1]*u[j+1];
01108
01109 delete gam;
01110
01111 return true;
01112 }