Added ability to pass callback to MADWF that is called every inner iteration and allows user to, for example, adjust the inner solver tolerance depending on residual

Added a general implementation of the Remez algorithm for producing arbitrary rational polynomial approximation with optional restriction to even/odd polynomials
Added implementation of computation of ZMobius parameters
Added Test_zMADWF_prec to test ZMobius in MADWF

Grid/algorithms/Algorithms.h

@@ -38,6 +38,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/approx/Remez.h>
 #include <Grid/algorithms/approx/MultiShiftFunction.h>
 #include <Grid/algorithms/approx/Forecast.h>
+#include <Grid/algorithms/approx/RemezGeneral.h>
+#include <Grid/algorithms/approx/ZMobius.h>
 
 #include <Grid/algorithms/iterative/Deflation.h>
 #include <Grid/algorithms/iterative/ConjugateGradient.h>

Grid/algorithms/approx/RemezGeneral.cc

@@ -0,0 +1,473 @@
+#include<math.h>
+#include<stdio.h>
+#include<stdlib.h>
+#include<string>
+#include<iostream>
+#include<iomanip>
+#include<cassert>
+
+#include<Grid/algorithms/approx/RemezGeneral.h>
+
+
+// Constructor
+AlgRemezGeneral::AlgRemezGeneral(double lower, double upper, long precision,
+				 bigfloat (*f)(bigfloat x, void *data), void *data): f(f), 
+										     data(data), 
+										     prec(precision),
+										     apstrt(lower), apend(upper), apwidt(upper - lower),
+										     n(0), d(0), pow_n(0), pow_d(0)
+{
+  bigfloat::setDefaultPrecision(prec);
+
+  std::cout<<"Approximation bounds are ["<<apstrt<<","<<apend<<"]\n";
+  std::cout<<"Precision of arithmetic is "<<precision<<std::endl;
+}
+
+//Determine the properties of the numerator and denominator polynomials
+void AlgRemezGeneral::setupPolyProperties(int num_degree, int den_degree, PolyType num_type_in, PolyType den_type_in){
+  pow_n = num_degree;
+  pow_d = den_degree;
+
+  if(pow_n % 2 == 0 && num_type_in == PolyType::Odd) assert(0);
+  if(pow_n % 2 == 1 && num_type_in == PolyType::Even) assert(0);
+
+  if(pow_d % 2 == 0 && den_type_in == PolyType::Odd) assert(0);
+  if(pow_d % 2 == 1 && den_type_in == PolyType::Even) assert(0);
+
+  num_type = num_type_in;
+  den_type = den_type_in;
+
+  num_pows.resize(pow_n+1);
+  den_pows.resize(pow_d+1);
+
+  int n_in = 0;
+  bool odd = num_type == PolyType::Full || num_type == PolyType::Odd;
+  bool even = num_type == PolyType::Full || num_type == PolyType::Even;
+  for(int i=0;i<=pow_n;i++){
+    num_pows[i] = -1;
+    if(i % 2 == 0 && even) num_pows[i] = n_in++;
+    if(i % 2 == 1 && odd) num_pows[i] = n_in++;
+  }
+
+  std::cout << n_in << " terms in numerator" << std::endl;
+  --n_in; //power is 1 less than the number of terms, eg  pow=1   a x^1  + b x^0
+
+  int d_in = 0;
+  odd = den_type == PolyType::Full || den_type == PolyType::Odd;
+  even = den_type == PolyType::Full || den_type == PolyType::Even;
+  for(int i=0;i<=pow_d;i++){
+    den_pows[i] = -1;
+    if(i % 2 == 0 && even) den_pows[i] = d_in++;
+    if(i % 2 == 1 && odd) den_pows[i] = d_in++;
+  }
+
+  std::cout << d_in << " terms in denominator" << std::endl;
+  --d_in;
+
+  n = n_in;
+  d = d_in;
+}
+
+//Setup algorithm
+void AlgRemezGeneral::reinitializeAlgorithm(){
+  spread = 1.0e37;
+  iter = 0;
+
+  neq = n + d + 1; //not +2 because highest-power term in denominator is fixed to 1
+
+  param.resize(neq);
+  yy.resize(neq+1);
+
+  //Initialize linear equation temporaries
+  A.resize(neq*neq);
+  B.resize(neq);
+  IPS.resize(neq);
+
+  //Initialize maximum and minimum errors
+  xx.resize(neq+2);
+  mm.resize(neq+1);
+  initialGuess();
+
+  //Initialize search steps
+  step.resize(neq+1);
+  stpini();
+}
+
+double AlgRemezGeneral::generateApprox(const int num_degree, const int den_degree, 
+				       const PolyType num_type_in, const PolyType den_type_in, 
+				       const double _tolerance, const int report_freq){
+  //Setup the properties of the polynomial
+  setupPolyProperties(num_degree, den_degree, num_type_in, den_type_in);
+
+  //Setup the algorithm
+  reinitializeAlgorithm();
+
+  bigfloat tolerance = _tolerance;
+
+  //Iterate until convergance
+  while (spread > tolerance) { 
+    if (iter++ % report_freq==0)
+      std::cout<<"Iteration " <<iter-1<<" spread "<<(double)spread<<" delta "<<(double)delta << std::endl; 
+
+    equations();
+    if (delta < tolerance) {
+      std::cout<<"Iteration " << iter-1 << " delta too small (" << delta << "<" << tolerance << "), try increasing precision\n";
+      assert(0);
+    };    
+    assert( delta>= tolerance );
+
+    search();
+  }
+
+  int sign;
+  double error = (double)getErr(mm[0],&sign);
+  std::cout<<"Converged at "<<iter<<" iterations; error = "<<error<<std::endl;
+
+  // Return the maximum error in the approximation
+  return error;
+}
+
+
+// Initial values of maximal and minimal errors
+void AlgRemezGeneral::initialGuess(){
+  // Supply initial guesses for solution points
+  long ncheb = neq;			// Degree of Chebyshev error estimate
+
+  // Find ncheb+1 extrema of Chebyshev polynomial
+  bigfloat a = ncheb;
+  bigfloat r;
+
+  mm[0] = apstrt;
+  for (long i = 1; i < ncheb; i++) {
+    r = 0.5 * (1 - cos((M_PI * i)/(double) a));
+    //r *= sqrt_bf(r);
+    r = (exp((double)r)-1.0)/(exp(1.0)-1.0);
+    mm[i] = apstrt + r * apwidt;
+  }
+  mm[ncheb] = apend;
+
+  a = 2.0 * ncheb;
+  for (long i = 0; i <= ncheb; i++) {
+    r = 0.5 * (1 - cos(M_PI * (2*i+1)/(double) a));
+    //r *= sqrt_bf(r); // Squeeze to low end of interval
+    r = (exp((double)r)-1.0)/(exp(1.0)-1.0);
+    xx[i] = apstrt + r * apwidt;
+  }
+}
+
+// Initialise step sizes
+void AlgRemezGeneral::stpini(){
+  xx[neq+1] = apend;
+  delta = 0.25;
+  step[0] = xx[0] - apstrt;
+  for (int i = 1; i < neq; i++) step[i] = xx[i] - xx[i-1];
+  step[neq] = step[neq-1];
+}
+
+// Search for error maxima and minima
+void AlgRemezGeneral::search(){
+  bigfloat a, q, xm, ym, xn, yn, xx1;
+  int emsign, ensign, steps;
+
+  int meq = neq + 1;
+
+  bigfloat eclose = 1.0e30;
+  bigfloat farther = 0l;
+
+  bigfloat xx0 = apstrt;
+
+  for (int i = 0; i < meq; i++) {
+    steps = 0;
+    xx1 = xx[i]; // Next zero
+    if (i == meq-1) xx1 = apend;
+    xm = mm[i];
+    ym = getErr(xm,&emsign);
+    q = step[i];
+    xn = xm + q;
+    if (xn < xx0 || xn >= xx1) {	// Cannot skip over adjacent boundaries
+      q = -q;
+      xn = xm;
+      yn = ym;
+      ensign = emsign;
+    } else {
+      yn = getErr(xn,&ensign);
+      if (yn < ym) {
+	q = -q;
+	xn = xm;
+	yn = ym;
+	ensign = emsign;
+      }
+    }
+  
+    while(yn >= ym) {		// March until error becomes smaller.
+      if (++steps > 10)
+      	break;
+      
+      ym = yn;
+      xm = xn;
+      emsign = ensign;
+      a = xm + q;
+      if (a == xm || a <= xx0 || a >= xx1)
+	break;// Must not skip over the zeros either side.      
+
+      xn = a;
+      yn = getErr(xn,&ensign);
+    }
+
+    mm[i] = xm;			// Position of maximum
+    yy[i] = ym;			// Value of maximum
+
+    if (eclose > ym) eclose = ym;
+    if (farther < ym) farther = ym;
+
+    xx0 = xx1; // Walk to next zero.
+  } // end of search loop
+
+  q = (farther - eclose);	// Decrease step size if error spread increased
+
+  if (eclose != 0.0) q /= eclose; // Relative error spread
+
+  if (q >= spread)
+    delta *= 0.5; // Spread is increasing; decrease step size
+  
+  spread = q;
+
+  for (int i = 0; i < neq; i++) {
+    q = yy[i+1];
+    if (q != 0.0) q = yy[i] / q  - (bigfloat)1l;
+    else q = 0.0625;
+    if (q > (bigfloat)0.25) q = 0.25;
+    q *= mm[i+1] - mm[i];
+    step[i] = q * delta;
+  }
+  step[neq] = step[neq-1];
+  
+  for (int i = 0; i < neq; i++) {	// Insert new locations for the zeros.
+    xm = xx[i] - step[i];
+
+    if (xm <= apstrt)
+      continue;
+
+    if (xm >= apend)
+      continue;
+
+    if (xm <= mm[i])
+      xm = (bigfloat)0.5 * (mm[i] + xx[i]);    
+
+    if (xm >= mm[i+1])
+      xm = (bigfloat)0.5 * (mm[i+1] + xx[i]);
+    
+    xx[i] = xm;
+  }
+}
+
+// Solve the equations
+void AlgRemezGeneral::equations(){
+  bigfloat x, y, z;
+  bigfloat *aa;
+  
+  for (int i = 0; i < neq; i++) {	// set up the equations for solution by simq()
+    int ip = neq * i;		// offset to 1st element of this row of matrix
+    x = xx[i];			// the guess for this row
+    y = func(x);		// right-hand-side vector
+
+    z = (bigfloat)1l;
+    aa = A.data()+ip;
+    int t = 0;
+    for (int j = 0; j <= pow_n; j++) {
+      if(num_pows[j] != -1){ *aa++ = z; t++; }
+      z *= x;
+    }
+    assert(t == n+1);
+
+    z = (bigfloat)1l;
+    t = 0;
+    for (int j = 0; j < pow_d; j++) {
+      if(den_pows[j] != -1){ *aa++ = -y * z; t++; }
+      z *= x;
+    }
+    assert(t == d);
+
+    B[i] = y * z;		// Right hand side vector
+  }
+
+  // Solve the simultaneous linear equations.
+  if (simq()){
+    std::cout<<"simq failed\n";
+    exit(0);
+  }
+}
+
+
+// Evaluate the rational form P(x)/Q(x) using coefficients
+// from the solution vector param
+bigfloat AlgRemezGeneral::approx(const bigfloat x) const{
+  // Work backwards toward the constant term.
+  int c = n;
+  bigfloat yn = param[c--];		// Highest order numerator coefficient
+  for (int i = pow_n-1; i >= 0; i--) yn = x * yn  +  (num_pows[i] != -1 ? param[c--] : bigfloat(0l));  
+
+  c = n+d;
+  bigfloat yd = 1l; //Highest degree coefficient is 1.0
+  for (int i = pow_d-1; i >= 0; i--) yd = x * yd  +  (den_pows[i] != -1 ? param[c--] : bigfloat(0l)); 
+
+  return(yn/yd);
+}
+
+// Compute size and sign of the approximation error at x
+bigfloat AlgRemezGeneral::getErr(bigfloat x, int *sign) const{
+  bigfloat f = func(x);
+  bigfloat e = approx(x) - f;
+  if (f != 0) e /= f;
+  if (e < (bigfloat)0.0) {
+    *sign = -1;
+    e = -e;
+  }
+  else *sign = 1;
+  
+  return(e);
+}
+
+// Solve the system AX=B
+int AlgRemezGeneral::simq(){
+
+  int ip, ipj, ipk, ipn;
+  int idxpiv;
+  int kp, kp1, kpk, kpn;
+  int nip, nkp;
+  bigfloat em, q, rownrm, big, size, pivot, sum;
+  bigfloat *aa;
+  bigfloat *X = param.data();
+
+  int n = neq;
+  int nm1 = n - 1;
+  // Initialize IPS and X
+  
+  int ij = 0;
+  for (int i = 0; i < n; i++) {
+    IPS[i] = i;
+    rownrm = 0.0;
+    for(int j = 0; j < n; j++) {
+      q = abs_bf(A[ij]);
+      if(rownrm < q) rownrm = q;
+      ++ij;
+    }
+    if (rownrm == (bigfloat)0l) {
+      std::cout<<"simq rownrm=0\n";
+      return(1);
+    }
+    X[i] = (bigfloat)1.0 / rownrm;
+  }
+  
+  for (int k = 0; k < nm1; k++) {
+    big = 0.0;
+    idxpiv = 0;
+    
+    for (int i = k; i < n; i++) {
+      ip = IPS[i];
+      ipk = n*ip + k;
+      size = abs_bf(A[ipk]) * X[ip];
+      if (size > big) {
+	big = size;
+	idxpiv = i;
+      }
+    }
+    
+    if (big == (bigfloat)0l) {
+      std::cout<<"simq big=0\n";
+      return(2);
+    }
+    if (idxpiv != k) {
+      int j = IPS[k];
+      IPS[k] = IPS[idxpiv];
+      IPS[idxpiv] = j;
+    }
+    kp = IPS[k];
+    kpk = n*kp + k;
+    pivot = A[kpk];
+    kp1 = k+1;
+    for (int i = kp1; i < n; i++) {
+      ip = IPS[i];
+      ipk = n*ip + k;
+      em = -A[ipk] / pivot;
+      A[ipk] = -em;
+      nip = n*ip;
+      nkp = n*kp;
+      aa = A.data()+nkp+kp1;
+      for (int j = kp1; j < n; j++) {
+	ipj = nip + j;
+	A[ipj] = A[ipj] + em * *aa++;
+      }
+    }
+  }
+  kpn = n * IPS[n-1] + n - 1;	// last element of IPS[n] th row
+  if (A[kpn] == (bigfloat)0l) {
+    std::cout<<"simq A[kpn]=0\n";
+    return(3);
+  }
+
+  
+  ip = IPS[0];
+  X[0] = B[ip];
+  for (int i = 1; i < n; i++) {
+    ip = IPS[i];
+    ipj = n * ip;
+    sum = 0.0;
+    for (int j = 0; j < i; j++) {
+      sum += A[ipj] * X[j];
+      ++ipj;
+    }
+    X[i] = B[ip] - sum;
+  }
+  
+  ipn = n * IPS[n-1] + n - 1;
+  X[n-1] = X[n-1] / A[ipn];
+  
+  for (int iback = 1; iback < n; iback++) {
+    //i goes (n-1),...,1
+    int i = nm1 - iback;
+    ip = IPS[i];
+    nip = n*ip;
+    sum = 0.0;
+    aa = A.data()+nip+i+1;
+    for (int j= i + 1; j < n; j++) 
+      sum += *aa++ * X[j];
+    X[i] = (X[i] - sum) / A[nip+i];
+  }
+  
+  return(0);
+}
+
+void AlgRemezGeneral::csv(std::ostream & os) const{
+  os << "Numerator" << std::endl;
+  for(int i=0;i<=pow_n;i++){
+    os << getCoeffNum(i) << "*x^" << i;
+    if(i!=pow_n) os << " + ";
+  }
+  os << std::endl;
+
+  os << "Denominator" << std::endl;
+  for(int i=0;i<=pow_d;i++){
+    os << getCoeffDen(i) << "*x^" << i;
+    if(i!=pow_d) os << " + ";
+  }
+  os << std::endl;
+
+  //For a true minimax solution the errors should all be equal and the signs should oscillate +-+-+- etc
+  int sign;
+  os << "Errors at maxima: coordinate, error, (sign)" << std::endl;
+  for(int i=0;i<neq+1;i++){ 
+    os << mm[i] << " " << getErr(mm[i],&sign) << " (" << sign << ")" << std::endl;
+  }
+
+  os << "Scan over range:" << std::endl;
+  int npt = 60;
+  bigfloat dlt = (apend - apstrt)/bigfloat(npt-1);
+
+  for (bigfloat x=apstrt; x<=apend; x = x + dlt) {
+    double f = evaluateFunc(x);
+    double r = evaluateApprox(x);
+    os<< x<<","<<r<<","<<f<<","<<r-f<<std::endl;
+  }
+  return;
+}

