AmIteratorTmpl.h

#ifndef PSCF_AM_ITERATOR_TMPL_H
#define PSCF_AM_ITERATOR_TMPL_H
 
/*
* PSCF - Polymer Self-Consistent Field Theory
*
* Copyright 2016 - 2022, 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/misc/Timer.h>            // member template
#include <util/accumulators/Average.h>  // member template
 
// 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();
 
      void readParameters(std::istream& in);
 
      int solve(bool isContinuation = false);
      
      void outputTimers(std::ostream& out);
      
      void clearTimers();
      
      std::string errorType() const;
 
   protected:
 
      std::string errorType_;
 
      // Members of parent classes with non-dependent names
      using Iterator::setClassName;
      using ParamComposite::read;
      using ParamComposite::readOptional;
 
      void setMaxItr(int maxItr);
 
      void setMaxHist(int maxHist);
 
      void setErrorType(std::string errorType);
 
      void readErrorType(std::istream& in);
 
      virtual bool isValidErrorType();
 
      virtual double norm(T const & hist);
 
      void allocateAM();
 
      virtual void clear();
 
      virtual void setup(bool isContinuation);
      
      virtual double computeError(T& residTrial, T& fieldTrial, 
                                  std::string errorType, 
                                  int verbose);
     
      double computeError(int verbose);
      
      virtual double computeLambda(double r);
 
      T const & residual() const;
 
      T const & field() const;
 
      int verbose() const;
      
      int totalItr();
      
      double timerTotal();
      
      double timerMDE();
 
      double timerAM();
 
      double timerResid();
 
      double timerError();
 
      double timerCoeff();
 
      double timerOmega();
 
      bool isAllocatedAM() const;
      
   private:
      
      // Private member variables
 
      double error_;
 
      double epsilon_;
 
      double lambda_;
      
      double r_;
 
      int maxItr_;
 
      int maxHist_;
 
      int nBasis_;
 
      int itr_;
      
      int totalItr_;
 
      int nElem_; 
 
      int verbose_;
 
      bool isAllocatedAM_;
 
      RingBuffer<T> fieldHists_;
 
      RingBuffer<T> fieldBasis_;
 
      RingBuffer<T> resHists_;
 
      RingBuffer<T> resBasis_;
 
      DMatrix<double> U_;
 
      DArray<double> coeffs_;
 
      DArray<double> v_;
 
      T fieldTrial_;
 
      T resTrial_;
 
      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
      
      // --- Non-virtual private functions (implemented here) ---- //
 
      void computeResidCoeff();
 
      void updateGuess();
 
      // --- Private virtual functions with default implementations --- //
 
      virtual void updateU(DMatrix<double> & U, 
                           RingBuffer<T> const & resBasis, int nHist);
 
      virtual void updateV(DArray<double> & v, T const & resCurrent, 
                           RingBuffer<T> const & resBasis, int nHist);
      
      // --- Pure virtual methods for doing AM iterator math --- //
 
      virtual void setEqual(T& a, T const & b) = 0;
 
      virtual double dotProduct(T const & a, T const & b) = 0;
 
      virtual double maxAbs(T const & hist) = 0;
 
      virtual 
      void updateBasis(RingBuffer<T> & basis, 
                       RingBuffer<T> const & hists) = 0;
 
      virtual 
      void addHistories(T& trial, 
                        RingBuffer<T> const & basis, 
                        DArray<double> coeffs, 
                        int nHist) = 0;
 
      virtual 
      void addPredictedError(T& fieldTrial, T const & resTrial, 
                             double lambda) = 0;
 
      // -- Pure virtual functions to exchange data with parent system -- //
     
      virtual bool hasInitialGuess() = 0;
     
      virtual int nElements() = 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;
 
   };
  
   /*
   * Return the current residual vector by const reference.
   */
   template <typename Iterator, typename T>
   T const & AmIteratorTmpl<Iterator,T>::residual() const
   {  return resHists_[0]; }
 
   /*
   * Return the current field/state vector by const reference.
   */
   template <typename Iterator, typename T>
   T const & AmIteratorTmpl<Iterator,T>::field() const
   {  return fieldHists_[0]; }
 
   /*
   * Integer level for verbosity of the log output (0-2).
   */
   template <typename Iterator, typename T>
   int AmIteratorTmpl<Iterator,T>::verbose() const
   {  return verbose_; }
 
   /*
   * Has memory required by AM algorithm been allocated?
   */
   template <typename Iterator, typename T>
   bool AmIteratorTmpl<Iterator,T>::isAllocatedAM() const
   {  return isAllocatedAM_; }
   
   /*
   * Return error type
   */ 
   template <typename Iterator, typename T>
   std::string AmIteratorTmpl<Iterator,T>::errorType() const
   {  return errorType_; }
   
   /*
   * 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