/* Ergo, version 3.8.2, a program for linear scaling electronic structure
 * calculations.
 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
 * and Anastasia Kruchinina.
 * 
 * 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 3 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, see <http://www.gnu.org/licenses/>.
 * 
 * Primary academic reference:
 * Ergo: An open-source program for linear-scaling electronic structure
 * calculations,
 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
 * Kruchinina,
 * SoftwareX 7, 107 (2018),
 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
 * 
 * For further information about Ergo, see <http://www.ergoscf.org>.
 */

/** @file LanczosLargestMagnitudeEig.h Class for computing the largest 
 * magnitude eigenvalue of a symmetric matrix with the Lanczos method.
 *
 * Copyright(c) Emanuel Rubensson 2007
 *
 * @author Emanuel Rubensson
 * @date February 2007
 *
 */
#ifndef MAT_LANCZOSLARGESTMAGNITUDEEIG
#define MAT_LANCZOSLARGESTMAGNITUDEEIG
#include <limits>
#include "Lanczos.h"
namespace mat { /* Matrix namespace */
  namespace arn { /* Arnoldi type methods namespace */
    template<typename Treal, typename Tmatrix, typename Tvector>
      class LanczosLargestMagnitudeEig 
      : public Lanczos<Treal, Tmatrix, Tvector> {
    public:
      LanczosLargestMagnitudeEig(Tmatrix const & AA, Tvector const & startVec, 
				 int maxIter = 100, int cap = 100) 
	: Lanczos<Treal, Tmatrix, Tvector>(AA, startVec, maxIter, cap), 
	eVal(0), acc(0), eigVectorTri(0), absTol( template_blas_get_num_limit_max<Treal>() ),
	relTol(template_blas_sqrt(template_blas_sqrt(getRelPrecision<Treal>()))), 
	eValTmp(0), accTmp(0) {}
      void setRelTol(Treal const newTol) { relTol = newTol; }
      void setAbsTol(Treal const newTol) { absTol = newTol; }

	inline void getLargestMagnitudeEig(Treal& ev, Treal& accuracy) {
	  ev = eVal;
	  accuracy = acc;
	}
	void getLargestMagnitudeEigPair(Treal& eValue, 
					Tvector& eVector, 
					Treal& accuracy);
	
	virtual void run() {
	  Lanczos<Treal, Tmatrix, Tvector>::run();
	  computeEigVec();
	  eVal = eValTmp;
	  acc = accTmp;
	  rerun();
	  /* FIXME! Elias added these extra Lanczos reruns for small
	     matrices to make the tests pass. This is bad. 
             Elias note 2011-01-19: now added one more rerun()
             because otherwise the test failed when run on Mac. 
             This is really bad. */
	  if(this->A.get_nrows() == 5) {
	    rerun();
	    rerun();
            rerun();
	  }
	}
	
	void rerun() {
#if 1
	  /* Because of the statistical chance of misconvergence: 
	   * Compute new eigenpair with eigenvector in direction orthogonal 
	   * to the first eigenvector.
	   */
	  Tvector newResidual(eVec);
	  newResidual.rand();
	  Treal sP = transpose(eVec) * newResidual;
	  newResidual += -sP * eVec;
	  this->restart(newResidual);
	  Lanczos<Treal, Tmatrix, Tvector>::run();
	  /* If the new eigenvalue has larger magnitude: 
	   * Use those instead
 	   */
	  if (template_blas_fabs(eValTmp) > template_blas_fabs(eVal)) {
	    eVal = eValTmp;
	    acc = accTmp;
	    computeEigVec();
	  }
	  // Guard against unrealistically small accuracies
	  Treal machine_epsilon = mat::getMachineEpsilon<Treal>();
	  acc = acc / template_blas_fabs(eVal) > 2 * machine_epsilon ? acc : 2 * template_blas_fabs(eVal) * machine_epsilon;
#endif
	}

	virtual ~LanczosLargestMagnitudeEig() {
	  delete[] eigVectorTri;
	}
    protected:
	Treal eVal;
	Tvector eVec;
	Treal acc;
	Treal* eigVectorTri; /** Eigenvector to the tridiagonal matrix
			      *  length: this->j
			      */	
	Treal absTol;
	Treal relTol;
	void computeEigenPairTri();
	void computeEigVec();
	virtual void update() {
	  computeEigenPairTri();
	}
	virtual bool converged() const;
	Treal eValTmp;
	Treal accTmp;
    private:
    };
    
    template<typename Treal, typename Tmatrix, typename Tvector>
      void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
      getLargestMagnitudeEigPair(Treal& eValue, 
				 Tvector& eVector, 
				 Treal& accuracy) {
      eValue = eVal;
      accuracy = acc;
      eVector = eVec;
    }


    template<typename Treal, typename Tmatrix, typename Tvector>
      void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
      computeEigenPairTri() {
      delete[] eigVectorTri;
      Treal* eigVectorMax = new Treal[this->j];
      Treal* eigVectorMin = new Treal[this->j];      
      
      /* Get largest eigenvalue. */
      int const lMax = this->j - 1;
      Treal eValMax;
      Treal accMax;
      this->Tri.getEigsByIndex(&eValMax, eigVectorMax, &accMax,  
			       lMax, lMax);
      /* Get smallest eigenvalue. */
      int const lMin = 0;
      Treal eValMin;
      Treal accMin;
      this->Tri.getEigsByIndex(&eValMin, eigVectorMin, &accMin,  
			       lMin, lMin);
      if (template_blas_fabs(eValMin) > template_blas_fabs(eValMax)) {
	eValTmp = eValMin;
	accTmp = accMin;
	delete[] eigVectorMax;
	eigVectorTri = eigVectorMin;
      }
      else {
	eValTmp = eValMax;
	accTmp = accMax;
      	delete[] eigVectorMin;
	eigVectorTri = eigVectorMax;
      }
    }
    