Grid/algorithms/approx/RemezGeneral.h

@@ -0,0 +1,170 @@
+/*
+  C.Kelly Jan 2020 based on implementation by M. Clark May 2005
+
+  AlgRemezGeneral is an implementation of the Remez algorithm for approximating an arbitrary function by a rational polynomial 
+  It includes optional restriction to odd/even polynomials for the numerator and/or denominator
+*/
+
+#ifndef INCLUDED_ALG_REMEZ_GENERAL_H
+#define INCLUDED_ALG_REMEZ_GENERAL_H
+
+#include <stddef.h>
+#include <Grid/GridStd.h>
+
+#ifdef HAVE_LIBGMP
+#include "bigfloat.h"
+#else
+#include "bigfloat_double.h"
+#endif
+
+
+class AlgRemezGeneral{
+ public:
+  enum PolyType { Even, Odd, Full };
+
+ private:
+
+  // In GSL-style, pass the function as a function pointer. Any data required to evaluate the function is passed in as a void pointer
+  bigfloat (*f)(bigfloat x, void *data);
+  void *data;
+
+  // The approximation parameters
+  std::vector<bigfloat> param;
+  bigfloat norm;
+
+  // The number of non-zero terms in the numerator and denominator
+  int n, d;
+  // The numerator and denominator degree (i.e.  the largest power)
+  int pow_n, pow_d;
+  
+  // Specify if the numerator and/or denominator are odd/even polynomials
+  PolyType num_type;
+  PolyType den_type;
+  std::vector<int> num_pows; //contains the mapping, with -1 if not present
+  std::vector<int> den_pows;
+
+  // The bounds of the approximation
+  bigfloat apstrt, apwidt, apend;
+
+  // Variables used to calculate the approximation
+  int nd1, iter;
+  std::vector<bigfloat> xx;
+  std::vector<bigfloat> mm;
+  std::vector<bigfloat> step;
+
+  bigfloat delta, spread;
+  
+  // Variables used in search
+  std::vector<bigfloat> yy;
+
+  // Variables used in solving linear equations
+  std::vector<bigfloat> A;
+  std::vector<bigfloat> B;
+  std::vector<int> IPS;
+
+  // The number of equations we must solve at each iteration (n+d+1)
+  int neq;
+
+  // The precision of the GNU MP library
+  long prec;
+
+  // Initialize member variables associated with the polynomial's properties
+  void setupPolyProperties(int num_degree, int den_degree, PolyType num_type_in, PolyType den_type_in);
+
+  // Initial values of maximal and minmal errors
+  void initialGuess();
+
+  // Initialise step sizes
+  void stpini();
+
+  // Initialize the algorithm
+  void reinitializeAlgorithm();
+
+  // Solve the equations
+  void equations();
+
+  // Search for error maxima and minima
+  void search(); 
+
+  // Calculate function required for the approximation
+  inline bigfloat func(bigfloat x) const{
+    return f(x, data);
+  }
+
+  // Compute size and sign of the approximation error at x
+  bigfloat getErr(bigfloat x, int *sign) const;
+
+  // Solve the system AX=B   where X = param
+  int simq();
+
+  // Evaluate the rational form P(x)/Q(x) using coefficients from the solution vector param
+  bigfloat approx(bigfloat x) const;
+
+ public:
+  
+  AlgRemezGeneral(double lower, double upper, long prec,
+		  bigfloat (*f)(bigfloat x, void *data), void *data);
+
+  inline int getDegree(void) const{ 
+    assert(n==d);
+    return n;
+  }
+  // Reset the bounds of the approximation
+  inline void setBounds(double lower, double upper) {
+    apstrt = lower;
+    apend = upper;
+    apwidt = apend - apstrt;
+  }
+
+  // Get the bounds of the approximation
+  inline void getBounds(double &lower, double &upper) const{ 
+    lower=(double)apstrt;
+    upper=(double)apend;
+  }
+
+  // Run the algorithm to generate the rational approximation
+  double generateApprox(int num_degree, int den_degree, 
+			PolyType num_type, PolyType den_type,
+			const double tolerance = 1e-15, const int report_freq = 1000);
+  
+  inline double generateApprox(int num_degree, int den_degree, 
+			       const double tolerance = 1e-15, const int report_freq = 1000){
+    return generateApprox(num_degree, den_degree, Full, Full, tolerance, report_freq);
+  }
+  
+  // Evaluate the rational form P(x)/Q(x) using coefficients from the
+  // solution vector param
+  inline double evaluateApprox(double x) const{
+    return (double)approx((bigfloat)x);
+  }
+
+  // Evaluate the rational form Q(x)/P(x) using coefficients from the solution vector param
+  inline double evaluateInverseApprox(double x) const{
+    return 1.0/(double)approx((bigfloat)x);
+  }  
+
+  // Calculate function required for the approximation
+  inline double evaluateFunc(double x) const{
+    return (double)func((bigfloat)x);
+  }
+
+  // Calculate inverse function required for the approximation
+  inline double evaluateInverseFunc(double x) const{
+    return 1.0/(double)func((bigfloat)x);
+  }
+
+  // Dump csv of function, approx and error
+  void csv(std::ostream &os = std::cout) const;
+
+  // Get the coefficient of the term x^i in the numerator
+  inline double getCoeffNum(const int i) const{    
+    return num_pows[i] == -1 ? 0. : double(param[num_pows[i]]);
+  }
+  // Get the coefficient of the term x^i in the denominator
+  inline double getCoeffDen(const int i) const{ 
+    if(i == pow_d) return 1.0;
+    else return den_pows[i] == -1 ? 0. : double(param[den_pows[i]+n+1]); 
+  }
+};
+
+#endif

