Davidson.h

/*
   CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
   Copyright (C) 2013-2016 Sebastian Wouters

   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.
*/

#ifndef DAVIDSON_CHEMPS2_H
#define DAVIDSON_CHEMPS2_H

namespace CheMPS2{
   class Davidson{

      public:


         Davidson( const int veclength, const int MAX_NUM_VEC, const int NUM_VEC_KEEP, const double RTOL, const double DIAG_CUTOFF, const bool debug_print, const char problem_type = 'E' );

         virtual ~Davidson();


         char FetchInstruction( double ** pointers );


         int GetNumMultiplications() const;

      private:

         int veclength; // The vector length
         int nMultiplications; // Current number of requested matrix-vector multiplications
         char state; // Current state of the algorithm --> based on this parameter the next instruction is given
         bool debug_print;
         char problem_type;

         // Davidson parameters
         int MAX_NUM_VEC;
         int NUM_VEC_KEEP;
         double DIAG_CUTOFF;
         double RTOL;

         // To store the vectors and the matrix x vectors
         int num_vec;
         double ** vecs;
         double ** Hvecs;
         int num_allocated;

         // The effective diagonalization problem
         double * mxM;
         double * mxM_eigs;
         double * mxM_vecs;
         int mxM_lwork;
         double * mxM_work;
         double * mxM_rhs;

         // Vector spaces
         double * t_vec;
         double * u_vec;
         double * work_vec;
         double * diag;
         double * RHS;

         // For the deflation
         double * Reortho_Lowdin;
         double * Reortho_Overlap_eigs;
         double * Reortho_Overlap;
         double * Reortho_Eigenvecs;

         // Control script functions
         double FrobeniusNorm( double * current_vector );
         void SafetyCheckGuess();
         void AddNewVec();
         double DiagonalizeSmallMatrixAndCalcResidual(); // Returns the residual norm
         void CalculateNewVec();
         void Deflation();
         void MxMafterDeflation();
         void SolveLinearSystemDeflation( const int NUM_SOLUTIONS );

   };
}

#endif