CRAMSolver.hxx

Go to the documentation of this file.

/*
  This file is part of MURE,
  Copyright (C) 2007-2021 MURE developers.
 
  MURE is free software: you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.
 
  MURE 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 Lesser General Public License for more details.
 
  You should have received a copy of the GNU Lesser General Public License
  along with MURE.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef _CRAMSolver_
#define _CRAMSolver_
 
#include <cstdlib>
#include <cstring>
#include <complex>
#include <vector>
#include <map>
#include <algorithm>
#include <iosfwd>
#include <new>
#include <utility>
 
#ifdef TEST
#include <iostream>
#else
#include "MUREGlobal.hxx"
#endif
#include "BatemanSolver.hxx"
#include "Timer.hxx"
 
#include <Eigen/Sparse>
#ifdef EIGEN_SuperLU_ColMajor
  #include <Eigen/SuperLUSupport>
  typedef Eigen::SparseMatrix<std::complex<double>, Eigen::ColMajor> t_SpMatrix;
  typedef vector<Eigen::SuperLU<t_SpMatrix > > t_SpMat_SolverVector;
#elif EIGEN_UmfPackLU_ColMajor
  #include <Eigen/UmfPackSupport>
  typedef Eigen::SparseMatrix<std::complex<double>, Eigen::ColMajor > t_SpMatrix;
  typedef vector< Eigen::UmfPackLU<t_SpMatrix > > t_SpMat_SolverVector;
#else
  // EIGEN_BiCGSTAB_RowMajor
  #include<Eigen/IterativeLinearSolvers>
  typedef Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor > t_SpMatrix;
  typedef vector< Eigen::BiCGSTAB<t_SpMatrix > > t_SpMat_SolverVector;
#endif
 
typedef Eigen::Triplet<double> t_triplet;
 
using namespace std;
 
//-----------------------------------------------------------------------------//
 
struct t_CRAMCoeff
{
    int CramOrderDiv2;
    double alpha0;
    vector<complex<double> > theta;
    vector<complex<double> > alpha;
};
// matrix exponential method
//________________________________________________________________________
 
class CRAMSolver : public BatemanSolver
{
  public:
    CRAMSolver(string type = "IPF", int order = 48, bool aut = true);   
    // CRAMSolver(const CRAMSolver& BS);     //!< Copy constructor
    // ~CRAMSolver();                   //!< Destructor
 
    void Solve(double *N0, double **Matrix, int N, double t1, double t2) override;
    void SetTheMatrix(double **f)
    {
        fTheMatrix = f;
    }
    void printCoeffs(int order = 0);
    void SetNSubsteps(int n)
    {
        nSubsteps = n;
        automatic = false;
    };
    void SetOrder(int n)
    {
        CRAMOrder = n;
        automatic = false;
    };
 
  protected:
 
  private:
    void BuildEqns(double *N0, double **Matrix); 
    t_CRAMCoeff GetCRAMCoefficients(int order = 0);
    t_SpMat_SolverVector GetSparseLUSolvers();
    // void IterateForSubsteps(double* N0, double** fMatrix);
    int GetNumberSubsteps(double &N0_i, double &A_ii, double dT, int order, double sigDigits);
    double GetSignificantDigits(double &N0_i, double &A_ii, double dT, int order, int niter);
    int GetNeededCRAMOrder(double &N0_i, double &A_ii, double dT, double sigDigits);
    // void GetNumberSubsteps(double& N0_i, double& A_ii, int& numberSubsteps, int& order);
    // int GetNumberSubsteps(double& N0_i, double& A_ii);
    void SetOptimalStrategy(double *N0, double **fMatrix);
    // int GetNeededCRAMOrder(double& N0_i, double& A_ii);
    int GetNearestCRAMOrder(int order);
    double GetLG10(double &n0, double &A_ii, double dT, int nT);
    void CleanUp();
 
    t_SpMat_SolverVector CRAMSpMatSolver;
    t_CRAMCoeff CRAMCoefficients;
    string CRAMType;
    string SparseSolver;
    int CRAMOrder;
    double      **fTheMatrix;       
    double CRAMDeltaT;
    double CRAMDeltaTsubstep;
    int nSubsteps;
    Timer HPTimer; // high precision C++ timer
    const double cutoff = 1e-20;
    const double significantDigits = 10;
    const int maxorder = 48;
    bool automatic = false;
    bool fBuildEqns;
    static const vector<int> fCramOrder;
    static const vector<t_CRAMCoeff> CRAMCoeff;
};
#endif