feature/initial-guesses (#55)

* [LinearSolver][INS]:  Implement NONE initial guess strategy.

This commit begins the port of the initial guess codes to mainline libP.

* [LinearSolver][INS] Implement ZERO initial guess strategy.

* [LinearSolver][INS] Implement the CLASSIC initial guess strategy.

* [LinearSolver][INS] Implement the QR initial guess strategy.

* [Mesh] Make members implementing basic codes for nodes, matrices, etc. static.

These don't have any associated state, and they're sometimes useful
outside of mesh_t, so this shouldn't hurt anything.

* [LinearSolver][INS] Implement the EXTRAP initial guess strategy.

* [LinearSolver][INS] Enable changing the PCG stopping criterion.

Two options:  stop when we've reduced by TOL relative to the initial residual
(ABS/REL-INITRESID, the default) or relative to the 2-norm of the RHS
(ABS/REL-RHS-2NORM).

Only implemented for PCG.

* [LinearSolver] Take at least 1 iteration in PCG.

* [Elliptic] Enforce continuity of the multigrid precon. output.

(Of course, this is only for continuous discretizations.)

This can make a difference for some of the EXTRAP initial guess methods
that are on the verge of instability.

Since this doesn't seem to matter in most cases, it would be better to
enforce the continuity only at the the output of the linear solver, but
there's an encapsulation problem:  the linear solver object does not
know about discretizations.

* [LinearSolver] Initialize linAlg.set() need by the initial guess methods.

* WIP: first prototype of testing for initial guess

* WIP: added initial guess to test.py

* WIP: adding CLASSIC test for initial guess

* [LinearSolver] Handle initial guess settings in the initial guess module.

In particular, get them out of elliptic.  (Doesn't really make sense for
elliptic to own them, since that's not a time-dependent solver.)

Co-authored-by: tcew <tim.warburton@gmail.com>

include/core.hpp

@@ -43,6 +43,8 @@ void parallelSort(int size, int rank, MPI_Comm comm,
 
 void matrixRightSolve(int NrowsA, int NcolsA, double *A, int NrowsB, int NcolsB, double *B, double *C);
 void matrixRightSolve(int NrowsA, int NcolsA, float *A, int NrowsB, int NcolsB, float *B, float *C);
+void matrixUnderdeterminedRightSolveMinNorm(int NrowsA, int NcolsA, dfloat *A, dfloat *b, dfloat *x);
+void matrixUnderdeterminedRightSolveCPQR(int NrowsA, int NcolsA, dfloat *A, dfloat *b, dfloat *x);
 
 void matrixEigenVectors(int N, double *A, double *VR, double *WR, double *WI);
 void matrixEigenVectors(int N, float *A, float *VR, float *WR, float *WI);
@@ -70,4 +72,4 @@ void matrixTranspose(const int M, const int N,
                      const long long int  *A, const int LDA,
                            long long int *AT, const int LDAT);
 
-#endif
+#endif

include/initialGuess.hpp

@@ -0,0 +1,175 @@
+/*
+
+The MIT License (MIT)
+
+Copyright (c) 2017 Tim Warburton, Noel Chalmers, Jesse Chan, Ali Karakus
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+*/
+
+#ifndef INITIALGUESS_HPP
+#define INITIALGUESS_HPP
+
+#include "linearSolver.hpp"
+
+// Abstract base class for different initial guess strategies.
+class initialGuessStrategy_t {
+protected:
+  platform_t& platform;
+  settings_t& settings;
+  MPI_Comm   comm;
+
+  dlong Ntotal;     // Degrees of freedom
+
+public:
+  initialGuessStrategy_t(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+  virtual ~initialGuessStrategy_t();
+
+  virtual void FormInitialGuess(occa::memory& o_x, occa::memory& o_rhs) = 0;
+  virtual void Update(solver_t& solver, occa::memory& o_x, occa::memory& o_rhs) = 0;
+};
+
+// Default initial guess strategy:  use whatever the user gave us.
+class igDefaultStrategy : public initialGuessStrategy_t {
+public:
+  igDefaultStrategy(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+
+  void FormInitialGuess(occa::memory& o_x, occa::memory& o_rhs);
+  void Update(solver_t &solver, occa::memory& o_x, occa::memory& o_rhs);
+};
+
+// Zero initial guess strategy:  use a zero initial guess.
+class igZeroStrategy : public initialGuessStrategy_t {
+public:
+  igZeroStrategy(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+
+  void FormInitialGuess(occa::memory& o_x, occa::memory& o_rhs);
+  void Update(solver_t &solver, occa::memory& o_x, occa::memory& o_rhs);
+};
+
+// Initial guess strategies based on RHS projection.
+class igProjectionStrategy : public initialGuessStrategy_t {
+protected:
+  dlong curDim;           // Current dimension of the initial guess space
+  dlong maxDim;           // Maximum dimension of the initial guess space
+
+  occa::memory o_btilde;  //  vector (e.g., to be added to space)
+  occa::memory o_xtilde;  // Solution vector corresponding to o_btilde
+  occa::memory o_Btilde;  //  space (orthogonalized)
+  occa::memory o_Xtilde;  // Solution space corresponding to  space
+
+  occa::memory o_w;       // Inverse-degree weights for inner products.
+
+  // temporary buffer for basis inner product output
+  dlong        ctmpNblocks;
+  dfloat       *ctmp;
+  occa::memory o_ctmp;
+
+  dfloat *alphas;         // Buffers for storing inner products.
+  dfloat *alphasThisRank;
+  occa::memory o_alphas;
+
+  occa::kernel igBasisInnerProductsKernel;
+  occa::kernel igReconstructKernel;
+  occa::kernel igScaleKernel;
+  occa::kernel igUpdateKernel;
+
+  void igBasisInnerProducts(occa::memory& o_x, occa::memory& o_Q, occa::memory& o_c, dfloat *c, dfloat *cThisRank);
+  void igReconstruct(occa::memory& o_u, dfloat a, occa::memory& o_c, occa::memory& o_Q, occa::memory& o_unew);
+
+public:
+  igProjectionStrategy(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+  virtual ~igProjectionStrategy();
+
+  virtual void FormInitialGuess(occa::memory& o_x, occa::memory& o_rhs);
+  virtual void Update(solver_t& solver, occa::memory& o_x, occa::memory& o_rhs) = 0;
+};
+
+// "Classic" initial guess strategy from Fischer's 1998 paper.
+class igClassicProjectionStrategy : public igProjectionStrategy {
+public:
+  igClassicProjectionStrategy(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+
+  void Update(solver_t &solver, occa::memory& o_x, occa::memory& o_rhs);
+};
+
+// Rolling QR update for projection history space a la Christensen's thesis.
+class igRollingQRProjectionStrategy : public igProjectionStrategy {
+private:
+  dfloat       *R;   // R factor in QR decomposition (row major)
+  occa::memory o_R;
+
+  occa::kernel igDropQRFirstColumnKernel;
+
+	void givensRotation(dfloat a, dfloat b, dfloat *c, dfloat *s);
+
+public:
+  igRollingQRProjectionStrategy(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+  ~igRollingQRProjectionStrategy();
+
+  void Update(solver_t &solver, occa::memory& o_x, occa::memory& o_rhs);
+};
+
+// Extrapolation initial guess strategy.
+class igExtrapStrategy : public initialGuessStrategy_t {
+private:
+  int Nhistory;
+  int shift;
+  int entry;
+  occa::memory o_xh;
+  occa::memory o_coeffs;
+  occa::kernel igExtrapKernel;
+  occa::kernel igExtrapSparseKernel;
+
+  int Nsparse;
+  occa::memory o_sparseIds;
+  occa::memory o_sparseCoeffs;
+
+  void extrapCoeffs(int m, int M, dfloat *c);
+
+public:
+  igExtrapStrategy(dlong _N, platform_t& _platform, settings_t& _settings, MPI_Comm _comm, int _weighted, occa::memory& _o_weight);
+
+  void FormInitialGuess(occa::memory& o_x, occa::memory& o_rhs);
+  void Update(solver_t &solver, occa::memory& o_x, occa::memory& o_rhs);
+};
+
+// Linear solver with successive-RHS initial-guess generation.
+class initialGuessSolver_t : public linearSolver_t {
+protected:
+  initialGuessStrategy_t* igStrategy;   // The initial guess strategy.
+  linearSolver_t*         linearSolver; // The linearSolver_t that does the solve.
+
+public:
+  initialGuessSolver_t(dlong _N, dlong _Nhalo, platform_t& _platform, settings_t& _settings, MPI_Comm _comm);
+  ~initialGuessSolver_t();
+
+  static initialGuessSolver_t* Setup(dlong _N, dlong _Nhalo,
+                                     platform_t& platform, settings_t& settings, MPI_Comm _comm,
+                                     int _weighted, occa::memory& _o_weight);
+
+  int Solve(solver_t& solver, precon_t& precon,
+            occa::memory& o_x, occa::memory& o_rhs,
+            const dfloat tol, const int MAXIT, const int verbose);
+};
+
+void initialGuessAddSettings(settings_t& settings, const string prefix = "");
+
+#endif /* INITIALGUESS_HPP */