Grid/algorithms/approx/ZMobius.cc

@@ -0,0 +1,183 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/approx/ZMobius.cc
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@phys.columbia.edu>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/algorithms/approx/ZMobius.h>
+#include <Grid/algorithms/approx/RemezGeneral.h>
+
+NAMESPACE_BEGIN(Grid);
+NAMESPACE_BEGIN(Approx);
+
+//Compute the tanh approximation
+inline double epsilonMobius(const double x, const std::vector<ComplexD> &w){
+  int Ls = w.size();
+
+  ComplexD fxp = 1., fmp = 1.;
+  for(int i=0;i<Ls;i++){
+    fxp = fxp * ( w[i] + x );
+    fmp = fmp * ( w[i] - x );
+  }
+  return ((fxp - fmp)/(fxp + fmp)).real();
+}
+inline double epsilonMobius(const double x, const std::vector<RealD> &w){
+  int Ls = w.size();
+
+  double fxp = 1., fmp = 1.;
+  for(int i=0;i<Ls;i++){
+    fxp = fxp * ( w[i] + x );
+    fmp = fmp * ( w[i] - x );
+  }
+  return (fxp - fmp)/(fxp + fmp);
+}
+
+
+
+//Compute the tanh approximation in a form suitable for the Remez
+bigfloat epsilonMobius(bigfloat x, void* data){
+  const std::vector<RealD> &omega = *( (std::vector<RealD> const*)data );
+  bigfloat fxp(1.0);
+  bigfloat fmp(1.0);
+
+  for(int i=0;i<omega.size();i++){
+    fxp = fxp * ( bigfloat(omega[i]) + x);
+    fmp = fmp * ( bigfloat(omega[i]) - x);
+  }
+  return (fxp - fmp)/(fxp + fmp);
+}
+
+//Compute the Zmobius Omega parameters suitable for eigenvalue range   -lambda_bound <= lambda <= lambda_bound
+//Note omega_i = 1/(b_i + c_i)   where b_i and c_i are the Mobius parameters
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out,
+			 const std::vector<RealD> &omega_in, const int Ls_in,
+			 const RealD lambda_bound){
+  assert(omega_in.size() == Ls_in);
+  omega_out.resize(Ls_out);
+
+  //Use the Remez algorithm to generate the appropriate rational polynomial
+  //For odd polynomial, to satisfy Haar condition must take either positive or negative half of range (cf https://arxiv.org/pdf/0803.0439.pdf page 6)  
+  AlgRemezGeneral remez(0, lambda_bound, 64, &epsilonMobius, (void*)&omega_in); 
+  remez.generateApprox(Ls_out-1, Ls_out,AlgRemezGeneral::Odd, AlgRemezGeneral::Even, 1e-15, 100);
+  remez.csv(std::cout);
+
+  //The rational approximation has the form  [ f(x) - f(-x) ] / [ f(x) + f(-x) ]  where  f(x) = \Prod_{i=0}^{L_s-1} ( \omega_i + x )
+  //cf https://academiccommons.columbia.edu/doi/10.7916/D8T72HD7  pg 102
+  //omega_i are therefore the negative of the complex roots of f(x)
+
+  //We can find the roots by recognizing that the eigenvalues of a matrix A are the roots of the characteristic polynomial
+  // \rho(\lambda) = det( A - \lambda I )    where I is the unit matrix
+  //The matrix whose characteristic polynomial is an arbitrary monic polynomial a0 + a1 x + a2 x^2 + ... x^n   is the companion matrix 
+  // A = | 0    1   0    0 0 .... 0 |
+  //     | 0    0   1    0 0 .... 0 |
+  //     | :    :   :    : :      : |
+  //     | 0    0   0    0 0      1
+  //     | -a0 -a1 -a2  ...  ... -an|
+
+
+  //Note the Remez defines the largest power to have unit coefficient
+  std::vector<RealD> coeffs(Ls_out+1);
+  for(int i=0;i<Ls_out+1;i+=2) coeffs[i] = coeffs[i] = remez.getCoeffDen(i); //even powers
+  for(int i=1;i<Ls_out+1;i+=2) coeffs[i] = coeffs[i] = remez.getCoeffNum(i); //odd powers
+
+  std::vector<std::complex<RealD> > roots(Ls_out);
+
+  //Form the companion matrix
+  Eigen::MatrixXd compn(Ls_out,Ls_out);
+  for(int i=0;i<Ls_out-1;i++) compn(i,0) = 0.;
+  compn(Ls_out - 1, 0) = -coeffs[0];
+  
+  for(int j=1;j<Ls_out;j++){
+    for(int i=0;i<Ls_out-1;i++) compn(i,j) = i == j-1 ? 1. : 0.;
+    compn(Ls_out - 1, j) = -coeffs[j];
+  }
+
+  //Eigensolve
+  Eigen::EigenSolver<Eigen::MatrixXd> slv(compn, false);
+
+  const auto & ev = slv.eigenvalues();
+  for(int i=0;i<Ls_out;i++)
+    omega_out[i] = -ev(i);
+
+  //Sort ascending (smallest at start of vector!)
+  std::sort(omega_out.begin(), omega_out.end(), 
+	    [&](const ComplexD &a, const ComplexD &b){ return a.real() < b.real() || (a.real() == b.real() && a.imag() < b.imag()); });
+
+  //McGlynn thesis pg 122 suggest improved iteration counts if magnitude of omega diminishes towards the center of the 5th dimension
+  std::vector<ComplexD> omega_tmp = omega_out;
+  int s_low=0, s_high=Ls_out-1, ss=0;
+  for(int s_from = Ls_out-1; s_from >= 0; s_from--){ //loop from largest omega
+    int s_to;
+    if(ss % 2 == 0){
+      s_to = s_low++;
+    }else{
+      s_to = s_high--;
+    }
+    omega_out[s_to] = omega_tmp[s_from];
+    ++ss;
+  }
+  
+  std::cout << "Resulting omega_i:" << std::endl;  
+  for(int i=0;i<Ls_out;i++)
+    std::cout << omega_out[i] << std::endl;
+
+  std::cout << "Test result matches the approximate polynomial found by the Remez" << std::endl;
+  std::cout << "<x> <remez approx> <poly approx> <diff poly approx remez approx> <exact> <diff poly approx exact>\n";
+  
+  int npt = 60;
+  double dlt = lambda_bound/double(npt-1);
+
+  for (int i =0; i<npt; i++){
+    double x = i*dlt;
+    double r = remez.evaluateApprox(x);
+    double p = epsilonMobius(x, omega_out);
+    double e = epsilonMobius(x, omega_in);
+
+    std::cout << x<< " " << r << " " << p <<" " <<r-p << " " << e << " " << e-p << std::endl;
+  }
+
+}
+  
+//mobius_param = b+c   with b-c=1
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound){
+  std::vector<RealD> omega_in(Ls_in, 1./mobius_param);
+  computeZmobiusOmega(omega_out, Ls_out, omega_in, Ls_in, lambda_bound);
+}
+
+//ZMobius class takes  gamma_i = (b+c) omega_i as its input, where b, c are factored out
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, 
+			 const RealD mobius_param_out, const int Ls_out, 
+			 const RealD mobius_param_in, const int Ls_in,
+			 const RealD lambda_bound){
+  computeZmobiusOmega(gamma_out, Ls_out, mobius_param_in, Ls_in, lambda_bound);
+  for(int i=0;i<Ls_out;i++) gamma_out[i] = gamma_out[i] * mobius_param_out;
+}
+//Assumes mobius_param_out == mobius_param_in
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound){
+  computeZmobiusGamma(gamma_out, mobius_param, Ls_out, mobius_param, Ls_in, lambda_bound);
+}
+
+NAMESPACE_END(Approx);
+NAMESPACE_END(Grid);

