mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Namespace and formatting changes
This commit is contained in:
paboyle
parent
fcf1ccf669
commit
21251f2e1b
@ -1,4 +1,4 @@
|
||||
/*************************************************************************************
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
@ -25,14 +25,14 @@ Author: Christoph Lehner <clehner@bnl.gov>
|
||||
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 */
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_CHEBYSHEV_H
|
||||
#define GRID_CHEBYSHEV_H
|
||||
|
||||
#include <Grid/algorithms/LinearOperator.h>
|
||||
|
||||
namespace Grid {
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
struct ChebyParams : Serializable {
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(ChebyParams,
|
||||
@ -41,337 +41,337 @@ struct ChebyParams : Serializable {
|
||||
int, Npoly);
|
||||
};
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Generic Chebyshev approximations
|
||||
////////////////////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field>
|
||||
class Chebyshev : public OperatorFunction<Field> {
|
||||
private:
|
||||
std::vector<RealD> Coeffs;
|
||||
int order;
|
||||
RealD hi;
|
||||
RealD lo;
|
||||
////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Generic Chebyshev approximations
|
||||
////////////////////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field>
|
||||
class Chebyshev : public OperatorFunction<Field> {
|
||||
private:
|
||||
std::vector<RealD> Coeffs;
|
||||
int order;
|
||||
RealD hi;
|
||||
RealD lo;
|
||||
|
||||
public:
|
||||
void csv(std::ostream &out){
|
||||
RealD diff = hi-lo;
|
||||
RealD delta = (hi-lo)*1.0e-9;
|
||||
for (RealD x=lo; x<hi; x+=delta) {
|
||||
delta*=1.1;
|
||||
RealD f = approx(x);
|
||||
out<< x<<" "<<f<<std::endl;
|
||||
}
|
||||
return;
|
||||
public:
|
||||
void csv(std::ostream &out){
|
||||
RealD diff = hi-lo;
|
||||
RealD delta = (hi-lo)*1.0e-9;
|
||||
for (RealD x=lo; x<hi; x+=delta) {
|
||||
delta*=1.1;
|
||||
RealD f = approx(x);
|
||||
out<< x<<" "<<f<<std::endl;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
// Convenience for plotting the approximation
|
||||
void PlotApprox(std::ostream &out) {
|
||||
out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
|
||||
for(RealD x=lo;x<hi;x+=(hi-lo)/50.0){
|
||||
out <<x<<"\t"<<approx(x)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
Chebyshev(){};
|
||||
Chebyshev(ChebyParams p){ Init(p.alpha,p.beta,p.Npoly);};
|
||||
Chebyshev(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD) ) {Init(_lo,_hi,_order,func);};
|
||||
Chebyshev(RealD _lo,RealD _hi,int _order) {Init(_lo,_hi,_order);};
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// c.f. numerical recipes "chebft"/"chebev". This is sec 5.8 "Chebyshev approximation".
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// CJ: the one we need for Lanczos
|
||||
void Init(RealD _lo,RealD _hi,int _order)
|
||||
{
|
||||
lo=_lo;
|
||||
hi=_hi;
|
||||
order=_order;
|
||||
|
||||
if(order < 2) exit(-1);
|
||||
Coeffs.resize(order);
|
||||
Coeffs.assign(0.,order);
|
||||
Coeffs[order-1] = 1.;
|
||||
};
|
||||
|
||||
void Init(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD))
|
||||
{
|
||||
lo=_lo;
|
||||
hi=_hi;
|
||||
order=_order;
|
||||
|
||||
if(order < 2) exit(-1);
|
||||
Coeffs.resize(order);
|
||||
for(int j=0;j<order;j++){
|
||||
RealD s=0;
|
||||
for(int k=0;k<order;k++){
|
||||
RealD y=std::cos(M_PI*(k+0.5)/order);
|
||||
RealD x=0.5*(y*(hi-lo)+(hi+lo));
|
||||
RealD f=func(x);
|
||||
s=s+f*std::cos( j*M_PI*(k+0.5)/order );
|
||||
}
|
||||
Coeffs[j] = s * 2.0/order;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
void JacksonSmooth(void){
|
||||
RealD M=order;
|
||||
RealD alpha = M_PI/(M+2);
|
||||
RealD lmax = std::cos(alpha);
|
||||
RealD sumUsq =0;
|
||||
std::vector<RealD> U(M);
|
||||
std::vector<RealD> a(M);
|
||||
std::vector<RealD> g(M);
|
||||
for(int n=0;n<=M;n++){
|
||||
U[n] = std::sin((n+1)*std::acos(lmax))/std::sin(std::acos(lmax));
|
||||
sumUsq += U[n]*U[n];
|
||||
}
|
||||
sumUsq = std::sqrt(sumUsq);
|
||||
|
||||
for(int i=1;i<=M;i++){
|
||||
a[i] = U[i]/sumUsq;
|
||||
}
|
||||
g[0] = 1.0;
|
||||
for(int m=1;m<=M;m++){
|
||||
g[m] = 0;
|
||||
for(int i=0;i<=M-m;i++){
|
||||
g[m]+= a[i]*a[m+i];
|
||||
}
|
||||
}
|
||||
for(int m=1;m<=M;m++){
|
||||
Coeffs[m]*=g[m];
|
||||
}
|
||||
}
|
||||
RealD approx(RealD x) // Convenience for plotting the approximation
|
||||
{
|
||||
RealD Tn;
|
||||
RealD Tnm;
|
||||
RealD Tnp;
|
||||
|
||||
RealD y=( x-0.5*(hi+lo))/(0.5*(hi-lo));
|
||||
|
||||
RealD T0=1;
|
||||
RealD T1=y;
|
||||
|
||||
RealD sum;
|
||||
sum = 0.5*Coeffs[0]*T0;
|
||||
sum+= Coeffs[1]*T1;
|
||||
|
||||
Tn =T1;
|
||||
Tnm=T0;
|
||||
for(int i=2;i<order;i++){
|
||||
Tnp=2*y*Tn-Tnm;
|
||||
Tnm=Tn;
|
||||
Tn =Tnp;
|
||||
sum+= Tn*Coeffs[i];
|
||||
}
|
||||
return sum;
|
||||
};
|
||||
|
||||
RealD approxD(RealD x)
|
||||
{
|
||||
RealD Un;
|
||||
RealD Unm;
|
||||
RealD Unp;
|
||||
|
||||
RealD y=( x-0.5*(hi+lo))/(0.5*(hi-lo));
|
||||
|
||||
RealD U0=1;
|
||||
RealD U1=2*y;
|
||||
|
||||
RealD sum;
|
||||
sum = Coeffs[1]*U0;
|
||||
sum+= Coeffs[2]*U1*2.0;
|
||||
|
||||
Un =U1;
|
||||
Unm=U0;
|
||||
for(int i=2;i<order-1;i++){
|
||||
Unp=2*y*Un-Unm;
|
||||
Unm=Un;
|
||||
Un =Unp;
|
||||
sum+= Un*Coeffs[i+1]*(i+1.0);
|
||||
}
|
||||
return sum/(0.5*(hi-lo));
|
||||
};
|
||||
|
||||
RealD approxInv(RealD z, RealD x0, int maxiter, RealD resid) {
|
||||
RealD x = x0;
|
||||
RealD eps;
|
||||
|
||||
int i;
|
||||
for (i=0;i<maxiter;i++) {
|
||||
eps = approx(x) - z;
|
||||
if (fabs(eps / z) < resid)
|
||||
return x;
|
||||
x = x - eps / approxD(x);
|
||||
}
|
||||
|
||||
return std::numeric_limits<double>::quiet_NaN();
|
||||
}
|
||||
|
||||
// Implement the required interface
|
||||
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
|
||||
|
||||
GridBase *grid=in._grid;
|
||||
|
||||
// std::cout << "Chevyshef(): in._grid="<<in._grid<<std::endl;
|
||||
//std::cout <<" Linop.Grid()="<<Linop.Grid()<<"Linop.RedBlackGrid()="<<Linop.RedBlackGrid()<<std::endl;
|
||||
|
||||
int vol=grid->gSites();
|
||||
|
||||
Field T0(grid); T0 = in;
|
||||
Field T1(grid);
|
||||
Field T2(grid);
|
||||
Field y(grid);
|
||||
|
||||
Field *Tnm = &T0;
|
||||
Field *Tn = &T1;
|
||||
Field *Tnp = &T2;
|
||||
|
||||
// Tn=T1 = (xscale M + mscale)in
|
||||
RealD xscale = 2.0/(hi-lo);
|
||||
RealD mscale = -(hi+lo)/(hi-lo);
|
||||
Linop.HermOp(T0,y);
|
||||
T1=y*xscale+in*mscale;
|
||||
|
||||
// sum = .5 c[0] T0 + c[1] T1
|
||||
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
|
||||
for(int n=2;n<order;n++){
|
||||
|
||||
Linop.HermOp(*Tn,y);
|
||||
|
||||
y=xscale*y+mscale*(*Tn);
|
||||
|
||||
*Tnp=2.0*y-(*Tnm);
|
||||
|
||||
out=out+Coeffs[n]* (*Tnp);
|
||||
|
||||
// Cycle pointers to avoid copies
|
||||
Field *swizzle = Tnm;
|
||||
Tnm =Tn;
|
||||
Tn =Tnp;
|
||||
Tnp =swizzle;
|
||||
|
||||
}
|
||||
// Convenience for plotting the approximation
|
||||
void PlotApprox(std::ostream &out) {
|
||||
out<<"Polynomial approx ["<<lo<<","<<hi<<"]"<<std::endl;
|
||||
for(RealD x=lo;x<hi;x+=(hi-lo)/50.0){
|
||||
out <<x<<"\t"<<approx(x)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
Chebyshev(){};
|
||||
Chebyshev(ChebyParams p){ Init(p.alpha,p.beta,p.Npoly);};
|
||||
Chebyshev(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD) ) {Init(_lo,_hi,_order,func);};
|
||||
Chebyshev(RealD _lo,RealD _hi,int _order) {Init(_lo,_hi,_order);};
|
||||
|
||||
template<class Field>
|
||||
class ChebyshevLanczos : public Chebyshev<Field> {
|
||||
private:
|
||||
std::vector<RealD> Coeffs;
|
||||
int order;
|
||||
RealD alpha;
|
||||
RealD beta;
|
||||
RealD mu;
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// c.f. numerical recipes "chebft"/"chebev". This is sec 5.8 "Chebyshev approximation".
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// CJ: the one we need for Lanczos
|
||||
void Init(RealD _lo,RealD _hi,int _order)
|
||||
{
|
||||
lo=_lo;
|
||||
hi=_hi;
|
||||
order=_order;
|
||||
|
||||
if(order < 2) exit(-1);
|
||||
Coeffs.resize(order);
|
||||
Coeffs.assign(0.,order);
|
||||
Coeffs[order-1] = 1.;
|
||||
};
|
||||
|
||||
public:
|
||||
ChebyshevLanczos(RealD _alpha,RealD _beta,RealD _mu,int _order) :
|
||||
void Init(RealD _lo,RealD _hi,int _order, RealD (* func)(RealD))
|
||||
{
|
||||
lo=_lo;
|
||||
hi=_hi;
|
||||
order=_order;
|
||||
|
||||
if(order < 2) exit(-1);
|
||||
Coeffs.resize(order);
|
||||
for(int j=0;j<order;j++){
|
||||
RealD s=0;
|
||||
for(int k=0;k<order;k++){
|
||||
RealD y=std::cos(M_PI*(k+0.5)/order);
|
||||
RealD x=0.5*(y*(hi-lo)+(hi+lo));
|
||||
RealD f=func(x);
|
||||
s=s+f*std::cos( j*M_PI*(k+0.5)/order );
|
||||
}
|
||||
Coeffs[j] = s * 2.0/order;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
void JacksonSmooth(void){
|
||||
RealD M=order;
|
||||
RealD alpha = M_PI/(M+2);
|
||||
RealD lmax = std::cos(alpha);
|
||||
RealD sumUsq =0;
|
||||
std::vector<RealD> U(M);
|
||||
std::vector<RealD> a(M);
|
||||
std::vector<RealD> g(M);
|
||||
for(int n=0;n<=M;n++){
|
||||
U[n] = std::sin((n+1)*std::acos(lmax))/std::sin(std::acos(lmax));
|
||||
sumUsq += U[n]*U[n];
|
||||
}
|
||||
sumUsq = std::sqrt(sumUsq);
|
||||
|
||||
for(int i=1;i<=M;i++){
|
||||
a[i] = U[i]/sumUsq;
|
||||
}
|
||||
g[0] = 1.0;
|
||||
for(int m=1;m<=M;m++){
|
||||
g[m] = 0;
|
||||
for(int i=0;i<=M-m;i++){
|
||||
g[m]+= a[i]*a[m+i];
|
||||
}
|
||||
}
|
||||
for(int m=1;m<=M;m++){
|
||||
Coeffs[m]*=g[m];
|
||||
}
|
||||
}
|
||||
RealD approx(RealD x) // Convenience for plotting the approximation
|
||||
{
|
||||
RealD Tn;
|
||||
RealD Tnm;
|
||||
RealD Tnp;
|
||||
|
||||
RealD y=( x-0.5*(hi+lo))/(0.5*(hi-lo));
|
||||
|
||||
RealD T0=1;
|
||||
RealD T1=y;
|
||||
|
||||
RealD sum;
|
||||
sum = 0.5*Coeffs[0]*T0;
|
||||
sum+= Coeffs[1]*T1;
|
||||
|
||||
Tn =T1;
|
||||
Tnm=T0;
|
||||
for(int i=2;i<order;i++){
|
||||
Tnp=2*y*Tn-Tnm;
|
||||
Tnm=Tn;
|
||||
Tn =Tnp;
|
||||
sum+= Tn*Coeffs[i];
|
||||
}
|
||||
return sum;
|
||||
};
|
||||
|
||||
RealD approxD(RealD x)
|
||||
{
|
||||
RealD Un;
|
||||
RealD Unm;
|
||||
RealD Unp;
|
||||
|
||||
RealD y=( x-0.5*(hi+lo))/(0.5*(hi-lo));
|
||||
|
||||
RealD U0=1;
|
||||
RealD U1=2*y;
|
||||
|
||||
RealD sum;
|
||||
sum = Coeffs[1]*U0;
|
||||
sum+= Coeffs[2]*U1*2.0;
|
||||
|
||||
Un =U1;
|
||||
Unm=U0;
|
||||
for(int i=2;i<order-1;i++){
|
||||
Unp=2*y*Un-Unm;
|
||||
Unm=Un;
|
||||
Un =Unp;
|
||||
sum+= Un*Coeffs[i+1]*(i+1.0);
|
||||
}
|
||||
return sum/(0.5*(hi-lo));
|
||||
};
|
||||
|
||||
RealD approxInv(RealD z, RealD x0, int maxiter, RealD resid) {
|
||||
RealD x = x0;
|
||||
RealD eps;
|
||||
|
||||
int i;
|
||||
for (i=0;i<maxiter;i++) {
|
||||
eps = approx(x) - z;
|
||||
if (fabs(eps / z) < resid)
|
||||
return x;
|
||||
x = x - eps / approxD(x);
|
||||
}
|
||||
|
||||
return std::numeric_limits<double>::quiet_NaN();
|
||||
}
|
||||
|
||||
// Implement the required interface
|
||||
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
|
||||
|
||||
GridBase *grid=in._grid;
|
||||
|
||||
// std::cout << "Chevyshef(): in._grid="<<in._grid<<std::endl;
|
||||
//std::cout <<" Linop.Grid()="<<Linop.Grid()<<"Linop.RedBlackGrid()="<<Linop.RedBlackGrid()<<std::endl;
|
||||
|
||||
int vol=grid->gSites();
|
||||
|
||||
Field T0(grid); T0 = in;
|
||||
Field T1(grid);
|
||||
Field T2(grid);
|
||||
Field y(grid);
|
||||
|
||||
Field *Tnm = &T0;
|
||||
Field *Tn = &T1;
|
||||
Field *Tnp = &T2;
|
||||
|
||||
// Tn=T1 = (xscale M + mscale)in
|
||||
RealD xscale = 2.0/(hi-lo);
|
||||
RealD mscale = -(hi+lo)/(hi-lo);
|
||||
Linop.HermOp(T0,y);
|
||||
T1=y*xscale+in*mscale;
|
||||
|
||||
// sum = .5 c[0] T0 + c[1] T1
|
||||
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
|
||||
for(int n=2;n<order;n++){
|
||||
|
||||
Linop.HermOp(*Tn,y);
|
||||
|
||||
y=xscale*y+mscale*(*Tn);
|
||||
|
||||
*Tnp=2.0*y-(*Tnm);
|
||||
|
||||
out=out+Coeffs[n]* (*Tnp);
|
||||
|
||||
// Cycle pointers to avoid copies
|
||||
Field *swizzle = Tnm;
|
||||
Tnm =Tn;
|
||||
Tn =Tnp;
|
||||
Tnp =swizzle;
|
||||
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template<class Field>
|
||||
class ChebyshevLanczos : public Chebyshev<Field> {
|
||||
private:
|
||||
std::vector<RealD> Coeffs;
|
||||
int order;
|
||||
RealD alpha;
|
||||
RealD beta;
|
||||
RealD mu;
|
||||
|
||||
public:
|
||||
ChebyshevLanczos(RealD _alpha,RealD _beta,RealD _mu,int _order) :
|
||||
alpha(_alpha),
|
||||
beta(_beta),
|
||||
mu(_mu)
|
||||
{
|
||||
order=_order;
|
||||
Coeffs.resize(order);
|
||||
for(int i=0;i<_order;i++){
|
||||
Coeffs[i] = 0.0;
|
||||
}
|
||||
Coeffs[order-1]=1.0;
|
||||
};
|
||||
|
||||
void csv(std::ostream &out){
|
||||
for (RealD x=-1.2*alpha; x<1.2*alpha; x+=(2.0*alpha)/10000) {
|
||||
RealD f = approx(x);
|
||||
out<< x<<" "<<f<<std::endl;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
RealD approx(RealD xx) // Convenience for plotting the approximation
|
||||
{
|
||||
RealD Tn;
|
||||
RealD Tnm;
|
||||
RealD Tnp;
|
||||
Real aa = alpha * alpha;
|
||||
Real bb = beta * beta;
|
||||
|
||||
RealD x = ( 2.0 * (xx-mu)*(xx-mu) - (aa+bb) ) / (aa-bb);
|
||||
|
||||
RealD y= x;
|
||||
|
||||
RealD T0=1;
|
||||
RealD T1=y;
|
||||
|
||||
RealD sum;
|
||||
sum = 0.5*Coeffs[0]*T0;
|
||||
sum+= Coeffs[1]*T1;
|
||||
|
||||
Tn =T1;
|
||||
Tnm=T0;
|
||||
for(int i=2;i<order;i++){
|
||||
Tnp=2*y*Tn-Tnm;
|
||||
Tnm=Tn;
|
||||
Tn =Tnp;
|
||||
sum+= Tn*Coeffs[i];
|
||||
}
|
||||
return sum;
|
||||
};
|
||||
|
||||
// shift_Multiply in Rudy's code
|
||||
void AminusMuSq(LinearOperatorBase<Field> &Linop, const Field &in, Field &out)
|
||||
{
|
||||
GridBase *grid=in._grid;
|
||||
Field tmp(grid);
|
||||
|
||||
RealD aa= alpha*alpha;
|
||||
RealD bb= beta * beta;
|
||||
|
||||
Linop.HermOp(in,out);
|
||||
out = out - mu*in;
|
||||
|
||||
Linop.HermOp(out,tmp);
|
||||
tmp = tmp - mu * out;
|
||||
|
||||
out = (2.0/ (aa-bb) ) * tmp - ((aa+bb)/(aa-bb))*in;
|
||||
};
|
||||
// Implement the required interface
|
||||
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
|
||||
|
||||
GridBase *grid=in._grid;
|
||||
|
||||
int vol=grid->gSites();
|
||||
|
||||
Field T0(grid); T0 = in;
|
||||
Field T1(grid);
|
||||
Field T2(grid);
|
||||
Field y(grid);
|
||||
|
||||
Field *Tnm = &T0;
|
||||
Field *Tn = &T1;
|
||||
Field *Tnp = &T2;
|
||||
|
||||
// Tn=T1 = (xscale M )*in
|
||||
AminusMuSq(Linop,T0,T1);
|
||||
|
||||
// sum = .5 c[0] T0 + c[1] T1
|
||||
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
|
||||
for(int n=2;n<order;n++){
|
||||
|
||||
AminusMuSq(Linop,*Tn,y);
|
||||
|
||||
*Tnp=2.0*y-(*Tnm);
|
||||
|
||||
out=out+Coeffs[n]* (*Tnp);
|
||||
|
||||
// Cycle pointers to avoid copies
|
||||
Field *swizzle = Tnm;
|
||||
Tnm =Tn;
|
||||
Tn =Tnp;
|
||||
Tnp =swizzle;
|
||||
|
||||
}
|
||||
beta(_beta),
|
||||
mu(_mu)
|
||||
{
|
||||
order=_order;
|
||||
Coeffs.resize(order);
|
||||
for(int i=0;i<_order;i++){
|
||||
Coeffs[i] = 0.0;
|
||||
}
|
||||
Coeffs[order-1]=1.0;
|
||||
};
|
||||
}
|
||||
|
||||
void csv(std::ostream &out){
|
||||
for (RealD x=-1.2*alpha; x<1.2*alpha; x+=(2.0*alpha)/10000) {
|
||||
RealD f = approx(x);
|
||||
out<< x<<" "<<f<<std::endl;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
RealD approx(RealD xx) // Convenience for plotting the approximation
|
||||
{
|
||||
RealD Tn;
|
||||
RealD Tnm;
|
||||
RealD Tnp;
|
||||
Real aa = alpha * alpha;
|
||||
Real bb = beta * beta;
|
||||
|
||||
RealD x = ( 2.0 * (xx-mu)*(xx-mu) - (aa+bb) ) / (aa-bb);
|
||||
|
||||
RealD y= x;
|
||||
|
||||
RealD T0=1;
|
||||
RealD T1=y;
|
||||
|
||||
RealD sum;
|
||||
sum = 0.5*Coeffs[0]*T0;
|
||||
sum+= Coeffs[1]*T1;
|
||||
|
||||
Tn =T1;
|
||||
Tnm=T0;
|
||||
for(int i=2;i<order;i++){
|
||||
Tnp=2*y*Tn-Tnm;
|
||||
Tnm=Tn;
|
||||
Tn =Tnp;
|
||||
sum+= Tn*Coeffs[i];
|
||||
}
|
||||
return sum;
|
||||
};
|
||||
|
||||
// shift_Multiply in Rudy's code
|
||||
void AminusMuSq(LinearOperatorBase<Field> &Linop, const Field &in, Field &out)
|
||||
{
|
||||
GridBase *grid=in._grid;
|
||||
Field tmp(grid);
|
||||
|
||||
RealD aa= alpha*alpha;
|
||||
RealD bb= beta * beta;
|
||||
|
||||
Linop.HermOp(in,out);
|
||||
out = out - mu*in;
|
||||
|
||||
Linop.HermOp(out,tmp);
|
||||
tmp = tmp - mu * out;
|
||||
|
||||
out = (2.0/ (aa-bb) ) * tmp - ((aa+bb)/(aa-bb))*in;
|
||||
};
|
||||
// Implement the required interface
|
||||
void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
|
||||
|
||||
GridBase *grid=in._grid;
|
||||
|
||||
int vol=grid->gSites();
|
||||
|
||||
Field T0(grid); T0 = in;
|
||||
Field T1(grid);
|
||||
Field T2(grid);
|
||||
Field y(grid);
|
||||
|
||||
Field *Tnm = &T0;
|
||||
Field *Tn = &T1;
|
||||
Field *Tnp = &T2;
|
||||
|
||||
// Tn=T1 = (xscale M )*in
|
||||
AminusMuSq(Linop,T0,T1);
|
||||
|
||||
// sum = .5 c[0] T0 + c[1] T1
|
||||
out = (0.5*Coeffs[0])*T0 + Coeffs[1]*T1;
|
||||
for(int n=2;n<order;n++){
|
||||
|
||||
AminusMuSq(Linop,*Tn,y);
|
||||
|
||||
*Tnp=2.0*y-(*Tnm);
|
||||
|
||||
out=out+Coeffs[n]* (*Tnp);
|
||||
|
||||
// Cycle pointers to avoid copies
|
||||
Field *swizzle = Tnm;
|
||||
Tnm =Tn;
|
||||
Tn =Tnp;
|
||||
Tnp =swizzle;
|
||||
|
||||
}
|
||||
}
|
||||
};
|
||||
NAMESPACE_END(Grid);
|
||||
#endif
|
||||
|
||||
@ -26,127 +26,127 @@ with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
/* END LEGAL */
|
||||
|
||||
#ifndef INCLUDED_FORECAST_H
|
||||
#define INCLUDED_FORECAST_H
|
||||
|
||||
namespace Grid {
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
// Abstract base class.
|
||||
// Takes a matrix (Mat), a source (phi), and a vector of Fields (chi)
|
||||
// and returns a forecasted solution to the system D*psi = phi (psi).
|
||||
template<class Matrix, class Field>
|
||||
class Forecast
|
||||
// Abstract base class.
|
||||
// Takes a matrix (Mat), a source (phi), and a vector of Fields (chi)
|
||||
// and returns a forecasted solution to the system D*psi = phi (psi).
|
||||
template<class Matrix, class Field>
|
||||
class Forecast
|
||||
{
|
||||
public:
|
||||
virtual Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& chi) = 0;
|
||||
};
|
||||
|
||||
// Implementation of Brower et al.'s chronological inverter (arXiv:hep-lat/9509012),
|
||||
// used to forecast solutions across poles of the EOFA heatbath.
|
||||
//
|
||||
// Modified from CPS (cps_pp/src/util/dirac_op/d_op_base/comsrc/minresext.C)
|
||||
template<class Matrix, class Field>
|
||||
class ChronoForecast : public Forecast<Matrix,Field>
|
||||
{
|
||||
public:
|
||||
Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& prev_solns)
|
||||
{
|
||||
public:
|
||||
virtual Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& chi) = 0;
|
||||
int degree = prev_solns.size();
|
||||
Field chi(phi); // forecasted solution
|
||||
|
||||
// Trivial cases
|
||||
if(degree == 0){ chi = zero; return chi; }
|
||||
else if(degree == 1){ return prev_solns[0]; }
|
||||
|
||||
RealD dot;
|
||||
ComplexD xp;
|
||||
Field r(phi); // residual
|
||||
Field Mv(phi);
|
||||
std::vector<Field> v(prev_solns); // orthonormalized previous solutions
|
||||
std::vector<Field> MdagMv(degree,phi);
|
||||
|
||||
// Array to hold the matrix elements
|
||||
std::vector<std::vector<ComplexD>> G(degree, std::vector<ComplexD>(degree));
|
||||
|
||||
// Solution and source vectors
|
||||
std::vector<ComplexD> a(degree);
|
||||
std::vector<ComplexD> b(degree);
|
||||
|
||||
// Orthonormalize the vector basis
|
||||
for(int i=0; i<degree; i++){
|
||||
v[i] *= 1.0/std::sqrt(norm2(v[i]));
|
||||
for(int j=i+1; j<degree; j++){ v[j] -= innerProduct(v[i],v[j]) * v[i]; }
|
||||
}
|
||||
|
||||
// Perform sparse matrix multiplication and construct rhs
|
||||
for(int i=0; i<degree; i++){
|
||||
b[i] = innerProduct(v[i],phi);
|
||||
Mat.M(v[i],Mv);
|
||||
Mat.Mdag(Mv,MdagMv[i]);
|
||||
G[i][i] = innerProduct(v[i],MdagMv[i]);
|
||||
}
|
||||
|
||||
// Construct the matrix
|
||||
for(int j=0; j<degree; j++){
|
||||
for(int k=j+1; k<degree; k++){
|
||||
G[j][k] = innerProduct(v[j],MdagMv[k]);
|
||||
G[k][j] = std::conj(G[j][k]);
|
||||
}}
|
||||
|
||||
// Gauss-Jordan elimination with partial pivoting
|
||||
for(int i=0; i<degree; i++){
|
||||
|
||||
// Perform partial pivoting
|
||||
int k = i;
|
||||
for(int j=i+1; j<degree; j++){ if(std::abs(G[j][j]) > std::abs(G[k][k])){ k = j; } }
|
||||
if(k != i){
|
||||
xp = b[k];
|
||||
b[k] = b[i];
|
||||
b[i] = xp;
|
||||
for(int j=0; j<degree; j++){
|
||||
xp = G[k][j];
|
||||
G[k][j] = G[i][j];
|
||||
G[i][j] = xp;
|
||||
}
|
||||
}
|
||||
|
||||
// Convert matrix to upper triangular form
|
||||
for(int j=i+1; j<degree; j++){
|
||||
xp = G[j][i]/G[i][i];
|
||||
b[j] -= xp * b[i];
|
||||
for(int k=0; k<degree; k++){ G[j][k] -= xp*G[i][k]; }
|
||||
}
|
||||
}
|
||||
|
||||
// Use Gaussian elimination to solve equations and calculate initial guess
|
||||
chi = zero;
|
||||
r = phi;
|
||||
for(int i=degree-1; i>=0; i--){
|
||||
a[i] = 0.0;
|
||||
for(int j=i+1; j<degree; j++){ a[i] += G[i][j] * a[j]; }
|
||||
a[i] = (b[i]-a[i])/G[i][i];
|
||||
chi += a[i]*v[i];
|
||||
r -= a[i]*MdagMv[i];
|
||||
}
|
||||
|
||||
RealD true_r(0.0);
|
||||
ComplexD tmp;
|
||||
for(int i=0; i<degree; i++){
|
||||
tmp = -b[i];
|
||||
for(int j=0; j<degree; j++){ tmp += G[i][j]*a[j]; }
|
||||
tmp = std::conj(tmp)*tmp;
|
||||
true_r += std::sqrt(tmp.real());
|
||||
}
|
||||
|
||||
RealD error = std::sqrt(norm2(r)/norm2(phi));
|
||||
std::cout << GridLogMessage << "ChronoForecast: |res|/|src| = " << error << std::endl;
|
||||
|
||||
return chi;
|
||||
};
|
||||
};
|
||||
|
||||
// Implementation of Brower et al.'s chronological inverter (arXiv:hep-lat/9509012),
|
||||
// used to forecast solutions across poles of the EOFA heatbath.
|
||||
//
|
||||
// Modified from CPS (cps_pp/src/util/dirac_op/d_op_base/comsrc/minresext.C)
|
||||
template<class Matrix, class Field>
|
||||
class ChronoForecast : public Forecast<Matrix,Field>
|
||||
{
|
||||
public:
|
||||
Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& prev_solns)
|
||||
{
|
||||
int degree = prev_solns.size();
|
||||
Field chi(phi); // forecasted solution
|
||||
|
||||
// Trivial cases
|
||||
if(degree == 0){ chi = zero; return chi; }
|
||||
else if(degree == 1){ return prev_solns[0]; }
|
||||
|
||||
RealD dot;
|
||||
ComplexD xp;
|
||||
Field r(phi); // residual
|
||||
Field Mv(phi);
|
||||
std::vector<Field> v(prev_solns); // orthonormalized previous solutions
|
||||
std::vector<Field> MdagMv(degree,phi);
|
||||
|
||||
// Array to hold the matrix elements
|
||||
std::vector<std::vector<ComplexD>> G(degree, std::vector<ComplexD>(degree));
|
||||
|
||||
// Solution and source vectors
|
||||
std::vector<ComplexD> a(degree);
|
||||
std::vector<ComplexD> b(degree);
|
||||
|
||||
// Orthonormalize the vector basis
|
||||
for(int i=0; i<degree; i++){
|
||||
v[i] *= 1.0/std::sqrt(norm2(v[i]));
|
||||
for(int j=i+1; j<degree; j++){ v[j] -= innerProduct(v[i],v[j]) * v[i]; }
|
||||
}
|
||||
|
||||
// Perform sparse matrix multiplication and construct rhs
|
||||
for(int i=0; i<degree; i++){
|
||||
b[i] = innerProduct(v[i],phi);
|
||||
Mat.M(v[i],Mv);
|
||||
Mat.Mdag(Mv,MdagMv[i]);
|
||||
G[i][i] = innerProduct(v[i],MdagMv[i]);
|
||||
}
|
||||
|
||||
// Construct the matrix
|
||||
for(int j=0; j<degree; j++){
|
||||
for(int k=j+1; k<degree; k++){
|
||||
G[j][k] = innerProduct(v[j],MdagMv[k]);
|
||||
G[k][j] = std::conj(G[j][k]);
|
||||
}}
|
||||
|
||||
// Gauss-Jordan elimination with partial pivoting
|
||||
for(int i=0; i<degree; i++){
|
||||
|
||||
// Perform partial pivoting
|
||||
int k = i;
|
||||
for(int j=i+1; j<degree; j++){ if(std::abs(G[j][j]) > std::abs(G[k][k])){ k = j; } }
|
||||
if(k != i){
|
||||
xp = b[k];
|
||||
b[k] = b[i];
|
||||
b[i] = xp;
|
||||
for(int j=0; j<degree; j++){
|
||||
xp = G[k][j];
|
||||
G[k][j] = G[i][j];
|
||||
G[i][j] = xp;
|
||||
}
|
||||
}
|
||||
|
||||
// Convert matrix to upper triangular form
|
||||
for(int j=i+1; j<degree; j++){
|
||||
xp = G[j][i]/G[i][i];
|
||||
b[j] -= xp * b[i];
|
||||
for(int k=0; k<degree; k++){ G[j][k] -= xp*G[i][k]; }
|
||||
}
|
||||
}
|
||||
|
||||
// Use Gaussian elimination to solve equations and calculate initial guess
|
||||
chi = zero;
|
||||
r = phi;
|
||||
for(int i=degree-1; i>=0; i--){
|
||||
a[i] = 0.0;
|
||||
for(int j=i+1; j<degree; j++){ a[i] += G[i][j] * a[j]; }
|
||||
a[i] = (b[i]-a[i])/G[i][i];
|
||||
chi += a[i]*v[i];
|
||||
r -= a[i]*MdagMv[i];
|
||||
}
|
||||
|
||||
RealD true_r(0.0);
|
||||
ComplexD tmp;
|
||||
for(int i=0; i<degree; i++){
|
||||
tmp = -b[i];
|
||||
for(int j=0; j<degree; j++){ tmp += G[i][j]*a[j]; }
|
||||
tmp = std::conj(tmp)*tmp;
|
||||
true_r += std::sqrt(tmp.real());
|
||||
}
|
||||
|
||||
RealD error = std::sqrt(norm2(r)/norm2(phi));
|
||||
std::cout << GridLogMessage << "ChronoForecast: |res|/|src| = " << error << std::endl;
|
||||
|
||||
return chi;
|
||||
};
|
||||
};
|
||||
|
||||
}
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
#endif
|
||||
|
||||
@ -27,7 +27,8 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
/* END LEGAL */
|
||||
#include <Grid/GridCore.h>
|
||||
|
||||
namespace Grid {
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
double MultiShiftFunction::approx(double x)
|
||||
{
|
||||
double a = norm;
|
||||
@ -53,4 +54,4 @@ void MultiShiftFunction::csv(std::ostream &out)
|
||||
}
|
||||
return;
|
||||
}
|
||||
}
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
@ -29,7 +29,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
#ifndef MULTI_SHIFT_FUNCTION
|
||||
#define MULTI_SHIFT_FUNCTION
|
||||
|
||||
namespace Grid {
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
class MultiShiftFunction {
|
||||
public:
|
||||
@ -63,5 +63,5 @@ public:
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
NAMESPACE_END(Grid);
|
||||
#endif
|
||||
|
||||
@ -58,8 +58,8 @@
|
||||
|
||||
/* Compute the partial fraction expansion coefficients (alpha) from the
|
||||
* factored form */
|
||||
namespace Grid {
|
||||
namespace Approx {
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
NAMESPACE_BEGIN(Approx);
|
||||
|
||||
static void construct_partfrac(izd *z) {
|
||||
int dn = z -> dn, dd = z -> dd, type = z -> type;
|
||||
@ -723,5 +723,6 @@ int main(int argc, char** argv) {
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
|
||||
#endif /* TEST */
|
||||
NAMESPACE_END(Approx);
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
@ -1,13 +1,13 @@
|
||||
/* -*- Mode: C; comment-column: 22; fill-column: 79; -*- */
|
||||
|
||||
#ifdef __cplusplus
|
||||
namespace Grid {
|
||||
namespace Approx {
|
||||
#include <Grid/Namespace.h>
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
NAMESPACE_BEGIN(Approx);
|
||||
#endif
|
||||
|
||||
#define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY>
|
||||
|
||||
|
||||
#ifndef ZOLOTAREV_INTERNAL
|
||||
#ifndef PRECISION
|
||||
#define PRECISION double
|
||||
@ -83,5 +83,6 @@ void zolotarev_free(zolotarev_data *zdata);
|
||||
#endif
|
||||
|
||||
#ifdef __cplusplus
|
||||
}}
|
||||
NAMESPACE_END(Approx);
|
||||
NAMESPACE_END(Grid);
|
||||
#endif
|
||||
|
||||
Loading…
Reference in New Issue
Block a user