diff options
author | hawke <> | 2003-11-19 09:40:00 +0000 |
---|---|---|
committer | hawke <> | 2003-11-19 09:40:00 +0000 |
commit | 271adecccf69c10a2b0d12a6e0a482ae24c36e51 (patch) | |
tree | 35bff7358d829cf1b8672392a1722e093d6c0326 /CarpetDev/CarpetCG/src | |
parent | bbbd74e612d0dbfba368ed9212360007de552683 (diff) |
Initial revision
darcs-hash:20031119094047-58737-ef9911a12e56712a462cb6e555f79d0fa4a7dd6a.gz
Diffstat (limited to 'CarpetDev/CarpetCG/src')
-rw-r--r-- | CarpetDev/CarpetCG/src/CG.cc | 934 |
1 files changed, 483 insertions, 451 deletions
diff --git a/CarpetDev/CarpetCG/src/CG.cc b/CarpetDev/CarpetCG/src/CG.cc index e27e465be..ea2a9f3d9 100644 --- a/CarpetDev/CarpetCG/src/CG.cc +++ b/CarpetDev/CarpetCG/src/CG.cc @@ -1,4 +1,4 @@ -/* $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/CG.cc,v 1.4 2004/01/25 14:57:28 schnetter Exp $ */ +/* $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/CG.cc,v 1.1 2003/11/19 10:40:47 hawke Exp $ */ #include <cassert> #include <cmath> @@ -55,14 +55,13 @@ namespace CarpetCG { int ierr; - vector<int> nboundaryzones(2*dim); + int * nboundaryzones; CCTK_REAL factor; - vector<CCTK_REAL> factors; - - vector<int> fromindex; // Can't pass things through CallLocalFunction. - vector<int> toindex; // Instead assume everything is going from a GF to a GF. + CCTK_REAL * restrict factors; + int * fromindex; // Can't pass things through CallLocalFunction. + int * toindex; // Instead assume everything is going from a GF to a GF. CCTK_REAL realconstant; // With only one required constant CCTK_REAL realoutput; // And one output CCTK_REAL realoutput_count; // And one _more_ output @@ -629,26 +628,45 @@ namespace CarpetCG { assert (calcres); assert (applybounds); - nelems = Util_TableGetIntArray - (options_table, 2*dim, &(common::nboundaryzones[0]), "nboundaryzones"); - if (nelems == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - for (d=0; d<dim; ++d) { - common::nboundaryzones[2*d ] = -1; - common::nboundaryzones[2*d+1] = -1; - } - } else if (nelems != 2*dim) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Options table key \"nboundaryzones\" is not an integer array of length %d", 2*dim); - return -1; + if (!common::nboundaryzones) + { + common::nboundaryzones = + (int*)malloc(2 * dim * sizeof *common::nboundaryzones); + nelems = Util_TableGetIntArray + (options_table, 2*dim, common::nboundaryzones, "nboundaryzones"); + if (nelems == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + for (d=0; d<dim; ++d) { + common::nboundaryzones[2*d ] = -1; + common::nboundaryzones[2*d+1] = -1; + } + } else if (nelems != 2*dim) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Options table key \"nboundaryzones\" is not an integer array of length %d", 2*dim); + return -1; + } + } - (common::fromindex).reserve(cg_maxsolvevars); - (common::toindex).reserve(cg_maxsolvevars); - (common::factors).reserve(nvars); - common::factor = 0.5; - for (int i = 0; i < nvars; i++) + + if (!common::fromindex) { - common::factors[i] = 1.0; + common::fromindex = (int*)malloc (nvars * sizeof *common::fromindex); + } + assert (common::fromindex); + if (!common::toindex) + { + common::toindex = (int*)malloc (nvars * sizeof *common::toindex); + } + assert (common::toindex); + + if (!common::factors) + { + common::factor = 0.5; + common::factors = (CCTK_REAL*)malloc (nvars * sizeof *common::factors); + for (int i = 0; i < nvars; i++) + { + common::factors[i] = 1.0; + } } maxiters = 1000; @@ -680,280 +698,343 @@ namespace CarpetCG { assert (minerror >= 0); assert (sminerror >= 0); assert (stepsize > 0); - - // Switch to global mode - BEGIN_GLOBAL_MODE(cctkGH) { - - // Fill common block - common::var = var; - common::res = res; - common::nvars = nvars; - common::options_table = options_table; - common::calcres = calcres; - common::applybnds = applybounds; - common::userdata = userdata; - - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Solving through nonlinear conjugate gradient iterations"); - } - - - - /* - * Literature: - * - * Jonathan Richard Shewchuk, "An Introduction to the Conjugate - * Gradient Method Without the Agonizing Pain", Technical Report - * CMU-CS-94-125, Aug. 1994 - * - * Available from http://www-2.cs.cmu.edu/~jrs/jrspapers.html - */ + + // Fill common block + common::var = var; + common::res = res; + common::nvars = nvars; + common::options_table = options_table; + common::calcres = calcres; + common::applybnds = applybounds; + common::userdata = userdata; + + // Check state + // (want global mode, or at least level mode) + assert (reflevel==-1 || component==-1); + const int saved_reflevel = reflevel; + const int saved_mglevel = mglevel; + set_mglevel ((cGH *)cctkGH, -1); + set_reflevel ((cGH *)cctkGH, -1); + + + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Solving through nonlinear conjugate gradient iterations"); + } + + + + /* + * Literature: + * + * Jonathan Richard Shewchuk, "An Introduction to the Conjugate + * Gradient Method Without the Agonizing Pain", Technical Report + * CMU-CS-94-125, Aug. 1994 + * + * Available from http://www-2.cs.cmu.edu/~jrs/jrspapers.html + */ - /* - * Notation: - * - * here there meaning - * - * iter i current cg iteration - * siter j current secant iteration - * iter2 k current cg iteration since last restart - * maxiters2 n PARA cg iterations before restart - * alpha secant step size - * beta cg step size - * delta cg distance - * epsilon PARA secant error tolerance (straight epsilon) - * minerror vareps PARA cg error tolerance (curved epsilon) - * eta secant distance - * stepsize sigma PARA secand method step parameter - * - * dir d cg direction - * -calcres f' PARA nonlinear operator - * res r cg residual - * prs s preconditioned cg residual - * var x INIT unknown variable - * - * -wgt M preconditioning matrix (not in this version) - */ + /* + * Notation: + * + * here there meaning + * + * iter i current cg iteration + * siter j current secant iteration + * iter2 k current cg iteration since last restart + * maxiters2 n PARA cg iterations before restart + * alpha secant step size + * beta cg step size + * delta cg distance + * epsilon PARA secant error tolerance (straight epsilon) + * minerror vareps PARA cg error tolerance (curved epsilon) + * eta secant distance + * stepsize sigma PARA secand method step parameter + * + * dir d cg direction + * -calcres f' PARA nonlinear operator + * res r cg residual + * prs s preconditioned cg residual + * var x INIT unknown variable + * + * -wgt M preconditioning matrix (not in this version) + */ - /* - * Algorithm: - * (Preconditioned nonlinear conjugate gradients with secant and - * Polak-Ribière) - * - * 01. i <= 0 - * 02. k <= 0 - * 03. r <= - f'(x) - * 04. calculate M \approx f''(x) - * 05. s <= M^-1 r - * 06. d <= s - * 07. delta_new <= r^T d - * 08. delta_0 <= delta_new - * 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO - * 10. j <= 0 - * 11. delta_d <= d^T d - * 12. alpha <= - sigma_0 - * 13. eta_prev <= [f'(x + sigma_0 d)]^T d - * 14. DO - * 15. eta <= [f'(x)]^T d - * 16. alpha <= alpha (eta) / (eta_prev - eta) - * 17. x <= x + alpha d - * 18. eta_prev <= eta - * 19. j <= j + 1 - * 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 - * 21. r <= - f'(x) - * 22. delta_old <= delta_new - * 23. delta_mid <= r^T s - * 24. calculate M \approx f''(x) - * 25. s <= M^-1 r - * 26. delta_new <= r^T s - * 27. beta <= (delta_new - delta_mid) / (delta_old) - * 28. k <= k + 1 - * 29. IF k = n OR beta <= 0 THEN - * 30. d <= s - * 31. k <= 0 - * 32. ELSE - * 33. d <= s + beta d - * 34. i <= i + 1 - */ + /* + * Algorithm: + * (Preconditioned nonlinear conjugate gradients with secant and + * Polak-Ribière) + * + * 01. i <= 0 + * 02. k <= 0 + * 03. r <= - f'(x) + * 04. calculate M \approx f''(x) + * 05. s <= M^-1 r + * 06. d <= s + * 07. delta_new <= r^T d + * 08. delta_0 <= delta_new + * 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO + * 10. j <= 0 + * 11. delta_d <= d^T d + * 12. alpha <= - sigma_0 + * 13. eta_prev <= [f'(x + sigma_0 d)]^T d + * 14. DO + * 15. eta <= [f'(x)]^T d + * 16. alpha <= alpha (eta) / (eta_prev - eta) + * 17. x <= x + alpha d + * 18. eta_prev <= eta + * 19. j <= j + 1 + * 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 + * 21. r <= - f'(x) + * 22. delta_old <= delta_new + * 23. delta_mid <= r^T s + * 24. calculate M \approx f''(x) + * 25. s <= M^-1 r + * 26. delta_new <= r^T s + * 27. beta <= (delta_new - delta_mid) / (delta_old) + * 28. k <= k + 1 + * 29. IF k = n OR beta <= 0 THEN + * 30. d <= s + * 31. k <= 0 + * 32. ELSE + * 33. d <= s + beta d + * 34. i <= i + 1 + */ - // Pointers to grid variables and temporary storage - vector<int> varptrs(nvars); - vector<int> resptrs(nvars); - // CCTK_REAL * restrict * restrict prsptrs; - // CCTK_REAL * restrict * restrict dirptrs; + // Pointers to grid variables and temporary storage + int varptrs[nvars]; + int resptrs[nvars]; + // CCTK_REAL * restrict * restrict prsptrs; + // CCTK_REAL * restrict * restrict dirptrs; - // Get storage pointers - for (n=0; n<nvars; ++n) { - varptrs[n] = var[n]; - resptrs[n] = res[n]; - } + // Get storage pointers + for (n=0; n<nvars; ++n) { + varptrs[n] = var[n]; + resptrs[n] = res[n]; + } - // Allocate temporary memory - // prsptrs = malloc (nvars * sizeof *prsptrs); - // dirptrs = malloc (nvars * sizeof *prsptrs); - assert (cg_maxsolvevars >= nvars); - vector<int> prsptrs(cg_maxsolvevars); - vector<int> dirptrs(cg_maxsolvevars); - for (n=0; n<nvars; ++n) { - char varname[1000]; - sprintf(varname,"CarpetCG::prsvars[%d]",n); - prsptrs[n] = CCTK_VarIndex(varname); - assert ((prsptrs[n] >= 0)&&(prsptrs[n] < CCTK_NumVars())); - sprintf(varname,"CarpetCG::dirvars[%d]",n); - dirptrs[n] = CCTK_VarIndex(varname); - assert ((dirptrs[n] >= 0)&&(dirptrs[n] < CCTK_NumVars())); - } + // Allocate temporary memory + // prsptrs = malloc (nvars * sizeof *prsptrs); + // dirptrs = malloc (nvars * sizeof *prsptrs); + assert (cg_maxsolvevars >= nvars); + int prsptrs[cg_maxsolvevars]; + int dirptrs[cg_maxsolvevars]; + assert (prsptrs); + assert (dirptrs); + for (n=0; n<nvars; ++n) { + char varname[1000]; + sprintf(varname,"CarpetCG::prsvars[%d]",n); + prsptrs[n] = CCTK_VarIndex(varname); + assert ((prsptrs[n] >= 0)&&(prsptrs[n] < CCTK_NumVars())); + sprintf(varname,"CarpetCG::dirvars[%d]",n); + dirptrs[n] = CCTK_VarIndex(varname); + assert ((dirptrs[n] >= 0)&&(dirptrs[n] < CCTK_NumVars())); + } - nexttime = time(0); + nexttime = time(0); - /* 01. i <= 0 */ - iter = 0; + /* 01. i <= 0 */ + iter = 0; - /* 02. k <= 0 */ - iter2 = 0; + /* 02. k <= 0 */ + iter2 = 0; - /* 03. r <= - f'(x) */ - /* 04. calculate M \approx f''(x) */ - common::ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_calcres); - ierr = common::ierr; - assert (! ierr); + /* 03. r <= - f'(x) */ + /* 04. calculate M \approx f''(x) */ + common::ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_calcres); + ierr = common::ierr; + assert (! ierr); + + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + } + CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); + + /* 05. s <= M^-1 r */ + // No preconditioning + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + common::toindex[n] = prsptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + // cout << "dirptrs " << dirptrs[0] << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); + + // cout << "dirptrs " << dirptrs[0] << endl; + + /* 06. d <= s */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = dirptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + // cout << "dirptrs " << dirptrs[0] << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_copy); + + // cout << "dirptrs " << dirptrs[0] << endl; + + /* 07. delta_new <= r^T d */ + // output (cctkGH, var, res, wgt, nvars, iter); + common::realoutput_count = 0; + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + // cout << "dirptrs " << dirptrs[0] << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_new = common::realoutput; + gsize = common::realoutput_count; + // cout << "delta_new " << delta_new << " gsize " << gsize << endl; + + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + epsilon = common::realoutput; +// cout << "epsilon " << epsilon << endl; + + /* 08. delta_0 <= delta_new */ + delta_0 = delta_new; + + /* 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO */ + while (iter < maxiters && epsilon / (nvars * gsize) > pow(minerror,2)) { + + if (verbose || veryverbose) { + currenttime = time(0); + if (veryverbose || (iter % outevery == 0 && currenttime >= nexttime)) { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d (%d since restart): residual is %g (%g,%d,%g)", + iter, iter2, (double)sqrt(epsilon / (nvars * gsize)),epsilon,nvars,gsize); + if (outeveryseconds > 0) { + while (nexttime <= currenttime) nexttime += outeveryseconds; + } + } + } + + /* 10. j <= 0 */ + siter = 0; + /* 11. delta_d <= d^T d */ for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = dirptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; + common::toindex[n] = -1; } - CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_d = common::realoutput; + +// cout << "delta_d " << delta_d << endl; + + /* 12. alpha <= - sigma_0 */ + alpha = - stepsize; + sum_alpha = alpha; + do_abort = 0; - /* 05. s <= M^-1 r */ - // No preconditioning + /* 13. eta_prev <= [f'(x + sigma_0 d)]^T d */ for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = prsptrs[n]; + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = varptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } + common::realconstant = -alpha; +// cout << "Add scaled, const " << -alpha << endl; + + CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); + common::ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_applybnds); + assert (! ierr); - // cout << "dirptrs " << dirptrs[0] << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); - - // cout << "dirptrs " << dirptrs[0] << endl; - - /* 06. d <= s */ + ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_calcres); + assert (! ierr); + for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = dirptrs[n]; + common::fromindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; - common::toindex[n] = -1; } + CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - // cout << "dirptrs " << dirptrs[0] << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_copy); - - // cout << "dirptrs " << dirptrs[0] << endl; - - /* 07. delta_new <= r^T d */ - // output (cctkGH, var, res, wgt, nvars, iter); - common::realoutput_count = 0; for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; + common::fromindex[n] = dirptrs[n]; common::toindex[n] = resptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - - // cout << "dirptrs " << dirptrs[0] << endl; - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_new = common::realoutput; - gsize = common::realoutput_count; - // cout << "delta_new " << delta_new << " gsize " << gsize << endl; - + eta = - common::realoutput; + +// cout << "eta " << eta << endl; + for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = resptrs[n]; + common::fromindex[n] = dirptrs[n]; + common::toindex[n] = varptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - epsilon = common::realoutput; - // cout << "epsilon " << epsilon << endl; - - /* 08. delta_0 <= delta_new */ - delta_0 = delta_new; - - /* 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO */ - while (iter < maxiters && epsilon / (nvars * gsize) > pow(minerror,2)) { - - if (verbose || veryverbose) { - currenttime = time(0); - if (veryverbose || (iter % outevery == 0 && currenttime >= nexttime)) { - CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d (%d since restart): residual is %g (%g,%d,%g)", - iter, iter2, (double)sqrt(epsilon / (nvars * gsize)),epsilon,nvars,gsize); - if (outeveryseconds > 0) { - while (nexttime <= currenttime) nexttime += outeveryseconds; - } - } - } - - /* 10. j <= 0 */ - siter = 0; + common::realconstant = alpha; +// cout << "Add scaled, const " << alpha << endl; - /* 11. delta_d <= d^T d */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = dirptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_d = common::realoutput; - - // cout << "delta_d " << delta_d << endl; + CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); + ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_applybnds); + assert (! ierr); - /* 12. alpha <= - sigma_0 */ - alpha = - stepsize; - sum_alpha = alpha; - do_abort = 0; + eta_prev = eta; - /* 13. eta_prev <= [f'(x + sigma_0 d)]^T d */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = varptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - common::realconstant = -alpha; - // cout << "Add scaled, const " << -alpha << endl; - - CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); - common::ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_applybnds); - assert (! ierr); + /* 14. DO */ + do { + if (veryverbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Secant iteration %d: step size is %g, orthogonality is %g", + siter, + (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); + } + + /* 15. eta <= [f'(x)]^T d */ ierr = 0; CallLocalFunction((cGH *)cctkGH, common::call_calcres); assert (! ierr); @@ -965,7 +1046,7 @@ namespace CarpetCG { common::fromindex[n] = -1; } CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - + for (int n = 0; n < nvars; n++) { common::fromindex[n] = dirptrs[n]; common::toindex[n] = resptrs[n]; @@ -974,10 +1055,34 @@ namespace CarpetCG { common::fromindex[n] = -1; common::toindex[n] = -1; } + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); eta = - common::realoutput; - // cout << "eta " << eta << endl; +// cout << "eta " << eta << endl; + + /* 16. alpha <= alpha (eta) / (eta_prev - eta) */ + alpha_old = alpha; + alpha *= eta / (eta_prev - eta); + sum_alpha += alpha; + if (veryverbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Changing step size, iteration %d: was %g, now %g (%g, %g)", + siter, alpha_old, alpha, eta, eta_prev); + } + assert (sum_alpha > - maxstepsize); + if (sum_alpha > maxstepsize) { + alpha = maxstepsize - alpha_old; + do_abort = 1; + if (verbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Secant iteration %d: limiting total step size", siter); + } + } + + /* 17. x <= x + alpha d */ for (int n = 0; n < nvars; n++) { common::fromindex[n] = dirptrs[n]; @@ -988,249 +1093,176 @@ namespace CarpetCG { common::toindex[n] = -1; } common::realconstant = alpha; - // cout << "Add scaled, const " << alpha << endl; +// cout << "Add scaled, const " << alpha << endl; CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); ierr = 0; CallLocalFunction((cGH *)cctkGH, common::call_applybnds); assert (! ierr); - + + /* 18. eta_prev <= eta */ eta_prev = eta; + + /* 19. j <= j + 1 */ + ++ siter; + + /* 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 */ + } while (siter < smaxiters + && pow(alpha,2) * delta_d > pow(sminerror,2) + && !do_abort); - /* 14. DO */ - do { + if (veryverbose) { + CCTK_VInfo + (CCTK_THORNSTRING, + " Secant iteration %d: step size is %g, orthogonality is %g", + siter, + (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); + } - if (veryverbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Secant iteration %d: step size is %g, orthogonality is %g", - siter, - (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); - } - - /* 15. eta <= [f'(x)]^T d */ - ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_calcres); - assert (! ierr); - - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - } - CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); - - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } + /* 21. r <= - f'(x) */ + /* 24. calculate M \approx f''(x) */ + ierr = 0; + CallLocalFunction((cGH *)cctkGH, common::call_calcres); + assert (! ierr); - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - eta = - common::realoutput; + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + } + CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); + + /* 23. delta_mid <= r^T s */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } - // cout << "eta " << eta << endl; + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_mid = common::realoutput; +// cout << "delta_mid " << delta_mid << endl; - /* 16. alpha <= alpha (eta) / (eta_prev - eta) */ - alpha_old = alpha; - alpha *= eta / (eta_prev - eta); - sum_alpha += alpha; - if (veryverbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Changing step size, iteration %d: was %g, now %g (%g, %g)", - siter, alpha_old, alpha, eta, eta_prev); - } - assert (sum_alpha > - maxstepsize); - if (sum_alpha > maxstepsize) { - alpha = maxstepsize - alpha_old; - do_abort = 1; - if (verbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Secant iteration %d: limiting total step size", siter); - } - } + /* 25. s <= M^-1 r */ + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + common::toindex[n] = prsptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } - /* 17. x <= x + alpha d */ + CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = dirptrs[n]; - common::toindex[n] = varptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - common::realconstant = alpha; - // cout << "Add scaled, const " << alpha << endl; + /* 22. delta_old <= delta_new */ + delta_old = delta_new; + + /* 26. delta_new <= r^T s */ +// output (cctkGH, var, res, wgt, nvars, iter+1); + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = prsptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } - CallLocalFunction((cGH *)cctkGH, common::call_add_scaled); - ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_applybnds); - assert (! ierr); - - /* 18. eta_prev <= eta */ - eta_prev = eta; - - /* 19. j <= j + 1 */ - ++ siter; - - /* 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 */ - } while (siter < smaxiters - && pow(alpha,2) * delta_d > pow(sminerror,2) - && !do_abort); + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + delta_new = common::realoutput; +// cout << "delta_new " << delta_new << endl; + + for (int n = 0; n < nvars; n++) { + common::fromindex[n] = resptrs[n]; + common::toindex[n] = resptrs[n]; + } + for (int n = nvars; n < cg_maxsolvevars; n++) { + common::fromindex[n] = -1; + common::toindex[n] = -1; + } + + CallLocalFunction((cGH *)cctkGH, common::call_dot_product); + epsilon = common::realoutput; +// cout << "epsilon " << epsilon << endl; - if (veryverbose) { - CCTK_VInfo - (CCTK_THORNSTRING, - " Secant iteration %d: step size is %g, orthogonality is %g", - siter, - (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize))); - } + /* 27. beta <= (delta_new - delta_mid) / (delta_old) */ + beta = (delta_new - delta_mid) / delta_old; - /* 21. r <= - f'(x) */ - /* 24. calculate M \approx f''(x) */ - ierr = 0; - CallLocalFunction((cGH *)cctkGH, common::call_calcres); - assert (! ierr); - - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - } - CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign); + /* 28. k <= k + 1 */ + ++ iter2; - /* 23. delta_mid <= r^T s */ + /* 29. IF k = n OR beta <= 0 THEN */ + if (iter2 >= maxiters2 || beta <= 0 || do_abort) { + /* Restart */ + + /* 30. d <= s */ for (int n = 0; n < nvars; n++) { common::fromindex[n] = prsptrs[n]; - common::toindex[n] = resptrs[n]; + common::toindex[n] = dirptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_mid = common::realoutput; - // cout << "delta_mid " << delta_mid << endl; + CallLocalFunction((cGH *)cctkGH, common::call_copy); - /* 25. s <= M^-1 r */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = prsptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } + /* 31. k <= 0 */ + iter2 = 0; - CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner); - - /* 22. delta_old <= delta_new */ - delta_old = delta_new; - - /* 26. delta_new <= r^T s */ - // output (cctkGH, var, res, wgt, nvars, iter+1); + /* 32. ELSE */ + } else { + /* Continue with CG iterations */ + + /* 33. d <= s + beta d */ for (int n = 0; n < nvars; n++) { common::fromindex[n] = prsptrs[n]; - common::toindex[n] = resptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - delta_new = common::realoutput; - // cout << "delta_new " << delta_new << endl; - - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = resptrs[n]; - common::toindex[n] = resptrs[n]; + common::toindex[n] = dirptrs[n]; } for (int n = nvars; n < cg_maxsolvevars; n++) { common::fromindex[n] = -1; common::toindex[n] = -1; } - - CallLocalFunction((cGH *)cctkGH, common::call_dot_product); - epsilon = common::realoutput; - // cout << "epsilon " << epsilon << endl; - - /* 27. beta <= (delta_new - delta_mid) / (delta_old) */ - beta = (delta_new - delta_mid) / delta_old; - - /* 28. k <= k + 1 */ - ++ iter2; - - /* 29. IF k = n OR beta <= 0 THEN */ - if (iter2 >= maxiters2 || beta <= 0 || do_abort) { - /* Restart */ - - /* 30. d <= s */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = dirptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - CallLocalFunction((cGH *)cctkGH, common::call_copy); - - /* 31. k <= 0 */ - iter2 = 0; - - /* 32. ELSE */ - } else { - /* Continue with CG iterations */ - - /* 33. d <= s + beta d */ - for (int n = 0; n < nvars; n++) { - common::fromindex[n] = prsptrs[n]; - common::toindex[n] = dirptrs[n]; - } - for (int n = nvars; n < cg_maxsolvevars; n++) { - common::fromindex[n] = -1; - common::toindex[n] = -1; - } - common::realconstant = beta; - // cout << "beta " << beta << endl; + common::realconstant = beta; + // cout << "beta " << beta << endl; - CallLocalFunction((cGH *)cctkGH, common::call_scale_and_add); + CallLocalFunction((cGH *)cctkGH, common::call_scale_and_add); - } /* if iter2 */ + } /* if iter2 */ /* 34. i <= i + 1 */ - ++ iter; + ++ iter; - } /* for iter */ + } /* for iter */ - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d (%d since restart): residual is %g", - iter, iter2, (double)sqrt(epsilon / (nvars * gsize))); - } + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d (%d since restart): residual is %g", + iter, iter2, (double)sqrt(epsilon / (nvars * gsize))); + } /* Free temporary memory */ - // for (n=0; n<nvars; ++n) { - // free (prsptrs[n]); - // free (dirptrs[n]); - // } - - // free (varptrs); - // free (resptrs); - // free (prsptrs); - // free (dirptrs); - - } END_GLOBAL_MODE; +// for (n=0; n<nvars; ++n) { +// free (prsptrs[n]); +// free (dirptrs[n]); +// } + +// free (varptrs); +// free (resptrs); +// free (prsptrs); +// free (dirptrs); + + free (common::fromindex); + free (common::toindex); + free (common::nboundaryzones); + + // Restore state + set_reflevel ((cGH *)cctkGH, saved_reflevel); + set_mglevel ((cGH *)cctkGH, saved_mglevel); return 0; } |