Grid/algorithms/approx/ZMobius.h

@@ -0,0 +1,57 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/approx/ZMobius.h
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@phys.columbia.edu>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_ZMOBIUS_APPROX_H
+#define GRID_ZMOBIUS_APPROX_H
+
+#include <Grid/GridCore.h>
+
+NAMESPACE_BEGIN(Grid);
+NAMESPACE_BEGIN(Approx);
+
+//Compute the Zmobius Omega parameters suitable for eigenvalue range   -lambda_bound <= lambda <= lambda_bound
+//Note omega_i = 1/(b_i + c_i)   where b_i and c_i are the Mobius parameters
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out,
+			 const std::vector<RealD> &omega_in, const int Ls_in,
+			 const RealD lambda_bound);
+  
+//mobius_param = b+c   with b-c=1
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound);
+
+//ZMobius class takes  gamma_i = (b+c) omega_i as its input, where b, c are factored out
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, 
+			 const RealD mobius_param_out, const int Ls_out, 
+			 const RealD mobius_param_in, const int Ls_in,
+			 const RealD lambda_bound);
+
+//Assumes mobius_param_out == mobius_param_in
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound);
+
+NAMESPACE_END(Approx);
+NAMESPACE_END(Grid);
+
+#endif

Grid/qcd/action/fermion/MADWF.h

