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
 
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);
 
   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(int verbose);
 
      T const & residual() const;
 
      T const & field() const;
 
      int verbose() const;
 
      bool isAllocatedAM() const;
 
   private:
 
      // Private member variables
 
      double epsilon_;
 
      double lambda_;
 
      int maxItr_;
 
      int maxHist_;
 
      int nBasis_;
 
      int itr_;
 
      int nElem_; 
 
      int verbose_;
 
      bool outputTime_;
 
      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_;
 
      // --- 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_; }
 
}
#include "AmIteratorTmpl.tpp"
#endif