Import the recently announced changes:

Import the recently announced changes:
1. Carpet has now an infrastructure for multiple maps (aka "grid
patches").  Instead of a single grid hierarchy there can now be
several.  This is largely untested, because the remainder of Cactus
cannot handle multiple coordinate systems.

2. The order in which the schedule bins are called has changed.  As Ian
Hawke pointed out, the previous order during time evolution was
inconsistent.  The initial data ordering did not allow for recovering
and was not usable for progressively solving elliptic equations for
initial data.

3. Carpet now supports convergence levels.  The convergence level
specifies by how many factors of two the resolution in the parameter
file should be coarsened (or refined, if negative).  This should make
convergence tests and test runs much easier.  It is, in principle, also
possible to run several convergence levels at once.  This has not been
tested because the remainder of Cactus cannot handle multiple
resolutions.  This will be necessary for a multigrid solver, and also
for having a shadow hierarchy to determine where to refine adaptively.

4. Carpet works together with the new CoordBase domain specification
parameters.  Without these, using convergence levels will lead to very
strange results.

5. The "modes" have changed.  There are now:
        meta mode:      the whole simulation
        global mode:    one convergence level
        level mode:     one refinement level
        singlemap mode: one map on one refinement level
        local mode:     as previously
The whole mode handling has been cleaned up.

6. The regridding thorn has been cleaned up.

7. The kind of prolongation stencil is now determined in Carpet, i.e. at
a fairly hight level, instead of in CarpetLib.

8. The low-order prolongation operators have been made much more
efficient (as have previously the higher-order ones).

9. Assorted smaller changes.

For Carpet users, there should be no major incompatibilities.  The major
improvements are 3 and 4 combined.  Here is an example:

CoordBase::domainsize    = extent
CoordBase::spacing       = gridspacing
CoordBase::zero_origin_x = yes
CoordBase::zero_origin_y = yes
CoordBase::zero_origin_z = yes
CoordBase::xextent       = 20.0
CoordBase::yextent       = 20.0
CoordBase::zextent       = 20.0
CoordBase::dx            = 1.0
CoordBase::dy            = 1.0
CoordBase::dz            = 1.0

CoordBase::boundary_shiftout_x_lower = 1
CoordBase::boundary_shiftout_y_lower = 1
CoordBase::boundary_shiftout_z_lower = 1

Carpet::domain_from_coordbase   = yes
Carpet::convergence_level       = 0

grid::type         = coordbase
grid::domain       = octant
grid::avoid_origin = no

This gives you a grid that extends from the origin ("zero_origin") up to
20.0 with a grid spacing of 1.0.  Symmetry zones and boundary zones are
added automatically.  The "shiftout" says that there is no boundary
point on the origin.  The staggering parameters (not shown) default to
"no".  In order to change the resolution, only the convergence level
has to be adjusted.  Note that the old way of specifying the domain
extent still works.

For Carpet developers, one major change is the new mode handling.  As
described in 5, the looping macros (that loop over all refinement
levels, or all components) have changed.

darcs-hash:20040125125727-07bb3-3368611314b2dcb8c8ae58ab3f501b683d7edb8f.gz

1 files changed, 429 insertions, 443 deletions

diff --git a/CarpetDev/CarpetCG/src/CG.cc b/CarpetDev/CarpetCG/src/CG.cc
index 37d45e8d6..e27e465be 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.3 2003/11/20 08:34:03 hawke Exp $ */
+/* $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/CG.cc,v 1.4 2004/01/25 14:57:28 schnetter Exp $ */
 
 #include <cassert>
 #include <cmath>
@@ -680,345 +680,280 @@ namespace CarpetCG {
     assert (minerror >= 0);
     assert (sminerror >= 0);
     assert (stepsize > 0);
-        
-    // 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;
-    if (reflevel!=-1) {
-      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
-     */
+
+    // Switch to global mode
+    BEGIN_GLOBAL_MODE(cctkGH) {
       
-    /*
-     * 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)
-     */
+      // 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;
       
-    /* 
-     * 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
-     */
+      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)
+       */
+      
+      /* 
+       * 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 
-    int varptrs[nvars];
-    int resptrs[nvars];
-    //    CCTK_REAL * restrict * restrict prsptrs;
-    //    CCTK_REAL * restrict * restrict dirptrs;  
+      // 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;  
       
-    // 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);
-    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()));
-    }
+      // 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()));
+      }
   
   
   
-    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);
-
-    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;
+      /* 03. r <= - f'(x) */
+      /* 04. calculate M \approx f''(x) */
+      common::ierr = 0;
+      CallLocalFunction((cGH *)cctkGH, common::call_calcres);
+      ierr = common::ierr;
+      assert (! ierr);
 
-      /* 11. delta_d <= d^T d */
       for (int n = 0; n < nvars; n++) {
-        common::fromindex[n] = dirptrs[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_dot_product);
-      delta_d = common::realoutput;
-
-//       cout << "delta_d " << delta_d << endl;
-    
-      /* 12. alpha <= - sigma_0 */
-      alpha = - stepsize;
-      sum_alpha = alpha;
-      do_abort = 0;
+      CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign);
     
-      /* 13. eta_prev <= [f'(x + sigma_0 d)]^T d */
+      /* 05. s <= M^-1 r */
+      // No preconditioning
       for (int n = 0; n < nvars; n++) {
-        common::fromindex[n] = dirptrs[n];
-        common::toindex[n] = varptrs[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;
       }
-      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);
     
-      ierr = 0;
-      CallLocalFunction((cGH *)cctkGH, common::call_calcres);
-      assert (! ierr);
-
+      //    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] = resptrs[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_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] = dirptrs[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);
-      eta = - common::realoutput;
-
-//       cout << "eta " << eta << endl;
-
+      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] = dirptrs[n];
-        common::toindex[n] = varptrs[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;
       }
-      common::realconstant = alpha;
-//       cout << "Add scaled, const " << alpha << endl;
-
-      CallLocalFunction((cGH *)cctkGH, common::call_add_scaled);
-      ierr = 0;
-      CallLocalFunction((cGH *)cctkGH, common::call_applybnds);
-      assert (! ierr);
     
-      eta_prev = eta;
+      CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+      epsilon = common::realoutput;
+      //     cout << "epsilon " << epsilon << endl;
     
-      /* 14. DO */
-      do {
+      /* 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 (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 (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;
+            }
+          }
         }
-      
-        /* 15. eta <= [f'(x)]^T d */
+    
+        /* 10. j <= 0 */
+        siter = 0;
+
+        /* 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;
+    
+        /* 12. alpha <= - sigma_0 */
+        alpha = - stepsize;
+        sum_alpha = alpha;
+        do_abort = 0;
+    
+        /* 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);
+    
         ierr = 0;
         CallLocalFunction((cGH *)cctkGH, common::call_calcres);
         assert (! ierr);
@@ -1030,7 +965,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];
@@ -1039,34 +974,10 @@ namespace CarpetCG {
           common::fromindex[n] = -1;
           common::toindex[n] = -1;
         }
-
         CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
         eta = - common::realoutput;
 
-//         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 */
+        //       cout << "eta " << eta << endl;
 
         for (int n = 0; n < nvars; n++) {
           common::fromindex[n] = dirptrs[n];
@@ -1077,174 +988,249 @@ 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);
     
-      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)));
-      }
+        /* 14. DO */
+        do {
     
-      /* 21. r <= - f'(x) */
-      /* 24. calculate M \approx f''(x) */
-      ierr = 0;
-      CallLocalFunction((cGH *)cctkGH, common::call_calcres);
-      assert (! ierr);
+          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);
-    
-      /* 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;
-      }
+          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);
 
-      CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
-      delta_mid = common::realoutput;
-//       cout << "delta_mid " << delta_mid << endl;
-      
-      /* 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;
-      }
-      
-      CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner);
+          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;
+          }
 
-      /* 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_dot_product);
+          eta = - common::realoutput;
 
-      CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
-      delta_new = common::realoutput;
-//       cout << "delta_new " << delta_new << 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] = resptrs[n];
-        common::toindex[n] = resptrs[n];
-      }
-      for (int n = nvars; n < cg_maxsolvevars; n++) {
-        common::fromindex[n] = -1;
-        common::toindex[n] = -1;
-      }
+          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_dot_product);
-      epsilon = common::realoutput;
-//       cout << "epsilon " << epsilon << 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);
     
-      /* 27. beta <= (delta_new - delta_mid) / (delta_old) */
-      beta = (delta_new - delta_mid) / delta_old;
+        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)));
+        }
     
