23 #if defined(_WIN32) || defined(_WIN64)
24 #define _CRTDBG_MAP_ALLOC
41 #include "../Parameters/Parameters.h"
44 #include "../VariousFunctions/variousFunctions.h"
45 #include "../Exceptions/error.h"
48 #include "../Lapack/Include/cpplapack.h"
52 #include "../GMRES/itlin.h"
55 void matvec(
int n,
double *x,
double *b);
56 void preconr(
int n,
double *x,
double *b);
59 extern "C" void gmres(
int n,
double *y, MATVEC *
matvec,
60 PRECON *
preconr, PRECON *preconl,
double *b,
61 struct ITLIN_OPT *opt,
struct ITLIN_INFO *info);
65 #include "../GMRES2/GMRES2.h"
86 string lower_border_condition_type,
87 string upper_border_condition_type,
88 string approximationMethod)
92 double Coef1, Dm_05, Dp_05, tau_cur;
96 double *
B =
new double[nx];
99 double *h =
new double[nx];
100 double *h_bar =
new double[nx];
104 double coef_05m, coef_05p;
107 for (j = 1; j <= nx-1; j++) {
108 h[j] = x[j] - x[j-1];
110 for (j = 0; j <= nx-2; j++) {
111 h_bar[j] = (h[j] + h[j+1])/2;
114 for (j = 1; j <= nx-2; j++) {
116 if (approximationMethod ==
"AM_Split_LR") {
118 Coef1 = pow(x[j], n)/h_bar[j];
120 Coef1 = Coef1/
VF::G(x[j]);
121 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));
122 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));
124 coef_05m = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j-1]*pow(x[j-1], -n));
125 coef_05p = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j+1]*pow(x[j+1], -n));
127 A[j] = Coef1 * coef_05m / h[j];
128 C[j] = Coef1 * coef_05p / h[j+1];
130 }
else if (approximationMethod ==
"AM_Split_LR_oldway") {
132 Coef1 = pow(x[j], n);
134 Coef1 = Coef1/
VF::G(x[j]);
135 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));
136 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));
138 coef_05m = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j-1]*pow(x[j-1], -n));
139 coef_05p = 0.5*(Dxx[j]*pow(x[j], -n) + Dxx[j+1]*pow(x[j+1], -n));
141 A[j] = Coef1 * coef_05m / h[j] / h[j+1];
142 C[j] = Coef1 * coef_05p / h[j+1] / h[j];
144 }
else if (approximationMethod ==
"AM_Split_C") {
146 Coef1 = pow(x[j], n)/h_bar[j];
148 Coef1 = Coef1/
VF::G(x[j]);
150 Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
151 Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
153 A[j] = Coef1 * (Dm_05 * pow(x[j] - h[j]/2, -n) / h[j]);
154 C[j] = Coef1 * (Dp_05 * pow(x[j] + h[j+1]/2, -n) / h[j+1]);
156 A[j] = Coef1 * (Dm_05 *
VF::G(x[j] - h[j]/2) * pow(x[j] - h[j]/2, -n) / h[j]);
157 C[j] = Coef1 * (Dp_05 *
VF::G(x[j] + h[j+1]/2) * pow(x[j] + h[j+1]/2, -n) / h[j+1]);
160 }
else if (approximationMethod ==
"AM_Split_C_oldway") {
162 Coef1 = pow(x[j], n)/h_bar[j];
164 Coef1 = Coef1/
VF::G(x[j]);
166 Dm_05 = 0.5*(Dxx[j] + Dxx[j-1]);
167 Dp_05 = 0.5*(Dxx[j] + Dxx[j+1]);
169 A[j] = Coef1 * (Dm_05 * pow(x[j] - h[j]/2, -n) / h[j]);
170 C[j] = Coef1 * (Dp_05 * pow(x[j] + h[j]/2, -n) / h[j+1]);
172 A[j] = Coef1 * (Dm_05 *
VF::G(x[j] - h[j]/2) * pow(x[j] - h[j]/2, -n) / h[j]);
173 C[j] = Coef1 * (Dp_05 *
VF::G(x[j] + h[j]/2) * pow(x[j] + h[j]/2, -n) / h[j+1]);
177 throw error_msg(
"APPRERR",
"Unknown approximation type: %d", approximationMethod.c_str());
184 if ((x[j] < alc || x[j] >
VC::pi - alc) && (g_flag == 1)) {
190 B1[j] = B[j] - (1.0/dt + 1/tau_cur);
193 for (i = 1; i <= nx-2; i++) {
201 if (lower_border_condition_type ==
"BCT_CONSTANT_VALUE") {
202 A[i] = 0; B1[i] = 1; C[i] = 0;
203 }
else if (lower_border_condition_type ==
"BCT_CONSTANT_DERIVATIVE") {
204 A[i] = 0; B1[i] = -1; C[i] = 1;
206 throw error_msg(
"1D_DIFF_BOUNDARY",
"1d_universal_solver: unknown lower boundary type: %d", lower_border_condition_type.c_str());
214 if (upper_border_condition_type ==
"BCT_CONSTANT_VALUE") {
215 A[i] = 0; B1[i] = 1; C[i] = 0;
216 }
else if (upper_border_condition_type ==
"BCT_CONSTANT_DERIVATIVE") {
217 A[i] = -1; B1[i] = 1; C[i] = 0;
219 throw error_msg(
"1D_DIFF_BOUNDARY",
"1d_universal_solver: unknown lower boundary type: %d", upper_border_condition_type.c_str());
252 double *R1 =
new double[nx];
257 for (i = 1; i <= nx-2; i++) {
271 tridag(A, B1, C, R1, f, nx);
284 if (type ==
"BCT_CONSTANT_VALUE") {
285 Matrix_A[id0][in] = 1;
286 Matrix_A[id1][in] = 0;
287 }
else if (type ==
"BCT_CONSTANT_DERIVATIVE") {
288 Matrix_A[id0][in] = -1/dh;
289 Matrix_A[id1][in] = 1/dh;
291 throw error_msg(
"2D_DIFF_BOUNDARY",
"1d_universal_solver: unknown boundary type: %d", type.c_str());
305 int L_size,
int pc_size,
int alpha_size,
312 string L_lowerBoundaryCondition_calculationType,
313 string L_upperBoundaryCondition_calculationType,
314 string pc_lowerBoundaryCondition_calculationType,
315 string pc_upperBoundaryCondition_calculationType,
316 string alpha_lowerBoundaryCondition_calculationType,
317 string alpha_upperBoundaryCondition_calculationType,
323 double tau,
double tauLpp) {
329 int il, im, ia, id, in;
350 DiagMatrix::iterator it;
351 for (it = matr_A.begin(); it != matr_A.end(); it++) it->second = 0;
352 for (it = matr_B.begin(); it != matr_B.end(); it++) it->second = 0;
353 for (it = matr_C.begin(); it != matr_C.end(); it++) it->second = 0;
360 for (il = 0; il < L_size; il++) {
361 for (im = 0; im < pc_size; im++) {
362 for (ia = 0; ia < alpha_size; ia++) {
364 in = matr_A.
index1d(ia, im, il);
369 if (il == 0 && L_size >= 3) {
373 matr_C[0][in] = L_lowerBoundaryCondition[im][ia];
374 id = matr_A.
index1d(ia, im, il+1) - in;
375 dh = L[il+1][im][ia] - L[il][im][ia];
376 AddBoundary(matr_A, L_lowerBoundaryCondition_calculationType, in,
id, dh);
378 }
else if (il == L_size-1 && L_size >= 3) {
380 matr_C[0][in] = L_upperBoundaryCondition[im][ia];
381 id = matr_A.
index1d(ia, im, il-1) - in;
382 dh = L[il][im][ia] - L[il-1][im][ia];
383 AddBoundary(matr_A, L_upperBoundaryCondition_calculationType, in,
id, dh);
385 }
else if (im == 0 && pc_size >= 3) {
387 matr_C[0][in] = pc_lowerBoundaryCondition[il][ia];
388 id = matr_A.
index1d(ia, im+1, il) - in;
389 dh = pc[il][im+1][ia] - pc[il][im][ia];
390 AddBoundary(matr_A, pc_lowerBoundaryCondition_calculationType, in,
id, dh);
392 }
else if (im == pc_size-1 && pc_size >= 3) {
394 matr_C[0][in] = pc_upperBoundaryCondition[il][ia];
395 id = matr_A.
index1d(ia, im-1, il) - in;
396 dh = pc[il][im][ia] - pc[il][im-1][ia];
397 AddBoundary(matr_A, pc_upperBoundaryCondition_calculationType, in,
id, dh);
398 }
else if (ia == 0 && alpha_size >= 3) {
400 matr_C[0][in] = alpha_lowerBoundaryCondition[il][im];
401 id = matr_A.
index1d(ia+1, im, il) - in;
402 dh = alpha[il][im][ia+1] - alpha[il][im][ia];
403 AddBoundary(matr_A, alpha_lowerBoundaryCondition_calculationType, in,
id, dh);
405 }
else if (ia == alpha_size-1 && alpha_size >= 3) {
407 matr_C[0][in] = alpha_upperBoundaryCondition[il][im];
408 id = matr_A.
index1d(ia-1, im, il) - in;
409 dh = alpha[il][im][ia] - alpha[il][im][ia-1];
410 AddBoundary(matr_A, alpha_upperBoundaryCondition_calculationType, in,
id, dh);
417 matr_A[0][in] += 1.0/dt;
419 matr_B[0][in] += 1.0/dt;
425 "DT_L_left",
"DT_L_right",
426 L, pc, alpha, DLL, Jacobian, -0.5);
428 "DT_L_right",
"DT_L_left",
429 L, pc, alpha, DLL, Jacobian, -0.5);
434 if (L[il][im][ia] >= Lpp) {
436 "DT_PC_left",
"DT_PC_right",
437 L, pc, alpha, Dpcpc, Jacobian, -0.5);
439 "DT_PC_right",
"DT_PC_left",
440 L, pc, alpha, Dpcpc, Jacobian, -0.5);
443 "DT_PC_left",
"DT_PC_right",
444 L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
446 "DT_PC_right",
"DT_PC_left",
447 L, pc, alpha, DpcpcLpp, Jacobian, -0.5);
451 if (alpha_size >=3) {
460 if ( (alpha[il][im][ia] <=
VF::alc(L[il][im][0]) || alpha[il][im][ia] >=
VC::pi -
VF::alc(L[il][im][0]))) {
465 matr_A[0][in] += 1.0/taulc;
469 if (L[il][im][ia] >= Lpp) {
471 "DT_ALPHA_left",
"DT_ALPHA_right",
472 L, pc, alpha, Daa, Jacobian, -0.5);
474 "DT_ALPHA_right",
"DT_ALPHA_left",
475 L, pc, alpha, Daa, Jacobian, -0.5);
478 "DT_ALPHA_left",
"DT_ALPHA_right",
479 L, pc, alpha, DaaLpp, Jacobian, -0.5);
481 "DT_ALPHA_right",
"DT_ALPHA_left",
482 L, pc, alpha, DaaLpp, Jacobian, -0.5);
486 if (alpha_size >= 3 && pc_size >= 3) {
489 if (L[il][im][ia] >= Lpp) {
493 "DT_ALPHA_left",
"DT_PC_right",
494 L, pc, alpha, Dpca, Jacobian, -0.25);
496 "DT_ALPHA_right",
"DT_PC_left",
497 L, pc, alpha, Dpca, Jacobian, -0.25);
501 "DT_PC_left",
"DT_ALPHA_right",
502 L, pc, alpha, Dpca, Jacobian, -0.25);
504 "DT_PC_right",
"DT_ALPHA_left",
505 L, pc, alpha, Dpca, Jacobian, -0.25);
509 "DT_ALPHA_left",
"DT_PC_left",
510 L, pc, alpha, Dpca, Jacobian, -0.25);
512 "DT_ALPHA_right",
"DT_PC_right",
513 L, pc, alpha, Dpca, Jacobian, -0.25);
517 "DT_PC_left",
"DT_ALPHA_left",
518 L, pc, alpha, Dpca, Jacobian, -0.25);
520 "DT_PC_right",
"DT_ALPHA_right",
521 L, pc, alpha, Dpca, Jacobian, -0.25);
527 "DT_ALPHA_left",
"DT_PC_right",
528 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
530 "DT_ALPHA_right",
"DT_PC_left",
531 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
533 "DT_ALPHA_left",
"DT_PC_left",
534 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
536 "DT_ALPHA_right",
"DT_PC_right",
537 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
542 "DT_PC_left",
"DT_ALPHA_right",
543 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
545 "DT_PC_right",
"DT_ALPHA_left",
546 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
548 "DT_PC_left",
"DT_ALPHA_left",
549 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
551 "DT_PC_right",
"DT_ALPHA_right",
552 L, pc, alpha, DpcaLpp, Jacobian, -0.25);
578 void matvec(
int m_size,
double *x,
double *b)
580 DiagMatrix::iterator it;
582 for (in = 0; in < m_size; in++) {
592 if (ik >= 0 && ik < m_size) {
598 b[in] += it->second[in] * x[ik];
620 void jacobi(
int m_size,
double *z,
double *w)
623 for (in = 0; in < m_size; in++) {
625 w[in] = z[in] /
A_copy[0][in];
640 void SOR(
int m_size,
double *z,
double *w)
643 DiagMatrix::iterator it;
646 for (in = 0; in < m_size; in++) {
647 diag =
A_copy[0][in] / rel;
648 w[in] = z[in] / diag;
650 for (it =
A_copy.begin(); (it->first < 0 && it !=
A_copy.end()); it++) {
651 if (in + it->first >= 0 && it->first < 0) {
653 w[in] -= w[in + it->first] * it->second[in] / diag;
665 for (j=0;j<n;j++) b[j] = x[j];
677 int m_size = A[0].size_x;
683 DiagMatrix::iterator it;
686 std::clock_t start_time, current_time;
714 for (in = 0; in < m_size; in++) {
716 B_copy[in] = B_copy[in] / A0;
718 it->second[in] = it->second[in]/A0;
745 struct ITLIN_OPT *opt =
new struct ITLIN_OPT;
746 struct ITLIN_INFO *info =
new struct ITLIN_INFO;
748 TERM_CHECK termcheck = CheckOnRestart;
753 opt->termcheck = termcheck;
755 opt->errorlevel = None;
756 opt->monitorlevel = None;
757 opt->datalevel = None;
758 opt->errorfile = stdout;
759 opt->monitorfile = stdout;
760 opt->datafile = NULL;
763 opt->iterfile = NULL;
765 opt->miscfile = NULL;
767 start_time = std::clock();
771 gmres(m_size, &X[0], &
matvec, NULL, NULL, &B_copy[0], opt, info);
781 current_time = std::clock();
782 Output::echo(5,
"Done (%.2f sec)\n", (( current_time - start_time ) / (
double)CLOCKS_PER_SEC ));
788 fprintf(stdout,
"\n Return code: %6i\n",info->rcode);
789 fprintf(stdout,
" Iterations: %6i\n",info->iter);
790 fprintf(stdout,
" Matvec calls: %6i\n",info->nomatvec);
791 fprintf(stdout,
" Right Precon calls: %6i\n",info->noprecr);
792 fprintf(stdout,
" Left Precon calls: %6i\n",info->noprecl);
793 fprintf(stdout,
" Estimated residual reduction: %e\n",info->precision);
795 error = 0; err_max = 0;
796 for (in = 0; in < m_size; in++) {
799 for (it = A.begin(); it != A.end(); it++) {
800 if (in + it->first >= 0 && in + it->first < m_size) {
801 error -= it->second[in] * X[in + it->first];
804 err_max =
max(err_max, fabs(error));
815 const CPPL::dgbmatrix *
A;
825 CPPL::dcovector
solve(
const CPPL::dcovector &RHS)
const {
826 int in, m_size = RHS.l;
827 CPPL::dcovector res(m_size);
828 for (in = 0; in < m_size; in++) {
829 if ((*A)(in, in) != 0) {
831 res(in) = RHS(in) / (*A)(in, in);
841 const CPPL::dgbmatrix *
A;
850 CPPL::dcovector
solve(
const CPPL::dcovector &RHS)
const {
851 int in, it, m_size = RHS.l;
852 CPPL::dcovector res(m_size);
855 for (in = 0; in < m_size; in++) {
856 diag = (*A)(in, in) / rel;
857 res(in) = RHS(in) / diag;
860 for (it = (*A).kl; it <= (*A).ku; it++) {
862 res(in) -= res(it) * (*A)(in, it) / diag;
877 CPPL::dcovector X_lp(X.
size_x);
879 for (i = 0; i < X.
size_x; i++) X_lp(i) = X(i);
881 CPPL::dcovector B_lp(B.
size_x);
883 for (i = 0; i < B.
size_x; i++) B_lp(i) =
B(i);
885 long m_size = A[0].size_x;
886 DiagMatrix::iterator it;
888 long kl = -it->first;
892 CPPL::dgbmatrix A_lp(m_size, m_size, kl, ku);
895 int modelMatrixLineNumber;
896 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
897 for (
int modelMatrixDiagNumber = -kl; modelMatrixDiagNumber <= ku; modelMatrixDiagNumber++) {
898 if (modelMatrixDiagNumber + modelMatrixLineNumber >= 0 && modelMatrixDiagNumber + modelMatrixLineNumber < m_size) {
899 A_lp(modelMatrixLineNumber, modelMatrixDiagNumber + modelMatrixLineNumber) = 0;
906 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
907 for (it = A.begin(); it != A.end(); it++) {
909 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
911 A_lp(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
922 GMRES(A_lp, X_lp, B_lp,
923 M, H_lp, GMRES_parameters.
SOL_i_max, iter,
927 GMRES(A_lp, X_lp, B_lp,
928 M, H_lp, GMRES_parameters.
SOL_i_max, iter,
941 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
942 X[modelMatrixLineNumber] = X_lp(modelMatrixLineNumber);
956 std::clock_t start_time, current_time;
959 DiagMatrix::iterator it;
962 int modelMatrixLineNumber;
965 long m_size = A[0].size_x;
967 long kl = -it->first;
974 CPPL::dgbmatrix lap_AB(m_size, m_size, kl, ku);
975 CPPL::dgbmatrix lap_AB_res(m_size, m_size, kl, ku);
978 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
979 for (
int modelMatrixDiagNumber = -kl; modelMatrixDiagNumber <= ku; modelMatrixDiagNumber++) {
980 if (modelMatrixDiagNumber + modelMatrixLineNumber >= 0 && modelMatrixDiagNumber + modelMatrixLineNumber < m_size) {
981 lap_AB(modelMatrixLineNumber, modelMatrixDiagNumber + modelMatrixLineNumber) = 0;
989 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
990 for (it = A.begin(); it != A.end(); it++) {
995 if (modelMatrixLineNumber + it->first >= 0 && modelMatrixLineNumber + it->first < m_size) {
997 lap_AB(modelMatrixLineNumber, modelMatrixLineNumber + it->first) = it->second[modelMatrixLineNumber];
1004 lap_AB_res = lap_AB;
1007 CPPL::dcovector lap_B(m_size);
1008 CPPL::dcovector lap_B_res(m_size);
1009 for(i = 0; i < lap_B.l; i++){
1014 start_time = std::clock();
1018 lap_AB.dgbsv(lap_B);
1020 current_time = std::clock();
1021 Output::echo(5,
"done (%.2f sec)\n", (( current_time - start_time ) / (
double)CLOCKS_PER_SEC ));
1024 lap_B_res = lap_AB_res*lap_B - lap_B_res;
1028 for (i = 0; i < lap_B_res.l; i++) {
1029 max = (max>lap_B_res(i))?max:lap_B_res(i);
1035 for (modelMatrixLineNumber = 0; modelMatrixLineNumber < m_size; modelMatrixLineNumber++) {
1036 X[modelMatrixLineNumber] = lap_B(modelMatrixLineNumber);
1049 double err_max = 1e99;
1050 int n = A[0].size_x;
1055 DiagMatrix::iterator it;
1058 while(err_max > EE && nsteps < max_steps) {
1060 for (i = 0; i < n; i++) {
1062 for (it = A.begin(); it != A.end(); it++) {
1064 if (j >= 0 && j < n && j != i) {
1065 sum = sum + it->second[i]/A[0][i] * X[j];
1068 X[i] = (1 - omega) * X[i] - omega * ( sum - B[i]/A[0][i]);
1072 DiagMatrix::iterator it;
1074 for (j = 0; j < n; j++) {
1077 for (it = A.begin(); it != A.end(); it++) {
1078 if (j + it->first >= 0 && j + it->first < n) {
1079 error -= it->second[j] * X[j + it->first];
1082 err_max =
max(err_max, fabs(error));
1085 Output::echo(0,
"[%d]Error_norm = %e \n", nsteps, err_max);
1096 if (derivativeType ==
"DT_ALPHA_left") {
1098 }
else if (derivativeType ==
"DT_ALPHA_right") {
1100 }
else if (derivativeType ==
"DT_PC_left") {
1102 }
else if (derivativeType ==
"DT_PC_right") {
1104 }
else if (derivativeType ==
"DT_L_left") {
1106 }
else if (derivativeType ==
"DT_L_right") {
1109 throw error_msg(
"DERIVATIVE_TYPE",
"Unknown derivative approximation: %s\n", derivativeType.c_str());
1123 int il,
int im,
int ia,
1124 string FirstDerivative,
string SecondDerivative,
1127 double multiplicator) {
1129 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1130 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1141 int in = matr_A.
index1d(ia, im, il);
1143 double dh1, dh11, dh12;
1148 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]);
1149 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]);
1150 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]);
1160 double common_part = multiplicator / Jacobian[il][im][ia] / dh1 ;
1163 id = matr_A.
index1d(ia, im, il) - in;
1165 matr_A[id][in] += + common_part / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
1168 id = matr_A.
index1d(ia+dAlpha2, im+dPc2, il+dL2) - in;
1169 matr_A[id][in] += - common_part / dh11 * Dxx[il][im][ia] * Jacobian[il][im][ia];
1173 id = matr_A.
index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1174 matr_A[id][in] += - common_part / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1177 id = matr_A.
index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
1178 matr_A[id][in] += + common_part / dh12 * Dxx[il+dL1][im+dPc1][ia+dAlpha1] * Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1187 #define I(l,k) (l)*n+(k)
1197 for (k = 0; k < nf; k++) {
1199 for(l=0; l < nf; l++)
1200 wf[k]+=af[( (k)*n+(l) )]*vf[l];
1211 dum = (
double*) malloc(n*
sizeof(
double));
1213 for (k = 0; k < (n-1); k++) {
1216 a_max = fabs ( a[
I(k,k)] );
1219 for(l = k + 1; l < n; l++)
1221 if( fabs( a[
I(l,k)] ) > a_max ) {
1222 a_max = fabs( a[
I(l,k)] );
1226 printf(
"GAUSS: zero element string=%d size=%d \n", k, n);
1227 printf(
" a(k,k)=%e \n",a[
I(k,k)]);
1231 memcpy( dum, &(a[
I(k,k)]), (n-k)*
sizeof(
double) );
1232 memcpy( &(a[
I(k,k)]), &(a[
I(l_max,k)]), (n-k)*
sizeof(
double) );
1233 memcpy( &(a[
I(l_max,k)]), dum, (n-k)*
sizeof(
double) );
1240 b[k] = b[k] / a[
I(k,k)];
1241 for (l = k + 1; l < n; l++)
1242 a[( (k)*n+(l) )] = a[( (k)*n+(l) )] / a[
I(k,k)];
1244 for (l = k + 1; l < n; l++) {
1245 for(j = k + 1; j < n; j++)
1246 a[
I(l,j)] = a[
I(l,j)] - a[
I(k,j)] * a[
I(l,k)];
1247 b[l] = b[l] - b[k] * a[
I(l,k)];
1256 bool tridag(
double a[],
double b[],
double c[],
double r[],
double u[],
long n) {
1259 gam =
new double[n];
1262 throw error_msg(
"TRIDAG_MATRIX_ERROR",
"tridag: error, b[0] = 0");
1267 u[0]=r[0]/(bet=b[0]);
1269 for (j=1;j<=n-1;j++) {
1271 bet=b[j]-a[j]*gam[j];
1274 throw error_msg(
"TRIDAG_MATRIX_ERROR",
"tridag: error, bet = 0");
1279 u[j]=(r[j]-a[j]*u[j-1])/bet;
1282 for (j=(n-2);j>=0;j--)
1283 u[j] -= gam[j+1]*u[j+1];
1302 int L_size,
int pc_size,
int alpha_size,
1309 string L_lowerBoundaryCondition_calculationType,
1310 string L_upperBoundaryCondition_calculationType,
1311 string pc_lowerBoundaryCondition_calculationType,
1312 string pc_upperBoundaryCondition_calculationType,
1313 string alpha_lowerBoundaryCondition_calculationType,
1314 string alpha_upperBoundaryCondition_calculationType,
1320 double tau,
double tauLpp) {
1326 int il, im, ia, id, in;
1347 DiagMatrix::iterator it;
1348 for (it = matr_A.begin(); it != matr_A.end(); it++) it->second = 0;
1349 for (it = matr_B.begin(); it != matr_B.end(); it++) it->second = 0;
1350 for (it = matr_C.begin(); it != matr_C.end(); it++) it->second = 0;
1357 for (il = 0; il < L_size; il++) {
1358 for (im = 0; im < pc_size; im++) {
1359 for (ia = 0; ia < alpha_size; ia++) {
1361 in = matr_A.
index1d(ia, im, il);
1366 if (il == 0 && L_size >= 3) {
1370 matr_C[0][in] = L_lowerBoundaryCondition[im][ia];
1371 id = matr_A.
index1d(ia, im, il+1) - in;
1372 dh = L[il+1][im][ia] - L[il][im][ia];
1373 AddBoundary(matr_A, L_lowerBoundaryCondition_calculationType, in,
id, dh);
1375 }
else if (il == L_size-1 && L_size >= 3) {
1377 matr_C[0][in] = L_upperBoundaryCondition[im][ia];
1378 id = matr_A.
index1d(ia, im, il-1) - in;
1379 dh = L[il][im][ia] - L[il-1][im][ia];
1380 AddBoundary(matr_A, L_upperBoundaryCondition_calculationType, in,
id, dh);
1382 }
else if (im == 0 && pc_size >= 3) {
1384 matr_C[0][in] = pc_lowerBoundaryCondition[il][ia];
1385 id = matr_A.
index1d(ia, im+1, il) - in;
1386 dh = pc[il][im+1][ia] - pc[il][im][ia];
1387 AddBoundary(matr_A, pc_lowerBoundaryCondition_calculationType, in,
id, dh);
1389 }
else if (im == pc_size-1 && pc_size >= 3) {
1391 matr_C[0][in] = pc_upperBoundaryCondition[il][ia];
1392 id = matr_A.
index1d(ia, im-1, il) - in;
1393 dh = pc[il][im][ia] - pc[il][im-1][ia];
1394 AddBoundary(matr_A, pc_upperBoundaryCondition_calculationType, in,
id, dh);
1395 }
else if (ia == 0 && alpha_size >= 3) {
1397 matr_C[0][in] = alpha_lowerBoundaryCondition[il][im];
1398 id = matr_A.
index1d(ia+1, im, il) - in;
1399 dh = alpha[il][im][ia+1] - alpha[il][im][ia];
1400 AddBoundary(matr_A, alpha_lowerBoundaryCondition_calculationType, in,
id, dh);
1402 }
else if (ia == alpha_size-1 && alpha_size >= 3) {
1404 matr_C[0][in] = alpha_upperBoundaryCondition[il][im];
1405 id = matr_A.
index1d(ia-1, im, il) - in;
1406 dh = alpha[il][im][ia] - alpha[il][im][ia-1];
1407 AddBoundary(matr_A, alpha_upperBoundaryCondition_calculationType, in,
id, dh);
1414 matr_A[0][in] += 1.0/dt;
1416 matr_B[0][in] += 1.0/dt;
1422 "DT_L_left",
"DT_L_right",
1423 L, pc, alpha, DLL, Jacobian, -0.5);
1425 "DT_L_right",
"DT_L_left",
1426 L, pc, alpha, DLL, Jacobian, -0.5);
1431 if (L[il][im][ia] >= Lpp) {
1441 "DT_PC_right",
"DT_PC_left",
1442 L, pc, alpha, Dpcpc, Jacobian, -1.0);
1457 "DT_PC_right",
"DT_PC_left",
1458 L, pc, alpha, DpcpcLpp, Jacobian, -1.0);
1463 if (alpha_size >=3) {
1480 if (L[il][im][ia] >= Lpp) {
1490 "DT_ALPHA_right",
"DT_ALPHA_left",
1491 L, pc, alpha, Daa, Jacobian, -1.0);
1504 "DT_ALPHA_right",
"DT_ALPHA_left",
1505 L, pc, alpha, DaaLpp, Jacobian, -1.0);
1510 if (alpha_size >= 3 && pc_size >= 3) {
1513 if (L[il][im][ia] >= Lpp) {
1536 "DT_ALPHA_right",
"DT_PC_right",
1537 L, pc, alpha, Dpca, Jacobian, -1.0);
1539 "DT_ALPHA_left",
"DT_PC_right",
1540 L, pc, alpha, Dpca, Jacobian, +1.0);
1542 "DT_ALPHA_right",
"DT_PC_left",
1543 L, pc, alpha, Dpca, Jacobian, +1.0);
1545 "DT_ALPHA_left",
"DT_PC_left",
1546 L, pc, alpha, Dpca, Jacobian, -1.0);
1552 "DT_PC_right",
"DT_ALPHA_right",
1553 L, pc, alpha, Dpca, Jacobian, -1.0);
1555 "DT_PC_left",
"DT_ALPHA_right",
1556 L, pc, alpha, Dpca, Jacobian, +1.0);
1558 "DT_PC_right",
"DT_ALPHA_left",
1559 L, pc, alpha, Dpca, Jacobian, +1.0);
1561 "DT_PC_left",
"DT_ALPHA_left",
1562 L, pc, alpha, Dpca, Jacobian, -1.0);
1603 "DT_ALPHA_right",
"DT_PC_right",
1604 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1606 "DT_ALPHA_right",
"DT_PC_left",
1607 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1609 "DT_ALPHA_left",
"DT_PC_right",
1610 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1612 "DT_ALPHA_left",
"DT_PC_left",
1613 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1619 "DT_PC_right",
"DT_ALPHA_right",
1620 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1622 "DT_PC_right",
"DT_ALPHA_left",
1623 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1625 "DT_PC_left",
"DT_ALPHA_right",
1626 L, pc, alpha, DpcaLpp, Jacobian, +1.0);
1628 "DT_PC_left",
"DT_ALPHA_left",
1629 L, pc, alpha, DpcaLpp, Jacobian, -1.0);
1658 int il,
int im,
int ia,
1659 string FirstDerivative,
string SecondDerivative,
1662 double multiplicator) {
1664 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1665 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1675 int in = matr_A.
index1d(ia, im, il);
1681 double Dxx_0, Dxx_1, Dxx_2, Dxx_11, Dxx_21;
1682 double J_0, J_1, J_2;
1685 dh_1 = abs((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]));
1686 dh_2 = abs((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]));
1688 Dxx_0 = Dxx[il][im][ia];
1689 Dxx_1 = Dxx[il+dL1][im+dPc1][ia+dAlpha1];
1690 Dxx_2 = Dxx[il+dL2][im+dPc2][ia+dAlpha2];
1691 Dxx_21 = Dxx[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2];
1693 J_0 = Jacobian[il][im][ia];
1694 J_1 = Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1695 J_2 = Jacobian[il+dL2][im+dPc2][ia+dAlpha2];
1700 common = multiplicator / J_0 / (dh_1+dh_2) * 2;
1702 id = matr_A.
index1d(ia, im, il) - in;
1704 matr_A[id][in] += -common / dh_1 * (Dxx_0+Dxx_1)/2 * (J_0+J_1)/2;
1706 id = matr_A.
index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1707 matr_A[id][in] += +common / dh_1 * (Dxx_0+Dxx_1)/2 * (J_0+J_1)/2;
1709 id = matr_A.
index1d(ia+dAlpha2, im+dPc2, il+dL2) - in;
1710 matr_A[id][in] += +common / dh_2 * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2;
1712 id = matr_A.
index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
1713 matr_A[id][in] += -common / dh_2 * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2;
1721 int il,
int im,
int ia,
1722 string FirstDerivative,
string SecondDerivative,
1725 double multiplicator) {
1727 int dL1 = 0, dPc1 = 0, dAlpha1 = 0;
1728 int dL2 = 0, dPc2 = 0, dAlpha2 = 0;
1738 int in = matr_A.
index1d(ia, im, il);
1743 double dh_1, dh_2, dh_11, dh_21;
1744 double Dxx_0, Dxx_1, Dxx_2, Dxx_11, Dxx_21;
1745 double Dxx_c1, Dxx_c2, Dxx_c3, Dxx_c4;
1746 double J_0, J_1, J_2, J_11, J_21;
1750 dh_1 = abs((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]));
1751 dh_2 = abs((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]));
1752 dh_11 = abs((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]));
1753 dh_21 = abs((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]));
1755 Dxx_0 = Dxx[il][im][ia];
1756 Dxx_1 = Dxx[il+dL1][im+dPc1][ia+dAlpha1];
1757 Dxx_2 = Dxx[il+dL2][im+dPc2][ia+dAlpha2];
1758 Dxx_11 = Dxx[il-dL1][im-dPc1][ia-dAlpha1];
1759 Dxx_21 = Dxx[il-dL2][im-dPc2][ia-dAlpha2];
1762 Dxx_c1 = Dxx[il+dL1+dL2][im+dPc1+dPc2][ia+dAlpha1+dAlpha2];
1763 Dxx_c2 = Dxx[il+dL1-dL2][im+dPc1-dPc2][ia+dAlpha1-dAlpha2];
1764 Dxx_c3 = Dxx[il-dL1+dL2][im-dPc1+dPc2][ia-dAlpha1+dAlpha2];
1765 Dxx_c4 = Dxx[il-dL1-dL2][im-dPc1-dPc2][ia-dAlpha1-dAlpha2];
1768 J_0 = Jacobian[il][im][ia];
1769 J_1 = Jacobian[il+dL1][im+dPc1][ia+dAlpha1];
1770 J_2 = Jacobian[il+dL2][im+dPc2][ia+dAlpha2];
1771 J_11 = Jacobian[il-dL1][im-dPc1][ia-dAlpha1];
1772 J_21 = Jacobian[il-dL2][im-dPc2][ia-dAlpha2];
1777 common = multiplicator / J_0 / (dh_1+dh_11) / (dh_2+dh_21);
1779 id = matr_A.
index1d(ia+dAlpha1, im+dPc1, il+dL1) - in;
1781 matr_A[id][in] += +common * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2;
1782 id = matr_A.
index1d(ia+dAlpha1+dAlpha2, im+dPc1+dPc2, il+dL1+dL2) - in;
1783 matr_A[id][in] += +common * (Dxx_0+Dxx_2)/2 * (J_0+J_2)/2;