GeneralEvenOddRationalRatioMixedPrec.h

Go to the documentation of this file.

    /*************************************************************************************
 
    Grid physics library, www.github.com/paboyle/Grid 
 
    Source file: ./lib/qcd/action/pseudofermion/GeneralEvenOddRationalRatioMixedPrec.h
 
    Copyright (C) 2015
 
    Author: Christopher Kelly <ckelly@bnl.gov>
    Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
    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.
 
    See the full license in the file "LICENSE" in the top level distribution directory
    *************************************************************************************/
    /*  END LEGAL */
#ifndef QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
#define QCD_PSEUDOFERMION_GENERAL_EVEN_ODD_RATIONAL_RATIO_MIXED_PREC_H
 
#include <Grid/algorithms/iterative/ConjugateGradientMultiShiftCleanup.h>
 
NAMESPACE_BEGIN(Grid);

    // Generic rational approximation for ratios of operators utilizing the mixed precision multishift algorithm
    // cf. GeneralEvenOddRational.h for details
      
    template<class ImplD, class ImplF>
    class GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction : public GeneralEvenOddRatioRationalPseudoFermionAction<ImplD> {
    private:
      typedef typename ImplD::FermionField FermionFieldD;
      typedef typename ImplF::FermionField FermionFieldF;
 
      FermionOperator<ImplD> & NumOpD;
      FermionOperator<ImplD> & DenOpD;
 
      FermionOperator<ImplF> & NumOpF;
      FermionOperator<ImplF> & DenOpF;
 
      Integer ReliableUpdateFreq;
    protected:
 
      //Action evaluation
      //Allow derived classes to override the multishift CG
      virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, FermionFieldD &out){
#if 1
    SchurDifferentiableOperator<ImplD> schurOp(numerator ? NumOpD : DenOpD);
    ConjugateGradientMultiShift<FermionFieldD> msCG(MaxIter, approx);
    msCG(schurOp,in, out);
#else
    SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
    SchurDifferentiableOperator<ImplF> schurOpF(numerator ? NumOpF : DenOpF);
    FermionFieldD inD(NumOpD.FermionRedBlackGrid());
    FermionFieldD outD(NumOpD.FermionRedBlackGrid());
 
    // Action better with higher precision?
    ConjugateGradientMultiShiftMixedPrec<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
    msCG(schurOpD, in, out);
#endif
      }
      //Force evaluation
      virtual void multiShiftInverse(bool numerator, const MultiShiftFunction &approx, const Integer MaxIter, const FermionFieldD &in, std::vector<FermionFieldD> &out_elems, FermionFieldD &out){
    SchurDifferentiableOperator<ImplD> schurOpD(numerator ? NumOpD : DenOpD);
    SchurDifferentiableOperator<ImplF>  schurOpF(numerator ? NumOpF  : DenOpF);
 
    FermionFieldD inD(NumOpD.FermionRedBlackGrid());
    FermionFieldD outD(NumOpD.FermionRedBlackGrid());
    std::vector<FermionFieldD> out_elemsD(out_elems.size(),NumOpD.FermionRedBlackGrid());
    ConjugateGradientMultiShiftMixedPrecCleanup<FermionFieldD, FermionFieldF> msCG(MaxIter, approx, NumOpF.FermionRedBlackGrid(), schurOpF, ReliableUpdateFreq);
    msCG(schurOpD, in, out_elems, out);
      }
      //Allow derived classes to override the gauge import
      virtual void ImportGauge(const typename ImplD::GaugeField &Ud){
 
    typename ImplF::GaugeField Uf(NumOpF.GaugeGrid());
    precisionChange(Uf, Ud);
 
    std::cout << "Importing "<<norm2(Ud)<<" "<< norm2(Uf)<<" " <<std::endl;
    
    NumOpD.ImportGauge(Ud);
    DenOpD.ImportGauge(Ud);
 
    NumOpF.ImportGauge(Uf);
    DenOpF.ImportGauge(Uf);
      }
      
    public:
      GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction(FermionOperator<ImplD>  &_NumOpD, FermionOperator<ImplD>  &_DenOpD, 
                                  FermionOperator<ImplF>  &_NumOpF, FermionOperator<ImplF>  &_DenOpF, 
                                  const RationalActionParams & p, Integer _ReliableUpdateFreq
                                  ) : GeneralEvenOddRatioRationalPseudoFermionAction<ImplD>(_NumOpD, _DenOpD, p),
                                  ReliableUpdateFreq(_ReliableUpdateFreq),
                                  NumOpD(_NumOpD), DenOpD(_DenOpD),
                                  NumOpF(_NumOpF), DenOpF(_DenOpF)
      {}
      
      virtual std::string action_name(){return "GeneralEvenOddRatioRationalMixedPrecPseudoFermionAction";}
    };
 
NAMESPACE_END(Grid);
 
#endif