    template<typename Treal, typename Tmatrix, typename Tvector>
      void LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
      computeEigVec() {
      /* Compute eigenvector as a linear combination of Krylov vectors */
      assert(eigVectorTri);
      this->getEigVector(eVec, eigVectorTri);      
    }

    
    template<typename Treal, typename Tmatrix, typename Tvector>
      bool LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::
      converged() const {
      bool conv = accTmp < absTol;                 /* Absolute accuracy */ 
      if (template_blas_fabs(eValTmp) > 0) {
	conv = conv && 
	  accTmp / template_blas_fabs(eValTmp) < relTol; /* Relative acc.*/
      }
      return conv;
    }
    



#if 1
    template<typename Treal, typename Tmatrix, typename Tvector>
      class LanczosLargestMagnitudeEigIfSmall 
      : public LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> {
    public:
      LanczosLargestMagnitudeEigIfSmall
	(Tmatrix const & AA, Tvector const & startVec, 
	 Treal const maxAbsVal,
	 int maxIter = 100, int cap = 100) 
	: LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector> 
	(AA, startVec, maxIter, cap), maxAbsValue(maxAbsVal) {
      }
        bool largestMagEigIsSmall() {return eigIsSmall;}
      
	virtual void run() {
	  Lanczos<Treal, Tmatrix, Tvector>::run();
	  this->computeEigVec();
	  this->eVal = this->eValTmp;
	  this->acc = this->accTmp;
	  if (eigIsSmall) /* No need to rerun if eigenvalue is large. */
	    this->rerun();
	}

    protected:
	Treal const maxAbsValue;
	bool eigIsSmall;
	virtual void update() {
	  LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::update();
	  eigIsSmall = template_blas_fabs(this->eValTmp) < maxAbsValue;  
	}
	virtual bool converged() const;
    private:
    };    

    template<typename Treal, typename Tmatrix, typename Tvector>
      bool LanczosLargestMagnitudeEigIfSmall<Treal, Tmatrix, Tvector>::
      converged() const {
      bool convAccuracy = 
	LanczosLargestMagnitudeEig<Treal, Tmatrix, Tvector>::converged();
      return convAccuracy || (!eigIsSmall); 
    }

#endif

  } /* end namespace arn */

  //// Utility function

  // EMANUEL COMMENT:
  // A function similar to euclIfSmall below lies/used to lie inside the MatrixSymmetric class.
  // There, the matrix was copied and truncated before the calculation.
  // It is unclear if that had a significant positive impact on the execution time - it definitely required more memory. 
  /** Returns interval containing the Euclidean norm of the matrix M.
   *  If it is smaller than 'maxAbsVal' it is computed with 
   *  the 'requestedAccuracy'.  
   *  If the norm is larger than 'maxAbsVal', the returned interval is based on 
   *  the Frobenius norm.
   */
  template<typename Tmatrix, typename Treal>
    Interval<Treal> euclIfSmall(Tmatrix const & M,
				Treal const requestedAbsAccuracy,
				Treal const requestedRelAccuracy,
				Treal const maxAbsVal,
				typename Tmatrix::VectorType * const eVecPtr = 0,
				int maxIter = -1) {
    assert(requestedAbsAccuracy >= 0);
    //    Treal frobTmp;
    /* Check if norm is really small, or can easily be determined to be 'large', in those cases quick return */
    Treal euclLowerBound;
    Treal euclUpperBound;
    M.quickEuclBounds(euclLowerBound, euclUpperBound);
    if ( euclUpperBound < requestedAbsAccuracy )
      // Norm is really small, quick return     
      return Interval<Treal>( euclLowerBound, euclUpperBound );
    if ( euclLowerBound > maxAbsVal )
      // Norm is not small, quick return
      return Interval<Treal>( euclLowerBound, euclUpperBound );
    if(maxIter == -1)
      maxIter = M.get_nrows() * 100;
    typename Tmatrix::VectorType guess;
    SizesAndBlocks cols;
    M.getCols(cols);
    guess.resetSizesAndBlocks(cols);
    guess.rand();
    mat::arn::LanczosLargestMagnitudeEigIfSmall<Treal, Tmatrix, typename Tmatrix::VectorType>
      lan(M, guess, maxAbsVal, maxIter);
    lan.setAbsTol( requestedAbsAccuracy );
    lan.setRelTol( requestedRelAccuracy );
    lan.run();
    Treal eVal = 0;
    Treal acc = 0;
    Treal eValMin = 0;
    if (lan.largestMagEigIsSmall()) {
      if (eVecPtr)
	lan.getLargestMagnitudeEigPair(eVal, *eVecPtr, acc); 
      else
	lan.getLargestMagnitudeEig(eVal, acc); 
      eVal = template_blas_fabs(eVal);
      eValMin = eVal - acc;
      eValMin = eValMin >= 0 ? eValMin : (Treal)0;
      return Interval<Treal>(eValMin, eVal + acc);
    }
    else {
      eValMin = euclLowerBound;
      eValMin = eValMin >= maxAbsVal ? eValMin : maxAbsVal; 
      return Interval<Treal>(eValMin, euclUpperBound);
    }
  }
  
} /* end namespace mat */
#endif