@@ -40,6 +40,11 @@ inline void convert(const Fieldi &from,Fieldo &to)
   to=from;
 }
 
+struct MADWFinnerIterCallbackBase{
+  virtual void operator()(const RealD current_resid){}
+  virtual ~MADWFinnerIterCallbackBase(){}
+};
+
 template<class Matrixo,class Matrixi,class PVinverter,class SchurSolver, class Guesser> 
 class MADWF 
 {
@@ -56,23 +61,29 @@ class MADWF
 
   RealD target_resid;
   int   maxiter;
- public:
 
+  //operator() is called on "callback" at the end of every inner iteration. This allows for example the adjustment of the inner
+  //tolerance to speed up subsequent iteration
+  MADWFinnerIterCallbackBase* callback;
+  
+ public:
   MADWF(Matrixo &_Mato,
 	Matrixi &_Mati,
 	PVinverter &_PauliVillarsSolvero,
 	SchurSolver &_SchurSolveri,
 	Guesser & _Guesseri,
 	RealD resid,
-	int _maxiter) :
+	int _maxiter,
+	MADWFinnerIterCallbackBase* _callback = NULL) :
 
   Mato(_Mato),Mati(_Mati),
     SchurSolveri(_SchurSolveri),
-    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri)
+    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri),
-  {   
+    callback(_callback)
-    target_resid=resid;
+    {
-    maxiter     =_maxiter; 
+      target_resid=resid;
-  };
+      maxiter     =_maxiter;
+    };
    
   void operator() (const FermionFieldo &src4,FermionFieldo &sol5)
   {
@@ -177,6 +188,8 @@ class MADWF
        std::cout << GridLogMessage << "Residual " << i << ": " << resid  << std::endl;
        std::cout << GridLogMessage << "***************************************" <<std::endl;
 
+       if(callback != NULL) (*callback)(resid);       
+       
        if (resid < target_resid) {
 	 return;
        }

tests/solver/Test_zMADWF_prec.cc

@@ -0,0 +1,303 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/solver/Test_zMADWF_prec.cc
+
+Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@phys.columbia.edu>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+
+//This test computes the zMobius approximation to the Mobius action and uses it within the MADWF context to accelerate an inversion
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+struct TestParams{
+  bool load_config;
+  std::string config_file;
+
+  double mass;
+
+  std::string outer_precon;
+  std::string inner_precon;
+  
+  int Ls_outer;
+  double b_plus_c_outer;
+  double resid_outer;
+  
+  int Ls_inner;
+  double b_plus_c_inner; //irrelevant for ZMobius
+  double resid_inner;
+  bool zmobius_inner;
+  double lambda_max; //upper bound of H_T eigenvalue range required to generate zMobius approximation
+  
+  TestParams(): load_config(true), config_file("ckpoint_lat.1000"), mass(0.01),
+		Ls_outer(24), b_plus_c_outer(2.0), resid_outer(1e-8),
+		Ls_inner(12), b_plus_c_inner(1.0), resid_inner(1e-8), zmobius_inner(true), lambda_max(1.42), outer_precon("Standard"), inner_precon("Standard")
+  {}
+  
+  void write(const std::string &file) const{
+    XmlWriter wr(file);
+#define DOIT(A) wr.writeDefault(#A, A)
+    DOIT(load_config);
+    DOIT(config_file);
+    DOIT(mass);
+    DOIT(outer_precon);
+    DOIT(inner_precon);
+    DOIT(Ls_outer);
+    DOIT(b_plus_c_outer);
+    DOIT(resid_outer);
+    DOIT(Ls_inner);
+    DOIT(b_plus_c_inner);
+    DOIT(resid_inner);
+    DOIT(zmobius_inner);
+    DOIT(lambda_max);
+#undef DOIT
+  }
+  void read(const std::string &file){
+    XmlReader rd(file);
+#define DOIT(A) rd.readDefault(#A, A)
+    DOIT(load_config);
+    DOIT(config_file);
+    DOIT(mass);
+    DOIT(outer_precon);
+    DOIT(inner_precon);
+    DOIT(Ls_outer);
+    DOIT(b_plus_c_outer);
+    DOIT(resid_outer);
+    DOIT(Ls_inner);
+    DOIT(b_plus_c_inner);
+    DOIT(resid_inner);
+    DOIT(zmobius_inner);
+    DOIT(lambda_max);
+#undef DOIT
+  }
+};
+
+struct RunParamsPrecStd{
+  typedef SchurRedBlackDiagMooeeSolve<LatticeFermionD> SchurSolverType;
+
+  template<typename Action>
+  using HermOpType = SchurDiagMooeeOperator<Action, LatticeFermionD>;
+};
+
+struct RunParamsPrecDiagTwo{
+  typedef SchurRedBlackDiagTwoSolve<LatticeFermionD> SchurSolverType;
+
+  template<typename Action>
+  using HermOpType = SchurDiagTwoOperator<Action, LatticeFermionD>;
+};
+
+
+struct CGincreaseTol : public MADWFinnerIterCallbackBase{
+  ConjugateGradient<LatticeFermionD> &cg_inner;  
+  RealD outer_resid;
+
+  CGincreaseTol(ConjugateGradient<LatticeFermionD> &cg_inner,
+		RealD outer_resid): cg_inner(cg_inner), outer_resid(outer_resid){}
+  
+  void operator()(const RealD current_resid){
+    std::cout << "CGincreaseTol with current residual " << current_resid << " changing inner tolerance " << cg_inner.Tolerance << " -> ";
+    while(cg_inner.Tolerance < current_resid) cg_inner.Tolerance *= 2;    
+    //cg_inner.Tolerance = outer_resid/current_resid;
+    std::cout << cg_inner.Tolerance << std::endl;
+  }
+};
+
+template<typename RunParamsOuter, typename RunParamsInner>
+void run(const TestParams &params){
+  RealD bmc = 1.0; //use Shamir kernel
+  std::vector<ComplexD> gamma_inner;
+
+  std::cout << "Compute parameters" << std::endl;
+  if(params.zmobius_inner){
+    Approx::computeZmobiusGamma(gamma_inner, params.b_plus_c_inner, params.Ls_inner, params.b_plus_c_outer, params.Ls_outer, params.lambda_max);
+  }else{
+    Approx::zolotarev_data *zdata = Approx::higham(1.0,params.Ls_inner);
+    gamma_inner.resize(params.Ls_inner);
+    for(int s=0;s<params.Ls_inner;s++) gamma_inner[s] = zdata->gamma[s];
+    Approx::zolotarev_free(zdata);
+  }
+  std::cout << "gamma:\n";
+  for(int s=0;s<params.Ls_inner;s++) std::cout << s << " " << gamma_inner[s] << std::endl;
+
+
+  GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
+      GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()),
+      GridDefaultMpi());
+  GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+
+  GridCartesian* FGrid_outer = SpaceTimeGrid::makeFiveDimGrid(params.Ls_outer, UGrid);
+  GridCartesian* FGrid_inner = SpaceTimeGrid::makeFiveDimGrid(params.Ls_inner, UGrid);
+
+  GridRedBlackCartesian* FrbGrid_outer = SpaceTimeGrid::makeFiveDimRedBlackGrid(params.Ls_outer, UGrid);
+  GridRedBlackCartesian* FrbGrid_inner = SpaceTimeGrid::makeFiveDimRedBlackGrid(params.Ls_inner, UGrid);
+
+
+  std::vector<int> seeds4({1, 2, 3, 4});
+  std::vector<int> seeds5({5, 6, 7, 8});
+
+  GridParallelRNG RNG5_outer(FGrid_outer);
+  RNG5_outer.SeedFixedIntegers(seeds5);
+
+  GridParallelRNG RNG4(UGrid);
+  RNG4.SeedFixedIntegers(seeds4);
+
+  LatticeFermionD src4(UGrid); random(RNG4,src4);
+
+  LatticeFermionD result_outer(FGrid_outer);
+  result_outer = Zero();
+  LatticeGaugeFieldD Umu(UGrid);
+
+  if(params.load_config){
+    FieldMetaData header;
+    NerscIO::readConfiguration(Umu, header, params.config_file);
+
+    for(int i=0;i<Nd;i++){
+      assert(header.dimension[i] == GridDefaultLatt()[i]);
+    }
+  }else{    
+    SU3::HotConfiguration(RNG4, Umu);
+  }
+    
+  std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt()
+            << "   Ls: " << params.Ls_outer << std::endl;
+
+  RealD M5 = 1.8;
+
+  RealD b_outer = (params.b_plus_c_outer + bmc)/2.;
+  RealD c_outer = (params.b_plus_c_outer - bmc)/2.;
+
+  RealD b_inner = (params.b_plus_c_inner + bmc)/2.;
+  RealD c_inner = (params.b_plus_c_inner - bmc)/2.;
+
+  MobiusFermionD D_outer(Umu, *FGrid_outer, *FrbGrid_outer, *UGrid, *UrbGrid, params.mass, M5, b_outer, c_outer);
+  ZMobiusFermionD D_inner(Umu, *FGrid_inner, *FrbGrid_inner, *UGrid, *UrbGrid, params.mass, M5, gamma_inner, b_inner, c_inner);
+
+  LatticeFermionD src_outer(FGrid_outer);
+  D_outer.ImportPhysicalFermionSource(src4,src_outer); //applies D_- 
+
+  //Solve using a regular even-odd preconditioned CG for the Hermitian operator
+  //M y = x
+  //Mprec y'_o = x'_o     where Mprec = Doo - Doe Dee^-1 Deo    and  x'_o = -Doe Dee^-1 x_e + x_o
+  //y_o = y'_o
+
+  //(Mprec^dag Mprec) y'_o = Mprec^dag x'_o 
+  //y'_o = (Mprec^dag Mprec)^-1 Mprec^dag x'_o 
+
+  //We can get Mprec^dag x'_o from x_o  from SchurRedBlackDiagMooeeSolve::RedBlackSource
+  ConjugateGradient<LatticeFermionD> CG_outer(params.resid_outer, 10000);
+  typename RunParamsOuter::SchurSolverType SchurSolver_outer(CG_outer);
+  
+  LatticeFermionD tmp_e_outer(FrbGrid_outer);
+  LatticeFermionD src_o_outer(FrbGrid_outer);
+  SchurSolver_outer.RedBlackSource(D_outer, src_outer, tmp_e_outer, src_o_outer);
+  
+  LatticeFermionD result_o_outer(FrbGrid_outer);
+  result_o_outer = Zero();
+
+  GridStopWatch CGTimer;
+  
+  typename RunParamsOuter::HermOpType<MobiusFermionD> HermOpEO_outer(D_outer);
+
+  CGTimer.Start();
+  CG_outer(HermOpEO_outer, src_o_outer, result_o_outer);
+  CGTimer.Stop();
+
+  std::cout << GridLogMessage << "Total outer CG time : " << CGTimer.Elapsed()
+            << std::endl;
+
+  CGTimer.Reset();
+
+  //Solve for y using MADWF with internal preconditioning
+
+  //typedef PauliVillarsSolverRBprec<LatticeFermionD, typename RunParamsOuter::SchurSolverType> PVtype;
+  //PVtype PV_outer(SchurSolver_outer);
+
+  typedef PauliVillarsSolverFourierAccel<LatticeFermionD, LatticeGaugeFieldD> PVtype;
+  PVtype PV_outer(Umu, CG_outer);
+
+  ConjugateGradient<LatticeFermionD> CG_inner(params.resid_inner, 10000, 0);
+
+  CGincreaseTol update(CG_inner, params.resid_outer);
+
+  typename RunParamsInner::SchurSolverType SchurSolver_inner(CG_inner);
+
+  ZeroGuesser<LatticeFermion> Guess;
+  MADWF<MobiusFermionD, ZMobiusFermionD, PVtype, typename RunParamsInner::SchurSolverType, ZeroGuesser<LatticeFermion> > madwf(D_outer, D_inner, PV_outer, SchurSolver_inner, Guess, params.resid_outer, 100, &update);
+  
+  LatticeFermionD result_MADWF(FGrid_outer);
+  result_MADWF = Zero();
+
+  CGTimer.Start();
+  madwf(src4, result_MADWF);
+  CGTimer.Stop();
+
+  LatticeFermionD result_o_MADWF(FrbGrid_outer);
+  pickCheckerboard(Odd, result_o_MADWF, result_MADWF);
+
+  std::cout << GridLogMessage << "Total MADWF time : " << CGTimer.Elapsed()
+            << std::endl;
+
+  LatticeFermionD diff = result_o_MADWF - result_o_outer;
+  std::cout <<GridLogMessage<< "Odd-parity MADWF result norm " << norm2(result_o_MADWF) 
+	    << " Regular result norm " << norm2(result_o_outer) 
+	    << " Norm of diff " << norm2(diff)<<std::endl;
+
+
+  //std::cout << GridLogMessage << "######## Dhop calls summary" << std::endl;
+  //D_outer.Report();
+}
+
+
+
+int main(int argc, char** argv) {
+  std::cout << "Init" << std::endl;
+  Grid_init(&argc, &argv);
+
+  TestParams params;
+  
+  if( GridCmdOptionExists(argv,argv+argc,"--params") ){
+    std::string pfile = GridCmdOptionPayload(argv,argv+argc,"--params");
+    if(pfile == "TEMPLATE"){
+      params.write("params.templ");
+      return 0;
+    }else{
+      params.read(pfile);
+    }
+  }
+
+  if(params.outer_precon == "Standard" && params.inner_precon == "Standard" ){
+    run<RunParamsPrecStd, RunParamsPrecStd>(params);
+  }else if(params.outer_precon == "DiagTwo" && params.inner_precon == "Standard"){
+    run<RunParamsPrecDiagTwo, RunParamsPrecStd>(params);
+  }else if(params.outer_precon == "Standard" && params.inner_precon == "DiagTwo"){
+    run<RunParamsPrecStd, RunParamsPrecDiagTwo>(params);
+  }else if(params.outer_precon == "DiagTwo" && params.inner_precon == "DiagTwo"){
+    run<RunParamsPrecDiagTwo, RunParamsPrecDiagTwo>(params);
+  }else assert(0);
+
+  Grid_finalize();
+}