AmIteratorTmpl.h

#ifndef PSCF_AM_ITERATOR_TMPL_H
#define PSCF_AM_ITERATOR_TMPL_H
 
/*
* PSCF - Polymer Self-Consistent Field
*
* Copyright 2015 - 2025, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <util/containers/DArray.h>       // member template
#include <util/containers/DMatrix.h>      // member template
#include <util/containers/RingBuffer.h>   // member template
#include <util/accumulators/Average.h>    // member template
#include <util/misc/Timer.h>              // member
 
// Uncomment to test details of Anderson-Mixing algorithm performance
//#define PSCF_AM_TEST
 
namespace Pscf {
 
   using namespace Util;

   template <typename Iterator, typename T>
   class AmIteratorTmpl : virtual public Iterator
   {
   public:

      AmIteratorTmpl();

      ~AmIteratorTmpl();

      virtual void readParameters(std::istream& in);

      int solve(bool isContinuation = false);

      void outputTimers(std::ostream& out) const;

      void clearTimers();
 
   protected:

      double epsilon_;

      int maxItr_;

      int maxHist_;

      int verbose_;

      std::string errorType_;
 
      // Parameter initialization

      void readErrorType(std::istream& in);

      virtual bool isValidErrorType();

      void readMixingParameters(std::istream& in, 
                                bool useLambdaRamp = true);
 
      // Protected AM mixing operations

      void allocateAM();

      bool isAllocatedAM() const;

      virtual void setup(bool isContinuation);

      virtual void clear();

      virtual 
      double computeError(T& residTrial, T& stateTrial,
                          std::string errorType, int verbose);

      double computeError(int verbose);

      double lambdaRampFactor();

      virtual double computeLambda();
 
      // Protected accessors for member variables

      int verbose() const;

      std::string errorType() const;

      T const & residual() const;

      T const & state() const;

      int totalItr();
 
      // --- Timer value accessors --- //

      double timerTotal();

      double timerMDE();

      double timerAM();

      double timerResid();

      double timerError();

      double timerCoeff();

      double timerOmega();
 
      // Inherited members of parent classes with non-dependent names
      using ParamComposite::read;
      using ParamComposite::readOptional;
 
   private:
 
      // Private member variables

      double error_;

      int nBasis_;

      int itr_;

      int totalItr_;

      int nElem_;

      double lambda_;

      double r_;

      bool useLambdaRamp_;

      bool isAllocatedAM_;

      RingBuffer<T> stateHistory_;

      RingBuffer<T> stateBasis_;

      RingBuffer<T> residualHistory_;

      RingBuffer<T> residualBasis_;

      DMatrix<double> U_;

      DArray<double> v_;

      DArray<double> coeffs_;

      T stateTrial_;

      T residualTrial_;

      T temp_;
 
      // Timers for analyzing performance
      Timer timerMDE_;
      Timer timerAM_;
      Timer timerResid_;
      Timer timerError_;
      Timer timerCoeff_;
      Timer timerOmega_;
      Timer timerTotal_;
 
      #ifdef PSCF_AM_TEST
      bool hasAmTest_{false};
      double preError_{0};
      double projectionError_{0};
      double correctionError_{0};
      double projectionRatio_{0};
      double correctionRatio_{0};
      int testCounter{0};
      #endif
 
      // Private non-virtual functions used in AM algorithm

      void updateBasis(RingBuffer<T> & basis,
                       RingBuffer<T> const & history);

      void updateU(DMatrix<double> & U,
                   RingBuffer<T> const & resBasis);

      void computeV(DArray<double> & v, 
                   T const & resCurrent,
                   RingBuffer<T> const & resBasis);

      void computeTrialCoeff();

      void updateTrial();

      void addEqVectors(T& v,
                        RingBuffer<T> const & basis,
                        DArray<double> coeffs);
 
      // Private virtual functions with a default implementation 

      virtual 
      void addCorrection(T& stateTrial, T const & residualTrial);
 
      // Private pure virtual functions that interact with parent system

      virtual int nElements() = 0;

      virtual bool hasInitialGuess() = 0;

      virtual void getCurrent(T& curr) = 0;

      virtual void evaluate() = 0;

      virtual void getResidual(T& resid) = 0;

      virtual void update(T& newGuess) = 0;

      virtual void outputToLog() = 0;
 
      // Private functions for vector math

      virtual void setEqual(T& a, T const & b) = 0;

      virtual double dotProduct(T const & a, T const & b) = 0;

      double norm(T const & a);

      virtual double maxAbs(T const & a) = 0;

      virtual void subVV(T& a, T const & b, T const & c) = 0;

      virtual void addEqVc(T& a, T const & b, double c) = 0;
 
   };
 
   /*
   * Return integer level for verbosity of the log output (0-2).
   */
   template <typename Iterator, typename T>
   int AmIteratorTmpl<Iterator,T>::verbose() const
   {  return verbose_; }
 
   /*
   * Return error type string.
   */
   template <typename Iterator, typename T>
   std::string AmIteratorTmpl<Iterator,T>::errorType() const
   {  return errorType_; }
 
   /*
   * Return the current state vector by const reference.
   */
   template <typename Iterator, typename T>
   T const & AmIteratorTmpl<Iterator,T>::state() const
   {  return stateHistory_[0]; }
 
   /*
   * Return the current residual vector by const reference.
   */
   template <typename Iterator, typename T>
   T const & AmIteratorTmpl<Iterator,T>::residual() const
   {  return residualHistory_[0]; }
 
   /*
   * Has memory required by the AM algorithm been allocated?
   */
   template <typename Iterator, typename T>
   bool AmIteratorTmpl<Iterator,T>::isAllocatedAM() const
   {  return isAllocatedAM_; }
 
   /*
   * Return total iteration counter.
   */
   template <typename Iterator, typename T>
   int AmIteratorTmpl<Iterator,T>::totalItr()
   {  return totalItr_; }
 
   /*
   * Return computing MDE time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerMDE()
   {  return timerMDE_.time(); }
 
   /*
   * Return computing AM time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerAM()
   {  return timerAM_.time(); }
 
   /*
   * Return computing Resid time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerResid()
   {  return timerResid_.time(); }
 
   /*
   * Return computing Error time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerError()
   {  return timerError_.time(); }
 
   /*
   * Return computing Coeff time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerCoeff()
   {  return timerCoeff_.time(); }
 
   /*
   * Return computing Omega time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerOmega()
   {  return timerOmega_.time(); }
 
   /*
   * Return total time cost.
   */
   template <typename Iterator, typename T>
   double AmIteratorTmpl<Iterator,T>::timerTotal()
   {  return timerTotal_.time(); }
 
}
#include "AmIteratorTmpl.tpp"
#endif