mirror of
https://github.com/paboyle/Grid.git
synced 2024-09-20 09:15:38 +01:00
216 lines
9.3 KiB
C++
216 lines
9.3 KiB
C++
#ifndef GRID_ALGORITHM_LINEAR_OP_H
|
|
#define GRID_ALGORITHM_LINEAR_OP_H
|
|
|
|
#include <Grid.h>
|
|
|
|
namespace Grid {
|
|
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
// LinearOperators Take a something and return a something.
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Hopefully linearity is satisfied and the AdjOp is indeed the Hermitian conjugate (transpose if real):
|
|
//
|
|
// i) F(a x + b y) = aF(x) + b F(y).
|
|
// ii) <x|Op|y> = <y|AdjOp|x>^\ast
|
|
//
|
|
// Would be fun to have a test linearity & Herm Conj function!
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
template<class Field> class LinearOperatorBase {
|
|
public:
|
|
virtual void Op (const Field &in, Field &out) = 0; // Abstract base
|
|
virtual void AdjOp (const Field &in, Field &out) = 0; // Abstract base
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Hermitian operators are self adjoint and only require Op to be defined, so refine the base
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
template<class Field> class HermitianOperatorBase : public LinearOperatorBase<Field> {
|
|
public:
|
|
void AdjOp(const Field &in, Field &out) {
|
|
return Op(in,out);
|
|
};
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Interface defining what I expect of a general sparse matrix
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
template<class Field> class SparseMatrixBase {
|
|
public:
|
|
// Full checkerboar operations
|
|
virtual void M (const Field &in, Field &out);
|
|
virtual void Mdag (const Field &in, Field &out);
|
|
virtual void MdagM(const Field &in, Field &out);
|
|
|
|
// half checkerboard operaions
|
|
virtual void Mpc (const Field &in, Field &out);
|
|
virtual void MpcDag (const Field &in, Field &out);
|
|
virtual void MpcDagMpc(const Field &in, Field &out);
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
// Whereas non hermitian takes a generic sparse matrix (e.g. lattice action)
|
|
// conforming to sparse matrix interface and builds the full checkerboard non-herm operator
|
|
// Op and AdjOp distinct.
|
|
// By sharing the class for Sparse Matrix across multiple operator wrappers, we can share code
|
|
// between RB and non-RB variants. Sparse matrix is like the fermion action def, and then
|
|
// the wrappers implement the specialisation of "Op" and "AdjOp" to the cases minimising
|
|
// replication of code.
|
|
/////////////////////////////////////////////////////////////////////////////////////////////
|
|
template<class SparseMatrix,class Field>
|
|
class NonHermitianOperator : public LinearOperatorBase<Field> {
|
|
SparseMatrix &_Mat;
|
|
public:
|
|
NonHermitianOperator(SparseMatrix &Mat): _Mat(Mat){};
|
|
void Op (const Field &in, Field &out){
|
|
_Mat.M(in,out);
|
|
}
|
|
void AdjOp (const Field &in, Field &out){
|
|
_Mat.Mdag(in,out);
|
|
}
|
|
};
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////
|
|
// Redblack Non hermitian wrapper
|
|
////////////////////////////////////////////////////////////////////////////////////
|
|
template<class SparseMatrix,class Field>
|
|
class NonHermitianRedBlackOperator : public LinearOperatorBase<Field> {
|
|
SparseMatrix &_Mat;
|
|
public:
|
|
NonHermitianRedBlackOperator(SparseMatrix &Mat): _Mat(Mat){};
|
|
void Op (const Field &in, Field &out){
|
|
this->Mpc(in,out);
|
|
}
|
|
void AdjOp (const Field &in, Field &out){ //
|
|
this->MpcDag(in,out);
|
|
}
|
|
};
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////
|
|
// Hermitian wrapper
|
|
////////////////////////////////////////////////////////////////////////////////////
|
|
template<class SparseMatrix,class Field>
|
|
class HermitianOperator : public HermitianOperatorBase<Field> {
|
|
SparseMatrix &_Mat;
|
|
public:
|
|
HermitianOperator(SparseMatrix &Mat): _Mat(Mat) {};
|
|
|
|
void Op (const Field &in, Field &out){
|
|
this->Mpc(in,out);
|
|
}
|
|
void AdjOp (const Field &in, Field &out){ //
|
|
this->MpcDag(in,out);
|
|
}
|
|
};
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////
|
|
// Hermitian RedBlack wrapper
|
|
////////////////////////////////////////////////////////////////////////////////////
|
|
template<class SparseMatrix,class Field>
|
|
class HermitianRedBlackOperator : public HermitianOperatorBase<Field> {
|
|
SparseMatrix &_Mat;
|
|
public:
|
|
HermitianRedBlackOperator(SparseMatrix &Mat): _Mat(Mat) {};
|
|
void Op (const Field &in, Field &out){
|
|
this->MpcDagMpc(in,out);
|
|
}
|
|
};
|
|
|
|
|
|
/////////////////////////////////////////////////////////////
|
|
// Base classes for functions of operators
|
|
/////////////////////////////////////////////////////////////
|
|
template<class Field> class OperatorFunction {
|
|
public:
|
|
virtual void operator() (LinearOperatorBase<Field> *Linop, const Field &in, Field &out) = 0;
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////
|
|
// Base classes for polynomial functions of operators ? needed?
|
|
/////////////////////////////////////////////////////////////
|
|
template<class Field> class OperatorPolynomial : public OperatorFunction<Field> {
|
|
public:
|
|
virtual void operator() (LinearOperatorBase<Field> *Linop,const Field &in, Field &out) = 0;
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////
|
|
// Base classes for iterative processes based on operators
|
|
// single input vec, single output vec.
|
|
/////////////////////////////////////////////////////////////
|
|
template<class Field> class IterativeProcess : public OperatorFunction<Field> {
|
|
public:
|
|
RealD Tolerance;
|
|
Integer MaxIterations;
|
|
virtual void operator() (LinearOperatorBase<Field> *Linop,const Field &in, Field &out) = 0;
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////
|
|
// Grand daddy iterative method
|
|
/////////////////////////////////////////////////////////////
|
|
template<class Field> class ConjugateGradient : public IterativeProcess<Field> {
|
|
public:
|
|
virtual void operator() (HermitianOperatorBase<Field> *Linop,const Field &in, Field &out) = 0;
|
|
};
|
|
|
|
/////////////////////////////////////////////////////////////
|
|
// A little more modern
|
|
/////////////////////////////////////////////////////////////
|
|
template<class Field> class PreconditionedConjugateGradient : public IterativeProcess<Field> {
|
|
public:
|
|
void operator() (HermitianOperatorBase<Field> *Linop,
|
|
OperatorFunction<Field> *Preconditioner,
|
|
const Field &in,
|
|
Field &out) = 0;
|
|
};
|
|
|
|
// FIXME : To think about
|
|
|
|
// Chroma functionality list defining LinearOperator
|
|
/*
|
|
virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const = 0;
|
|
virtual void operator() (T& chi, const T& psi, enum PlusMinus isign, Real epsilon) const
|
|
virtual const Subset& subset() const = 0;
|
|
virtual unsigned long nFlops() const { return 0; }
|
|
virtual void deriv(P& ds_u, const T& chi, const T& psi, enum PlusMinus isign) const
|
|
class UnprecLinearOperator : public DiffLinearOperator<T,P,Q>
|
|
const Subset& subset() const {return all;}
|
|
};
|
|
*/
|
|
|
|
// Chroma interface defining GaugeAction
|
|
/*
|
|
template<typename P, typename Q> class GaugeAction
|
|
virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
|
|
virtual GaugeState<P,Q>* createState(const Q& q) const
|
|
virtual const GaugeBC<P,Q>& getGaugeBC() const
|
|
virtual const Set& getSet(void) const = 0;
|
|
virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const
|
|
virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
|
|
|
|
class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
|
|
typedef multi1d<LatticeColorMatrix> P;
|
|
typedef multi1d<LatticeColorMatrix> Q;
|
|
virtual void staple(LatticeColorMatrix& result,
|
|
const Handle< GaugeState<P,Q> >& state,
|
|
int mu, int cb) const = 0;
|
|
*/
|
|
|
|
// Chroma interface defining FermionAction
|
|
/*
|
|
template<typename T, typename P, typename Q> class FermAct4D : public FermionAction<T,P,Q>
|
|
virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
|
|
virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
|
|
virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
|
|
virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
|
|
virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
|
|
virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
|
|
virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
|
|
virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
|
|
class DiffFermAct4D : public FermAct4D<T,P,Q>
|
|
virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
|
|
virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
|
|
*/
|
|
}
|
|
|
|
#endif
|