From 41f73ec0836fd5ec3093edba4174b466815cf799 Mon Sep 17 00:00:00 2001
From: David Murphy <djm2131@columbia.edu>
Date: Wed, 16 Aug 2017 12:37:38 -0400
Subject: [PATCH] Add ChronoForecast class for forecasting solutions across
 poles in the EOFA heatbath

---
 lib/algorithms/Algorithms.h      |   3 +-
 lib/algorithms/approx/Forecast.h | 152 +++++++++++++++++++++++++++++++
 2 files changed, 154 insertions(+), 1 deletion(-)
 create mode 100644 lib/algorithms/approx/Forecast.h

diff --git a/lib/algorithms/Algorithms.h b/lib/algorithms/Algorithms.h
index 5123c7a1..361ccd9c 100644
--- a/lib/algorithms/Algorithms.h
+++ b/lib/algorithms/Algorithms.h
@@ -1,6 +1,6 @@
     /*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+    Grid physics library, www.github.com/paboyle/Grid
 
     Source file: ./lib/Algorithms.h
 
@@ -37,6 +37,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/approx/Chebyshev.h>
 #include <Grid/algorithms/approx/Remez.h>
 #include <Grid/algorithms/approx/MultiShiftFunction.h>
+#include <Grid/algorithms/approx/Forecast.h>
 
 #include <Grid/algorithms/iterative/ConjugateGradient.h>
 #include <Grid/algorithms/iterative/ConjugateResidual.h>
diff --git a/lib/algorithms/approx/Forecast.h b/lib/algorithms/approx/Forecast.h
new file mode 100644
index 00000000..87eb84a6
--- /dev/null
+++ b/lib/algorithms/approx/Forecast.h
@@ -0,0 +1,152 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/algorithms/approx/Forecast.h
+
+Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: David Murphy <dmurphy@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 INCLUDED_FORECAST_H
+#define INCLUDED_FORECAST_H
+
+namespace Grid {
+
+  // Abstract base class.
+  // Takes a matrix (Mat), a source (phi), and a vector of Fields (chi)
+  // and returns a forecasted solution to the system D*psi = phi (psi).
+  template<class Matrix, class Field>
+  class Forecast
+  {
+    public:
+      virtual Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& chi) = 0;
+  };
+
+  // Implementation of Brower et al.'s chronological inverter (arXiv:hep-lat/9509012),
+  // used to forecast solutions across poles of the EOFA heatbath.
+  //
+  // Modified from CPS (cps_pp/src/util/dirac_op/d_op_base/comsrc/minresext.C)
+  template<class Matrix, class Field>
+  class ChronoForecast : public Forecast<Matrix,Field>
+  {
+    public:
+      Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& prev_solns)
+      {
+        int degree = prev_solns.size();
+        Field chi(phi); // forecasted solution
+
+        // Trivial cases
+        if(degree == 0){ chi = zero; return chi; }
+        else if(degree == 1){ return prev_solns[0]; }
+
+        RealD dot;
+        ComplexD xp;
+        Field r(phi); // residual
+        Field Mv(phi);
+        std::vector<Field> v(prev_solns); // orthonormalized previous solutions
+        std::vector<Field> MdagMv(degree,phi);
+
+        // Array to hold the matrix elements
+        std::vector<std::vector<ComplexD>> G(degree, std::vector<ComplexD>(degree));
+
+        // Solution and source vectors
+        std::vector<ComplexD> a(degree);
+        std::vector<ComplexD> b(degree);
+
+        // Orthonormalize the vector basis
+        for(int i=0; i<degree; i++){
+          v[i] *= 1.0/std::sqrt(norm2(v[i]));
+          for(int j=i+1; j<degree; j++){ v[j] -= innerProduct(v[i],v[j]) * v[i]; }
+        }
+
+        // Perform sparse matrix multiplication and construct rhs
+        for(int i=0; i<degree; i++){
+          b[i] = innerProduct(v[i],phi);
+          Mat.M(v[i],Mv);
+          Mat.Mdag(Mv,MdagMv[i]);
+          G[i][i] = innerProduct(v[i],MdagMv[i]);
+        }
+
+        // Construct the matrix
+        for(int j=0; j<degree; j++){
+        for(int k=j+1; k<degree; k++){
+          G[j][k] = innerProduct(v[j],MdagMv[k]);
+          G[k][j] = std::conj(G[j][k]);
+        }}
+
+        // Gauss-Jordan elimination with partial pivoting
+        for(int i=0; i<degree; i++){
+
+          // Perform partial pivoting
+          int k = i;
+          for(int j=i+1; j<degree; j++){ if(std::abs(G[j][j]) > std::abs(G[k][k])){ k = j; } }
+          if(k != i){
+            xp = b[k];
+            b[k] = b[i];
+            b[i] = xp;
+            for(int j=0; j<degree; j++){
+              xp = G[k][j];
+              G[k][j] = G[i][j];
+              G[i][j] = xp;
+            }
+          }
+
+          // Convert matrix to upper triangular form
+          for(int j=i+1; j<degree; j++){
+            xp = G[j][i]/G[i][i];
+            b[j] -= xp * b[i];
+            for(int k=0; k<degree; k++){ G[j][k] -= xp*G[i][k]; }
+          }
+        }
+
+        // Use Gaussian elimination to solve equations and calculate initial guess
+        chi = zero;
+        r = phi;
+        for(int i=degree-1; i>=0; i--){
+          a[i] = 0.0;
+          for(int j=i+1; j<degree; j++){ a[i] += G[i][j] * a[j]; }
+          a[i] = (b[i]-a[i])/G[i][i];
+          chi += a[i]*v[i];
+          r -= a[i]*MdagMv[i];
+        }
+
+        RealD true_r(0.0);
+        ComplexD tmp;
+        for(int i=0; i<degree; i++){
+          tmp = -b[i];
+          for(int j=0; j<degree; j++){ tmp += G[i][j]*a[j]; }
+          tmp = std::conj(tmp)*tmp;
+          true_r += std::sqrt(tmp.real());
+        }
+
+        RealD error = std::sqrt(norm2(r)/norm2(phi));
+        std::cout << GridLogMessage << "ChronoForecast: |res|/|src| = " << error << std::endl;
+
+        return chi;
+      };
+  };
+
+}
+
+#endif