-      /* 28. k <= k + 1 */
-      ++ iter2;
+        /* 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);
     
-      /* 29. IF k = n OR beta <= 0 THEN */
-      if (iter2 >= maxiters2 || beta <= 0 || do_abort) {
-        /* Restart */
-      
-        /* 30. d <= s */
+        /* 23. delta_mid <= r^T s */
         for (int n = 0; n < nvars; n++) {
           common::fromindex[n] = prsptrs[n];
-          common::toindex[n] = dirptrs[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_copy);
-      
-        /* 31. k <= 0 */
-        iter2 = 0;
+
+        CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+        delta_mid = common::realoutput;
+        //       cout << "delta_mid " << delta_mid << endl;
       
-        /* 32. ELSE */
-      } else {
-        /* Continue with CG iterations */
+        /* 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;
+        }
       
-        /* 33. d <= s + beta d */
+        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);
         for (int n = 0; n < nvars; n++) {
           common::fromindex[n] = prsptrs[n];
-          common::toindex[n] = dirptrs[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];
         }
         for (int n = nvars; n < cg_maxsolvevars; n++) {
           common::fromindex[n] = -1;
           common::toindex[n] = -1;
         }
-        common::realconstant = beta;
-        //        cout << "beta " << beta << endl;
+
+        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;
         
-        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);
-
-    // Restore state
-    if (reflevel!=saved_reflevel) {
-      set_reflevel ((cGH *)cctkGH, saved_reflevel);
-      set_mglevel ((cGH *)cctkGH, saved_mglevel);
-    }
+      //     for (n=0; n<nvars; ++n) {
+      //       free (prsptrs[n]);
+      //       free (dirptrs[n]);
+      //     }
+      
+      //     free (varptrs);
+      //     free (resptrs);
+      //     free (prsptrs);
+      //     free (dirptrs);
+      
+    } END_GLOBAL_MODE;
     
     return 0;
   }