mirror of
https://github.com/paboyle/Grid.git
synced 2026-10-27 16:13:19 +00:00
Compare commits
3
Commits
d7c9076228
...
5f8f396e9b
| Author | SHA1 | Date | |
|---|---|---|---|
 Peter Boyle
|
5f8f396e9b | ||
|
Peter Boyle
|
841e59e8c0 | ||
|
Peter Boyle
|
b55ca3d50d |
No files matched your search
@@ -0,0 +1,59 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./Grid/algorithms/iterative/GCRCoefficients.h
|
||||
|
||||
Copyright (C) 2026
|
||||
|
||||
Author: Peter Boyle <pboyle@bnl.gov>
|
||||
|
||||
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.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#pragma once
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// Recorded GCR coefficients: per-step means over calls of the step length
|
||||
// a_k and the orthogonalisation coefficients b_kj (already scaled and
|
||||
// signed as applied: p_{k+1} = r + sum_j b_kj p_{k-j}).
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
struct GCRCoefficients {
|
||||
int mmax = 0;
|
||||
std::vector<ComplexD> a_sum; // [k]
|
||||
std::vector<int> a_n;
|
||||
std::vector<std::vector<ComplexD> > b_sum; // [k][j]
|
||||
std::vector<std::vector<int> > b_n;
|
||||
void RecordA(int k, ComplexD a){
|
||||
if ( (int)a_sum.size() <= k ) { a_sum.resize(k+1,ComplexD(0.0)); a_n.resize(k+1,0); }
|
||||
a_sum[k] += a; a_n[k]++;
|
||||
}
|
||||
void RecordB(int k, const std::vector<ComplexD> &b){
|
||||
if ( (int)b_sum.size() <= k ) { b_sum.resize(k+1); b_n.resize(k+1); }
|
||||
if ( b_sum[k].size() < b.size() ) { b_sum[k].resize(b.size(),ComplexD(0.0)); b_n[k].resize(b.size(),0); }
|
||||
for(int j=0;j<(int)b.size();j++){ b_sum[k][j] += b[j]; b_n[k][j]++; }
|
||||
}
|
||||
int Steps(void) const { return a_sum.size(); }
|
||||
int Calls(void) const { return a_n.size() ? a_n[0] : 0; }
|
||||
ComplexD A(int k) const { return a_sum[k]/(double)a_n[k]; }
|
||||
int NB(int k) const { return (k<(int)b_sum.size()) ? b_sum[k].size() : 0; }
|
||||
ComplexD B(int k,int j) const { return b_sum[k][j]/(double)b_n[k][j]; }
|
||||
void Report(const std::string &name) const {
|
||||
std::cout << GridLogMessage << "GCRCoefficients " << name << ": " << Calls() << " calls, " << Steps() << " steps, mmax " << mmax << std::endl;
|
||||
for(int k=0;k<Steps();k++){
|
||||
std::cout << GridLogMessage << " step " << k << " a=(" << real(A(k)) << "," << imag(A(k)) << ")";
|
||||
for(int j=0;j<NB(k);j++) std::cout << " b[" << j << "]=(" << real(B(k,j)) << "," << imag(B(k,j)) << ")";
|
||||
std::cout << std::endl;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
@@ -28,6 +28,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
/* END LEGAL */
|
||||
#ifndef GRID_PREC_GCR_NON_HERM_H
|
||||
#define GRID_PREC_GCR_NON_HERM_H
|
||||
#include <Grid/algorithms/iterative/GCRCoefficients.h>
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
//VPGCR Abe and Zhang, 2005.
|
||||
@@ -79,6 +80,11 @@ public:
|
||||
// same applies and no reductions. Off by default; boss rank prints.
|
||||
int LogCoeffs = 0;
|
||||
void LogCoefficients(int l) { LogCoeffs = l; };
|
||||
// Optional recorder of the per-step coefficients (means over calls), for
|
||||
// replay by GCRReplaySmoother (Smoothers.h). Records only; no effect on
|
||||
// the iteration.
|
||||
GCRCoefficients *Recorder = nullptr;
|
||||
void SetCoefficientRecorder(GCRCoefficients *r) { Recorder = r; if(r) r->mmax = mmax; };
|
||||
|
||||
PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,LinearFunction<Field> &Prec,int _mmax,int _nstep) :
|
||||
Tolerance(tol),
|
||||
@@ -223,6 +229,7 @@ public:
|
||||
LinalgTimer.Start();
|
||||
rq= innerProduct(q[peri_k],r); // what if rAr not real?
|
||||
a = rq/qq[peri_k];
|
||||
if ( Recorder ) Recorder->RecordA(k,a);
|
||||
|
||||
axpy(psi,a,p[peri_k],psi);
|
||||
|
||||
@@ -274,6 +281,7 @@ public:
|
||||
bcoef[back] = -bcoef[back]/qq[peri_back];
|
||||
if ( LogCoeffs ) bs<<" b["<<back<<"]="<<bcoef[back];
|
||||
}
|
||||
if ( Recorder && northog ) Recorder->RecordB(k,bcoef);
|
||||
if ( northog ) {
|
||||
axpyMulti(p[peri_kp],bcoef,pwin);
|
||||
qq[peri_kp]=axpyMultiNorm(q[peri_kp],bcoef,qwin);
|
||||
|
||||
@@ -140,6 +140,7 @@ public:
|
||||
///////////////////////////////////////////////////////////////////////////
|
||||
void Leaf(BlockCyclicMatrix &A, int64_t b)
|
||||
{
|
||||
GRID_TRACE("SchurLeaf");
|
||||
BlockCyclicLayout &L = A.layout;
|
||||
nLeaf++;
|
||||
if ( (int)(b % L.Pr) != L.prow ) return;
|
||||
@@ -211,6 +212,7 @@ public:
|
||||
int64_t span = b1-b0;
|
||||
GRID_ASSERT( span >= 1 );
|
||||
if ( span == 1 ) { Leaf(A, b0); return; }
|
||||
GRID_TRACE("SchurNode");
|
||||
nNode++;
|
||||
|
||||
int64_t bm = b0 + span/2;
|
||||
@@ -244,8 +246,10 @@ public:
|
||||
tGemm += usecond();
|
||||
|
||||
// 9. Off-diagonal signs
|
||||
{ GRID_TRACE("SchurCopy");
|
||||
WindowCopyScale(mone, Ut, A, c0,m, m,c1);
|
||||
WindowCopyScale(mone, Tt, A, m,c1, c0,m);
|
||||
}
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////
|
||||
@@ -285,6 +289,22 @@ public:
|
||||
<< " leaf " << tLeaf/1.0e6
|
||||
<< " copy " << tCopy/1.0e6 << ")"
|
||||
<< std::endl;
|
||||
// SUMMA breakdown: boss-rank seconds, plus the ring/gemm time spread over
|
||||
// ranks (min/max) -- skew shows as max >> min.
|
||||
RealD ring = (SUMMA.tRingA+SUMMA.tRingB)/1.0e6;
|
||||
RealD rmin = -ring, rmax = ring; grid->GlobalMax(rmin); grid->GlobalMax(rmax); rmin = -rmin; // no GlobalMin: max of negation
|
||||
RealD gmin = -SUMMA.tGemm/1.0e6, gmax = SUMMA.tGemm/1.0e6; grid->GlobalMax(gmin); grid->GlobalMax(gmax); gmin = -gmin;
|
||||
double gb = SUMMA.bytesRing/1.0e9;
|
||||
std::cout << GridLogMessage << "BlockCyclicSumma:"
|
||||
<< " multiplies " << SUMMA.nMultiply << " gemms " << SUMMA.nGemm << " ring msgs " << SUMMA.nRingMsg
|
||||
<< " | boss secs: alloc " << SUMMA.tAlloc/1.0e6
|
||||
<< " pack " << SUMMA.tPack/1.0e6
|
||||
<< " ringA " << SUMMA.tRingA/1.0e6 << " ringB " << SUMMA.tRingB/1.0e6
|
||||
<< " gemm " << SUMMA.tGemm/1.0e6
|
||||
<< " | ring min/max over ranks " << rmin << "/" << rmax
|
||||
<< " gemm min/max " << gmin << "/" << gmax
|
||||
<< " | ring bytes/rank " << gb << " GB -> " << (ring>0 ? gb/ring : 0.0) << " GB/s/rank (boss)"
|
||||
<< std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@@ -127,6 +127,15 @@ class BlockCyclicSumma
|
||||
public:
|
||||
GridBLAS BLAS;
|
||||
|
||||
// Per-rank telemetry (boss-rank seconds when printed; no reductions here).
|
||||
// Every Multiply is: buffer alloc, pack panels, ring A along the process
|
||||
// row, ring B along the process column, then the local GEMMs. The rings
|
||||
// are synchronous SendToRecvFrom, so tRing is time the GPU is idle unless
|
||||
// a future version overlaps them with the GEMMs.
|
||||
double tAlloc=0, tPack=0, tRingA=0, tRingB=0, tGemm=0;
|
||||
uint64_t bytesRing=0, nRingMsg=0, nMultiply=0, nGemm=0;
|
||||
void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; }
|
||||
|
||||
static int Overlap(int64_t a0,int64_t a1,int64_t b0,int64_t b1)
|
||||
{ return (a0 < b1) && (b0 < a1); }
|
||||
|
||||
@@ -188,8 +197,11 @@ public:
|
||||
const uint64_t slotA = (uint64_t)mloc_i*nb;
|
||||
const uint64_t slotB1 = (uint64_t)nb*nloc_j; // one panel
|
||||
const uint64_t slotB = (uint64_t)S*slotB1;
|
||||
nMultiply++;
|
||||
tAlloc -= usecond();
|
||||
deviceVector<ComplexD> Abuf( slotA*Pc ? slotA*Pc : 1 );
|
||||
deviceVector<ComplexD> Bbuf( slotB*Pr ? slotB*Pr : 1 );
|
||||
tAlloc += usecond();
|
||||
|
||||
deviceVector<ComplexD *> ap(1), bp(1), cp(1);
|
||||
std::vector<ComplexD *> ptr(1);
|
||||
@@ -201,6 +213,8 @@ public:
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
// Pack MY panels of this round into my origin slots.
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
tPack -= usecond();
|
||||
{ GRID_TRACE("SummaPack");
|
||||
for(int64_t s=r0; s<r1; s++){
|
||||
int64_t nb_s = L.BlockSize(s);
|
||||
if ( (int)(s%Pc) == pcol && mloc_i ){ // my A panel: block-column s
|
||||
@@ -235,12 +249,16 @@ public:
|
||||
}
|
||||
}
|
||||
accelerator_barrier();
|
||||
}
|
||||
tPack += usecond();
|
||||
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
// Ring allgather along my process ROW: Pc-1 symmetric steps. At
|
||||
// step t send the slot of origin (pcol-t+1), receive origin (pcol-t).
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
if ( Pc > 1 && slotA ){
|
||||
GRID_TRACE("SummaRingA");
|
||||
tRingA -= usecond();
|
||||
int dest = prow*Pc + (pcol+1)%Pc;
|
||||
int src = prow*Pc + (pcol-1+Pc)%Pc;
|
||||
for(int t=1;t<Pc;t++){
|
||||
@@ -249,12 +267,16 @@ public:
|
||||
grid->SendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest,
|
||||
(void *)(&Abuf[0]+slotA*cr), src,
|
||||
slotA*sizeof(ComplexD));
|
||||
bytesRing += slotA*sizeof(ComplexD); nRingMsg++;
|
||||
}
|
||||
tRingA += usecond();
|
||||
}
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
// Ring allgather along my process COLUMN: Pr-1 symmetric steps.
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
if ( Pr > 1 && slotB ){
|
||||
GRID_TRACE("SummaRingB");
|
||||
tRingB -= usecond();
|
||||
int dest = ((prow+1)%Pr)*Pc + pcol;
|
||||
int src = ((prow-1+Pr)%Pr)*Pc + pcol;
|
||||
for(int t=1;t<Pr;t++){
|
||||
@@ -263,12 +285,16 @@ public:
|
||||
grid->SendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest,
|
||||
(void *)(&Bbuf[0]+slotB*rr), src,
|
||||
slotB*sizeof(ComplexD));
|
||||
bytesRing += slotB*sizeof(ComplexD); nRingMsg++;
|
||||
}
|
||||
tRingB += usecond();
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
// Local update, ascending s: fixed order, bitwise-reproducible.
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
tGemm -= usecond();
|
||||
{ GRID_TRACE("SummaGEMM");
|
||||
for(int64_t s=r0; s<r1; s++){
|
||||
int64_t nb_s = L.BlockSize(s);
|
||||
if ( !(mloc_i && nloc_j && nb_s) ) { firstblock = 0; continue; }
|
||||
@@ -291,7 +317,10 @@ public:
|
||||
bp, (int)nb_s,
|
||||
beta_use, cp, (int)C.layout.mloc);
|
||||
BLAS.synchronise();
|
||||
nGemm++;
|
||||
}
|
||||
}
|
||||
tGemm += usecond();
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
@@ -34,3 +34,4 @@ Author: Peter Boyle <pboyle@bnl.gov>
|
||||
#include <Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHS.h>
|
||||
#include <Grid/algorithms/multigrid/GeneralCoarsenedMatrixMultiRHSV2.h>
|
||||
#include <Grid/algorithms/multigrid/MrhsPromotedOperator.h>
|
||||
#include <Grid/algorithms/multigrid/Smoothers.h>
|
||||
@@ -0,0 +1,269 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./Grid/algorithms/multigrid/Smoothers.h
|
||||
|
||||
Copyright (C) 2026
|
||||
|
||||
Author: Peter Boyle <pboyle@bnl.gov>
|
||||
|
||||
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.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#pragma once
|
||||
#include <Grid/algorithms/iterative/GCRCoefficients.h>
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// Multigrid smoothers as LinearFunction<Field> objects.
|
||||
//
|
||||
// Gathered here from ~15 copies in tests/solver and tests/debug (HDCR /
|
||||
// HDCG era), plus the fixed-polynomial smoothers of 2026:
|
||||
//
|
||||
// ChebyshevSmoother Chebyshev approx to 1/x on [lo,hi], applied
|
||||
// through HermOp (the original HDCG smoother on
|
||||
// a Hermitian shifted operator).
|
||||
// ChebyshevNonHermitianSmoother same polynomial applied through Op(), for
|
||||
// a non-Hermitian (near-normal, real-spectrum)
|
||||
// smoother operator such as shifted PVdagM.
|
||||
// ChebyshevInverter one Chebyshev-corrected step with residual print.
|
||||
// MirsSmoother shifted-MdagM CG, HDCG arXiv:1402.2585.
|
||||
// GCRReplaySmoother replays a GCR's recorded step lengths a_k and
|
||||
// orthogonalisation coefficients b_kj with NO
|
||||
// inner products: one matvec per step, zero
|
||||
// reductions. The "PreconditionerMirsPoly"
|
||||
// idea of 1402.2585 p.13 applied to GCR.
|
||||
//
|
||||
// Recording: PrecGeneralisedConjugateResidualNonHermitian::SetCoefficientRecorder
|
||||
// accumulates per-step means over calls into a GCRCoefficients; construct a
|
||||
// GCRReplaySmoother from it.
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
inline RealD InverseApproximation(RealD x){ return 1.0/x; }
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// HermOp-based Chebyshev smoother. The second template parameter and the
|
||||
// 5-argument constructor exist only so the historical call sites
|
||||
// ChebyshevSmoother<LatticeFermion,DomainWallFermionD> S(lo,hi,ord,HermOp,Ddwf);
|
||||
// compile unchanged; the Matrix argument was never used.
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field,class Matrix=void> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _SmootherOperator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
template<class M>
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator, M &) :
|
||||
ChebyshevSmoother(_lo,_hi,_ord,SmootherOperator) {};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Cheby(_SmootherOperator,in,out);
|
||||
}
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// Op()-based Chebyshev smoother: x = S(A) r with S the Chebyshev fit to 1/x
|
||||
// on [lo,hi]. Same three-term recurrence as Chebyshev<Field>::operator()
|
||||
// but through Op, for the non-Hermitian smoother operators of the PVdagM
|
||||
// multigrid (real coefficients / near-normal, as the recorded GCR
|
||||
// coefficients show). `order` matvecs, no reductions.
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field> class ChebyshevNonHermitianSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
LinearOperatorBase<Field> &Linop;
|
||||
RealD lo, hi; int order;
|
||||
std::vector<RealD> Coeffs;
|
||||
ChebyshevNonHermitianSmoother(RealD _lo,RealD _hi,int _order,LinearOperatorBase<Field> &Op)
|
||||
: Linop(Op), lo(_lo), hi(_hi), order(_order)
|
||||
{
|
||||
GRID_ASSERT(order>=2);
|
||||
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));
|
||||
s=s+InverseApproximation(x)*std::cos( j*M_PI*(k+0.5)/order );
|
||||
}
|
||||
Coeffs[j] = s * 2.0/order;
|
||||
}
|
||||
std::cout << GridLogMessage<<" ChebyshevNonHermitian smoother order "<<order<<" ["<<lo<<","<<hi<<"]"<<std::endl;
|
||||
}
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
GRID_TRACE("ChebyshevSmoother");
|
||||
GridBase *grid=in.Grid();
|
||||
Field T0(grid); T0 = in;
|
||||
Field T1(grid), T2(grid), y(grid);
|
||||
Field *Tnm=&T0, *Tn=&T1, *Tnp=&T2;
|
||||
RealD xscale = 2.0/(hi-lo);
|
||||
RealD mscale = -(hi+lo)/(hi-lo);
|
||||
Linop.Op(T0,y);
|
||||
axpby(T1,xscale,mscale,y,in);
|
||||
axpby(out,0.5*Coeffs[0],Coeffs[1],T0,T1);
|
||||
for(int n=2;n<order;n++){
|
||||
Linop.Op(*Tn,y);
|
||||
axpby(y,xscale,mscale,y,(*Tn));
|
||||
axpby(*Tnp,2.0,-1.0,y,(*Tnm));
|
||||
if ( Coeffs[n] != 0.0 ) axpy(out,Coeffs[n],*Tnp,out);
|
||||
Field *swizzle=Tnm; Tnm=Tn; Tn=Tnp; Tnp=swizzle;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field> class ChebyshevInverter : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _Operator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) :
|
||||
_Operator(Operator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field r(in.Grid());
|
||||
Field AinvR(in.Grid());
|
||||
_Operator.HermOp(out,r);
|
||||
r = in - r; // b - A x
|
||||
Cheby(_Operator,r,AinvR); // A^{-1} ( b - A x ) ~ A^{-1} b - x
|
||||
out = out + AinvR;
|
||||
_Operator.HermOp(out,r);
|
||||
r = in - r; // b - A x
|
||||
RealD rr = norm2(r);
|
||||
RealD ss = norm2(in);
|
||||
std::cout << GridLogMessage << "ChebshevInverse resid " <<::sqrt(rr/ss)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// MIRS: CG on the infra-red shifted MdagM (HDCG, arXiv:1402.2585).
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
Field src(in.Grid());
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// Replay of a recorded GCR with a trivial preconditioner:
|
||||
// p_0 = r_0 ; x_{k+1} = x_k + a_k p_k ; r_{k+1} = r_k - a_k A p_k ;
|
||||
// p_{k+1} = r_{k+1} + sum_{j<northog(k)} b_kj p_{k-j}
|
||||
// One matvec per step, no reductions, history of mmax p vectors held
|
||||
// persistently (allocated on first use, per grid).
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field> class GCRReplaySmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
LinearOperatorBase<Field> &Linop;
|
||||
int mmax, nstep;
|
||||
std::vector<ComplexD> a;
|
||||
std::vector<std::vector<ComplexD> > b;
|
||||
GridBase *hist_grid = nullptr;
|
||||
std::vector<Field> p;
|
||||
GCRReplaySmoother(LinearOperatorBase<Field> &Op, const GCRCoefficients &c) : Linop(Op)
|
||||
{
|
||||
mmax = c.mmax; GRID_ASSERT(mmax>=1);
|
||||
nstep = c.Steps(); GRID_ASSERT(nstep>=1);
|
||||
a.resize(nstep); b.resize(nstep);
|
||||
for(int k=0;k<nstep;k++){
|
||||
a[k] = c.A(k);
|
||||
b[k].resize(c.NB(k));
|
||||
for(int j=0;j<c.NB(k);j++) b[k][j] = c.B(k,j);
|
||||
}
|
||||
std::cout << GridLogMessage << " GCRReplaySmoother: " << nstep << " steps, mmax " << mmax
|
||||
<< ", from " << c.Calls() << " recorded calls" << std::endl;
|
||||
}
|
||||
void operator() (const Field &src, Field &psi)
|
||||
{
|
||||
GRID_TRACE("GCRReplaySmoother");
|
||||
GridBase *grid = src.Grid();
|
||||
if ( hist_grid != grid ) {
|
||||
p.clear(); p.reserve(mmax);
|
||||
for(int i=0;i<mmax;i++) p.emplace_back(grid);
|
||||
hist_grid = grid;
|
||||
}
|
||||
Field r(grid), q(grid);
|
||||
r = src;
|
||||
psi = Zero();
|
||||
p[0] = r;
|
||||
for(int k=0;k<nstep;k++){
|
||||
int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
|
||||
Linop.Op(p[peri_k],q); // q_k = A p_k
|
||||
axpy(psi, a[k], p[peri_k], psi);
|
||||
if ( k==nstep-1 ) break;
|
||||
axpy(r, -a[k], q, r);
|
||||
p[peri_kp] = r;
|
||||
for(int j=0;j<(int)b[k].size();j++){
|
||||
int peri_back=(k-j)%mmax;
|
||||
axpy(p[peri_kp], b[k][j], p[peri_back], p[peri_kp]);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
// A LinearFunction that forwards to a replaceable target: lets a V-cycle be
|
||||
// built once around a smoother slot whose implementation is swapped at run
|
||||
// time (record with the adaptive GCR, then replay the polynomial).
|
||||
//////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field> class SwitchableSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
LinearFunction<Field> *current;
|
||||
std::string label;
|
||||
SwitchableSmoother(LinearFunction<Field> &initial, const std::string &l="initial") : current(&initial), label(l) {}
|
||||
void Set(LinearFunction<Field> &f, const std::string &l)
|
||||
{
|
||||
current = &f; label = l;
|
||||
std::cout << GridLogMessage << " SwitchableSmoother -> " << l << std::endl;
|
||||
}
|
||||
void operator() (const Field &in, Field &out) { (*current)(in,out); }
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
@@ -93,6 +93,18 @@ int FineSmootherOrder = 6;
|
||||
int FineSmootherMmax = 6;
|
||||
RealD CoarseSmootherShift = 0.1;
|
||||
int PowerIterations = 0; // >0: power-iterate the smoother operators before the solves (spectral edge)
|
||||
// Smoother implementation per level (Smoothers.h):
|
||||
// gcr : the adaptive PGCR (default, as always)
|
||||
// replay : run the PGCR with a coefficient recorder for the first
|
||||
// PolyRecordIters outer steps, then switch to GCRReplaySmoother
|
||||
// (same polynomial, no inner products) -- 1402.2585 p.13 revisited
|
||||
// cheb : ChebyshevNonHermitianSmoother, 1/x on [ChebLo,ChebHi], order
|
||||
// = the GCR step count of that level
|
||||
std::string FineSmootherMode = "gcr";
|
||||
std::string CoarseSmootherMode = "gcr";
|
||||
int PolyRecordIters = 4;
|
||||
RealD FineChebLo = 3.0, FineChebHi = 137.0; // from the harvested polynomial and PowerIterations edge
|
||||
RealD CoarseChebLo = 8.0, CoarseChebHi = 45.0;
|
||||
int CoarseSmootherNstep = 2;
|
||||
int CoarseSmootherMmax = 2;
|
||||
RealD CoarseSolverTol = 0.05;
|
||||
@@ -137,6 +149,13 @@ void ParseEnvironment(void)
|
||||
if(getenv("FineSmootherMmax")) FineSmootherMmax = atoi(getenv("FineSmootherMmax"));
|
||||
if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift"));
|
||||
if(getenv("PowerIterations")) PowerIterations = atoi(getenv("PowerIterations"));
|
||||
if(getenv("FineSmootherMode")) FineSmootherMode = getenv("FineSmootherMode");
|
||||
if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode");
|
||||
if(getenv("PolyRecordIters")) PolyRecordIters = atoi(getenv("PolyRecordIters"));
|
||||
if(getenv("FineChebLo")) FineChebLo = atof(getenv("FineChebLo"));
|
||||
if(getenv("FineChebHi")) FineChebHi = atof(getenv("FineChebHi"));
|
||||
if(getenv("CoarseChebLo")) CoarseChebLo = atof(getenv("CoarseChebLo"));
|
||||
if(getenv("CoarseChebHi")) CoarseChebHi = atof(getenv("CoarseChebHi"));
|
||||
if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep"));
|
||||
if(getenv("CoarseSmootherMmax")) CoarseSmootherMmax = atoi(getenv("CoarseSmootherMmax"));
|
||||
if(getenv("CoarseSolverTol")) CoarseSolverTol = atof(getenv("CoarseSolverTol"));
|
||||
@@ -346,6 +365,7 @@ public:
|
||||
std::string name = "Level 1";
|
||||
LinearOperatorBase<Field> &Linop;
|
||||
MrhsLinearFunction<Field> &Preconditioner;
|
||||
std::function<void(int)> OnStep; // called with the outer step count after every step (smoother switching)
|
||||
void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; }
|
||||
void Name(std::string n){ name = n; }
|
||||
void SetZeroGuess(int z){ ZeroGuess=z; }
|
||||
@@ -387,6 +407,7 @@ public:
|
||||
Preconditioner(r,p[0]); vOp(p[0],q[0]); qq[0]=vnorm2(q[0]); cp=vnorm2(r);
|
||||
for(int k=0;k<nstep;k++){
|
||||
steps++; int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
|
||||
if ( OnStep ) OnStep(steps);
|
||||
rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k];
|
||||
vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
|
||||
std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"] resid "<<cp<<" target "<<rsq<<std::endl;
|
||||
@@ -938,8 +959,9 @@ int main (int argc, char ** argv)
|
||||
CoarseSmootherGCR(0.01,1,ShiftedC,simpleC,CoarseSmootherMmax,CoarseSmootherNstep);
|
||||
CoarseSmootherGCR.Level(2); CoarseSmootherGCR.Name("Csmoother"); CoarseSmootherGCR.SetZeroGuess(1);
|
||||
|
||||
SwitchableSmoother<CoarseVector> CoarseSmootherSlot(CoarseSmootherGCR,"Csmoother GCR");
|
||||
MrhsCoarseThreeLevelPrec<CoarseVector,CoarseCoarseVector>
|
||||
L2to3Precon(LinOpC, CoarseSmootherGCR, MrhsProjectorL2, ccSolve,
|
||||
L2to3Precon(LinOpC, CoarseSmootherSlot, MrhsProjectorL2, ccSolve,
|
||||
Coarse5d, CoarseCoarse5d, CCMrhs, nr);
|
||||
|
||||
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>
|
||||
@@ -951,13 +973,51 @@ int main (int argc, char ** argv)
|
||||
SmootherGCR.LogCoefficients(SmootherCoeffLog);
|
||||
CoarseSmootherGCR.LogCoefficients(SmootherCoeffLog);
|
||||
|
||||
MrhsTwoLevelMG<LatticeFermionD,CoarseVector,FineSmoother_t>
|
||||
ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CMrhs);
|
||||
SwitchableSmoother<LatticeFermionD> FineSmootherSlot(SmootherGCR,"Fsmoother GCR");
|
||||
MrhsTwoLevelMG<LatticeFermionD,CoarseVector,SwitchableSmoother<LatticeFermionD> >
|
||||
ThreeLevelPrecon(PVdagM, FineSmootherSlot, MrhsProjector, L2PGCR, Coarse5d, CMrhs);
|
||||
|
||||
MrhsPGCRNonHermitian<LatticeFermionD>
|
||||
L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep);
|
||||
L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1);
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
// Smoother modes. Objects live for the duration of this RunSolve.
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
std::unique_ptr<ChebyshevNonHermitianSmoother<LatticeFermionD> > FineCheb;
|
||||
std::unique_ptr<ChebyshevNonHermitianSmoother<CoarseVector> > CoarseCheb;
|
||||
std::unique_ptr<GCRReplaySmoother<LatticeFermionD> > FineReplay;
|
||||
std::unique_ptr<GCRReplaySmoother<CoarseVector> > CoarseReplay;
|
||||
GCRCoefficients recF, recC;
|
||||
if ( FineSmootherMode == "cheb" ) {
|
||||
FineCheb.reset(new ChebyshevNonHermitianSmoother<LatticeFermionD>(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM));
|
||||
FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev");
|
||||
}
|
||||
if ( CoarseSmootherMode == "cheb" ) {
|
||||
CoarseCheb.reset(new ChebyshevNonHermitianSmoother<CoarseVector>(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC));
|
||||
CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev");
|
||||
}
|
||||
if ( FineSmootherMode == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF);
|
||||
if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC);
|
||||
L1PGCR.OnStep = [&](int step){
|
||||
if ( step != PolyRecordIters ) return;
|
||||
if ( FineSmootherMode == "replay" ) {
|
||||
SmootherGCR.SetCoefficientRecorder(nullptr);
|
||||
recF.Report("Fsmoother");
|
||||
FineReplay.reset(new GCRReplaySmoother<LatticeFermionD>(ShiftedPVdagM,recF));
|
||||
FineSmootherSlot.Set(*FineReplay,"Fsmoother replay");
|
||||
}
|
||||
if ( CoarseSmootherMode == "replay" ) {
|
||||
CoarseSmootherGCR.SetCoefficientRecorder(nullptr);
|
||||
recC.Report("Csmoother");
|
||||
CoarseReplay.reset(new GCRReplaySmoother<CoarseVector>(ShiftedC,recC));
|
||||
CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay");
|
||||
}
|
||||
};
|
||||
std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << " coarse " << CoarseSmootherMode
|
||||
<< (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? " (record for "+std::to_string(PolyRecordIters)+" outer steps)" : "")
|
||||
<< std::endl;
|
||||
|
||||
std::vector<LatticeFermionD> src(nr,FGrid), sol(nr,FGrid);
|
||||
for(int r=0;r<nr;r++){ gaussian(RNG5,src[r]); sol[r]=Zero(); }
|
||||
|
||||
|
||||
@@ -0,0 +1,96 @@
|
||||
#!/bin/bash -l
|
||||
#SBATCH --job-name=batched-linalg
|
||||
#SBATCH --nodes=1
|
||||
#SBATCH --ntasks-per-node=8
|
||||
#SBATCH --cpus-per-task=7
|
||||
#SBATCH --gpus-per-node=8
|
||||
#SBATCH --time=0:20:00
|
||||
#SBATCH --account=phy157_dwf
|
||||
#SBATCH --gpu-bind=none
|
||||
#SBATCH --exclusive
|
||||
#SBATCH --mem=0
|
||||
#SBATCH -S 0
|
||||
#SBATCH -q debug
|
||||
|
||||
##############################################################################
|
||||
# GPU gate for the batched Krylov linear algebra (Lattice_reduction.h:
|
||||
# innerProductMulti / axpyMulti / axpyMultiNorm with the by-value ViewPack
|
||||
# kernel argument) and the PGCR that sits on it. Nothing here has run on a
|
||||
# GPU before this job.
|
||||
#
|
||||
# B1 : Test_multi_reduction 8 ranks, 2x4 -- every batch width {2,4,8,16},
|
||||
# >16 chunking, fine (SpinColour) and coarse (iVector<CComplex,8> on a
|
||||
# 5D grid) fields, against the sequential reference, tol 1e-12.
|
||||
# B2 : Test_multi_reduction 1 rank -- same, no comms: separates a
|
||||
# kernel fault from a reduction/allreduce fault if B1 fails.
|
||||
# B3 : Test_pgcr_history 8 ranks -- solver on top: bitwise-
|
||||
# repeatable history reuse, true residual, second mmax.
|
||||
# B4 : Test_ring_allreduce 8 ranks, 2x4 -- P2P ring allreduce vs
|
||||
# MPI_Allreduce, all types/sizes, bitwise repeat, 16 MB timing.
|
||||
# (Also not yet run on Frontier; cheap to include.)
|
||||
#
|
||||
# Each stage is a separate srun with unbuffered output and its own timeout;
|
||||
# the summary at the end lists pass/fail per stage from the binaries' own
|
||||
# "ALL PASS"/"FAILURES" lines, and any GRID_ASSERT or signal.
|
||||
##############################################################################
|
||||
|
||||
cat << EOF > select_gpu
|
||||
#!/bin/bash
|
||||
export GPU_MAP=(0 1 2 3 7 6 5 4)
|
||||
export NUMA_MAP=(3 3 1 1 2 2 0 0)
|
||||
export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
|
||||
export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]}
|
||||
export HIP_VISIBLE_DEVICES=\$GPU
|
||||
unset ROCR_VISIBLE_DEVICES
|
||||
if [ \$SLURM_PROCID = "0" ]; then echo \$*; fi
|
||||
exec numactl -m \$NUMA -N \$NUMA \$*
|
||||
EOF
|
||||
chmod +x ./select_gpu
|
||||
|
||||
root=$HOME/ParallelIO/systems/Frontier
|
||||
source $root/sourceme-rocm7.2.sh
|
||||
|
||||
export OMP_NUM_THREADS=1
|
||||
export MPICH_GPU_SUPPORT_ENABLED=1
|
||||
export MPICH_SMP_SINGLE_COPY_MODE=CMA
|
||||
export MPICH_OFI_NIC_POLICY=GPU
|
||||
module load libfabric
|
||||
|
||||
OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000"
|
||||
T=$root/tests/debug
|
||||
|
||||
echo "========================================================="
|
||||
echo "binaries (git hash line is libGrid's build, not the test's)"
|
||||
echo "========================================================="
|
||||
ls -l $T/Test_multi_reduction $T/Test_pgcr_history $T/Test_ring_allreduce
|
||||
|
||||
run_stage () {
|
||||
name=$1; nranks=$2; mpi=$3; grid=$4; bin=$5
|
||||
echo "========================================================="
|
||||
echo "$name : $bin -n$nranks --mpi $mpi --grid $grid"
|
||||
echo "========================================================="
|
||||
timeout 300 srun -u -N1 -n$nranks ./select_gpu $T/$bin --mpi $mpi --grid $grid $OPTS1 2>&1 | tee $name.out \
|
||||
| grep -E "PASS|FAIL|T4 timing|GRID_ASSERT|Signal|signal|abort|Bus error|Segmentation"
|
||||
rc=${PIPESTATUS[0]}
|
||||
echo "$name exit code $rc (124 = timeout)"
|
||||
}
|
||||
|
||||
run_stage B1_multi_reduction_8 8 1.1.2.4 16.16.16.32 Test_multi_reduction
|
||||
run_stage B2_multi_reduction_1 1 1.1.1.1 16.16.16.16 Test_multi_reduction
|
||||
run_stage B3_pgcr_history_8 8 1.1.2.4 16.16.16.16 Test_pgcr_history
|
||||
run_stage B4_ring_allreduce_8 8 1.1.2.4 16.16.16.32 Test_ring_allreduce
|
||||
|
||||
echo "========================================================="
|
||||
echo "summary"
|
||||
echo "========================================================="
|
||||
for s in B1_multi_reduction_8 B2_multi_reduction_1 B3_pgcr_history_8 B4_ring_allreduce_8; do
|
||||
verdict=$(grep -oE "ALL PASS|FAILURES" $s.out | head -1)
|
||||
fails=$(grep -c "\*\* FAIL \*\*" $s.out)
|
||||
asserts=$(grep -cE "GRID_ASSERT|Signal|Bus error|Segmentation" $s.out)
|
||||
printf "%-24s %-10s fails=%s asserts/signals=%s\n" $s "${verdict:-NO-VERDICT}" $fails $asserts
|
||||
done
|
||||
echo "--- any FAIL lines ---"
|
||||
grep -h "\*\* FAIL \*\*" B*.out | head -20
|
||||
echo "--- ring allreduce timing (16 MB ComplexF) ---"
|
||||
grep -h "T4 timing" B4_ring_allreduce_8.out
|
||||
echo "========================================================="
|
||||
@@ -0,0 +1,149 @@
|
||||
#!/bin/bash -l
|
||||
#SBATCH --job-name=schur2d-ladder
|
||||
#SBATCH --nodes=36
|
||||
#SBATCH --ntasks-per-node=8
|
||||
#SBATCH --cpus-per-task=7
|
||||
#SBATCH --gpus-per-node=8
|
||||
#SBATCH --time=2:00:00
|
||||
#SBATCH --account=phy157_dwf
|
||||
#SBATCH --gpu-bind=none
|
||||
#SBATCH --exclusive
|
||||
#SBATCH --mem=0
|
||||
#SBATCH -S 0
|
||||
##SBATCH -q debug
|
||||
|
||||
##############################################################################
|
||||
# Staged shakeout of the 2D block-cyclic dense inverse (DENSE_SCHUR2D).
|
||||
#
|
||||
# F1 : unit tests, 1 node / 8 ranks (seconds)
|
||||
# F2 : scale rehearsal, 1 node, N=13824 (seconds)
|
||||
# F3 : scale rehearsal, 36 nodes, N=138240 (the timing number)
|
||||
# F4 : production example, 1D baseline then 2D, same nodes
|
||||
#
|
||||
# F4's 2D leg is gated on F3 passing: if the synthetic production-shape
|
||||
# inverse fails or is slow, do not spend the 400 s multigrid setup on it.
|
||||
##############################################################################
|
||||
|
||||
cat << EOF > select_gpu
|
||||
#!/bin/bash
|
||||
export GPU_MAP=(0 1 2 3 7 6 5 4)
|
||||
export NUMA_MAP=(3 3 1 1 2 2 0 0)
|
||||
export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
|
||||
export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]}
|
||||
export HIP_VISIBLE_DEVICES=\$GPU
|
||||
unset ROCR_VISIBLE_DEVICES
|
||||
|
||||
if [ \$SLURM_PROCID = "0" ]
|
||||
then
|
||||
echo \$*
|
||||
fi
|
||||
|
||||
exec numactl -m \$NUMA -N \$NUMA \$*
|
||||
EOF
|
||||
|
||||
chmod +x ./select_gpu
|
||||
|
||||
root=$HOME/ParallelIO/systems/Frontier
|
||||
source $root/sourceme-rocm7.2.sh
|
||||
|
||||
export OMP_NUM_THREADS=7
|
||||
export MPICH_GPU_SUPPORT_ENABLED=1
|
||||
export MPICH_SMP_SINGLE_COPY_MODE=CMA
|
||||
export MPICH_OFI_NIC_POLICY=GPU
|
||||
|
||||
module load libfabric
|
||||
|
||||
# Stale-environment protection: none of the gather knobs apply to the 2D
|
||||
# path, and DENSE_GATHER=1 without FORCE would (correctly) abort at start.
|
||||
unset DENSE_GATHER DENSE_GATHER_FORCE DENSE_GATHER_DEBUG DENSE_GATHER_MIN_BYTES
|
||||
|
||||
OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000"
|
||||
|
||||
##############################################################################
|
||||
echo "========================================================="
|
||||
echo "F1: unit tests, 1 node, 8 ranks"
|
||||
echo "========================================================="
|
||||
##############################################################################
|
||||
for t in Test_blockcyclic Test_summa Test_schur2d Test_schur2d_redist
|
||||
do
|
||||
echo "----- F1 $t -----"
|
||||
srun -N1 -n8 ./select_gpu $root/tests/debug/$t --mpi 1.1.2.4 --grid 16.16.16.16 $OPTS1
|
||||
done
|
||||
|
||||
##############################################################################
|
||||
echo "========================================================="
|
||||
echo "F2: scale rehearsal, 1 node, N=13824 (nb=1728, grid 2x4)"
|
||||
echo "========================================================="
|
||||
##############################################################################
|
||||
S2D_N=13824 srun -N1 -n8 ./select_gpu $root/tests/debug/Test_schur2d_scale \
|
||||
--mpi 1.1.2.4 --grid 16.16.16.16 $OPTS1
|
||||
|
||||
##############################################################################
|
||||
echo "========================================================="
|
||||
echo "F3: scale rehearsal, 36 nodes, N=138240 (nb=480, grid 16x18)"
|
||||
echo " THE number: invert phase vs the 1D baseline 328 s"
|
||||
echo "========================================================="
|
||||
##############################################################################
|
||||
S2D_N=138240 srun -N36 -n288 ./select_gpu $root/tests/debug/Test_schur2d_scale \
|
||||
--mpi 3.6.4.4 --grid 48.48.48.96 $OPTS1
|
||||
F3RC=$?
|
||||
echo "F3 exit code $F3RC"
|
||||
|
||||
##############################################################################
|
||||
echo "========================================================="
|
||||
echo "F4: production example (1D baseline first, then 2D if F3 passed)"
|
||||
echo "========================================================="
|
||||
##############################################################################
|
||||
|
||||
vol=48.48.48.96
|
||||
MPI_GEOM=3.6.4.4
|
||||
|
||||
export CoarseSmootherNstep=2
|
||||
export FineSmootherOrder=6
|
||||
|
||||
export SUBSPACE_FILE=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/subspace_nb64.scidac
|
||||
unset SLAB_FILE # force a fresh import in both legs
|
||||
|
||||
export DENSE_SCHUR=1
|
||||
export MASS=0.00078
|
||||
|
||||
export BLOCK=2.2.3.3 # L1->L2 blocking
|
||||
export BLOCK2=4.4.2.4 # L2->L3 blocking
|
||||
|
||||
export FineSmootherShift=0.1
|
||||
export CoarseSmootherShift=0.1
|
||||
export CoarseSolverTol=0.05
|
||||
export CoarseSolverOrder=200
|
||||
export L3_TOL=3.0e-1
|
||||
export L3_MAXIT=2
|
||||
export L3_NSTEP=50
|
||||
export OuterTol=1e-8
|
||||
export OuterMmax=4
|
||||
export OuterNstep=8
|
||||
export DENSE_CC=1
|
||||
export DENSE_APPLY_PROFILE=1
|
||||
unset DENSE_CC_CHECK
|
||||
export DENSE_SPLITK=128
|
||||
export DENSE_DEVICE_SUM=1
|
||||
export GRID_ALLOC_NCACHE_LARGE=64
|
||||
export NRHS=4
|
||||
|
||||
echo "----- F4a: 1D baseline (DENSE_SCHUR2D unset) -----"
|
||||
unset DENSE_SCHUR2D
|
||||
srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
|
||||
--mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap
|
||||
|
||||
if [ $F3RC -eq 0 ]
|
||||
then
|
||||
echo "----- F4b: 2D block-cyclic (DENSE_SCHUR2D=1) -----"
|
||||
export DENSE_SCHUR2D=1
|
||||
srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
|
||||
--mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap
|
||||
else
|
||||
echo "----- F4b SKIPPED: F3 failed (rc=$F3RC), not spending the setup on it -----"
|
||||
fi
|
||||
|
||||
echo "========================================================="
|
||||
echo "ladder complete; binaries were:"
|
||||
grep -h "git commit hash" slurm-$SLURM_JOB_ID.out 2>/dev/null | sort -u
|
||||
echo "========================================================="
|
||||
@@ -0,0 +1,125 @@
|
||||
#!/bin/bash -l
|
||||
#SBATCH --job-name=smoother-modes
|
||||
#SBATCH --nodes=36
|
||||
#SBATCH --ntasks-per-node=8
|
||||
#SBATCH --cpus-per-task=7
|
||||
#SBATCH --gpus-per-node=8
|
||||
#SBATCH --time=1:00:00
|
||||
#SBATCH --account=phy157_dwf
|
||||
#SBATCH --gpu-bind=none
|
||||
#SBATCH --exclusive
|
||||
#SBATCH --mem=0
|
||||
#SBATCH -S 0
|
||||
|
||||
##############################################################################
|
||||
# The 1402.2585 p.13 comparison on this machine: adaptive GCR smoothers vs
|
||||
# the same polynomials frozen (GCRReplaySmoother: recorded for the first
|
||||
# PolyRecordIters outer steps, then replayed with no inner products), vs a
|
||||
# Chebyshev 1/x fit on the measured interval. Same banked point otherwise.
|
||||
#
|
||||
# M1 gcr / gcr reference
|
||||
# M2 replay / replay both levels frozen after PolyRecordIters steps
|
||||
# M3 replay / gcr fine frozen only (the fine smoother is the reduction-
|
||||
# heavy one at Nrhs=1)
|
||||
# M4 cheb / gcr fine Chebyshev [FineChebLo,FineChebHi] order Fso
|
||||
# M5 cheb / cheb
|
||||
#
|
||||
# Laptop 8^4 findings (hot config, Ls=4, NBASIS=8): replay/replay converges
|
||||
# (28 vs 23 outer); cheb on the FINE level diverges there while cheb on the
|
||||
# coarse level is fine -- the tiny operator has modes the fixed polynomial
|
||||
# amplifies (left of / off the axis) that GCR handles adaptively. Production
|
||||
# spectrum is near-normal with edge 130.5 (shift 0.1) so M4/M5 may behave
|
||||
# differently; they are cheap to include and cheap to discard.
|
||||
#
|
||||
# Readouts per mode: Fouter count, s/RHS, and the FINAL exact-halo residual.
|
||||
##############################################################################
|
||||
|
||||
cat << EOF > select_gpu
|
||||
#!/bin/bash
|
||||
export GPU_MAP=(0 1 2 3 7 6 5 4)
|
||||
export NUMA_MAP=(3 3 1 1 2 2 0 0)
|
||||
export GPU=\${GPU_MAP[\$SLURM_LOCALID]}
|
||||
export NUMA=\${NUMA_MAP[\$SLURM_LOCALID]}
|
||||
export HIP_VISIBLE_DEVICES=\$GPU
|
||||
unset ROCR_VISIBLE_DEVICES
|
||||
if [ \$SLURM_PROCID = "0" ]; then echo \$*; fi
|
||||
exec numactl -m \$NUMA -N \$NUMA \$*
|
||||
EOF
|
||||
chmod +x ./select_gpu
|
||||
|
||||
root=$HOME/ParallelIO/systems/Frontier
|
||||
source $root/sourceme-rocm7.2.sh
|
||||
|
||||
export OMP_NUM_THREADS=7
|
||||
export MPICH_GPU_SUPPORT_ENABLED=1
|
||||
export MPICH_SMP_SINGLE_COPY_MODE=CMA
|
||||
export MPICH_OFI_NIC_POLICY=GPU
|
||||
module load libfabric
|
||||
unset DENSE_GATHER DENSE_GATHER_FORCE DENSE_GATHER_DEBUG DENSE_GATHER_MIN_BYTES
|
||||
|
||||
OPTS1="--accelerator-threads 8 --shm 4096 --shm-mpi 1 --device-mem 32000"
|
||||
vol=48.48.48.96
|
||||
MPI_GEOM=3.6.4.4
|
||||
|
||||
# banked point (2026-08-26 sweep): Css 2.0, Nstep 2, Fso 6, mmax 4, svm 8
|
||||
export FineSmootherShift=0.1
|
||||
export FineSmootherOrder=6
|
||||
export FineSmootherMmax=4
|
||||
export CoarseSmootherShift=2.0
|
||||
export CoarseSmootherNstep=2
|
||||
export CoarseSmootherMmax=2
|
||||
export CoarseSolverTol=0.05
|
||||
export CoarseSolverOrder=200
|
||||
export CoarseSolverMmax=8
|
||||
export OuterTol=1e-8
|
||||
export OuterMmax=6
|
||||
export OuterNstep=12
|
||||
export SUBSPACE_FILE=/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/subspace_nb64.scidac
|
||||
unset SLAB_FILE
|
||||
export DENSE_SCHUR=1
|
||||
export DENSE_SCHUR2D=1
|
||||
export MASS=0.00078
|
||||
export BLOCK=2.2.3.3
|
||||
export BLOCK2=4.4.2.4
|
||||
export L3_TOL=3.0e-1
|
||||
export L3_MAXIT=2
|
||||
export L3_NSTEP=50
|
||||
export DENSE_CC=1
|
||||
export DENSE_APPLY_PROFILE=1
|
||||
unset DENSE_CC_CHECK
|
||||
export DENSE_SPLITK=128
|
||||
export DENSE_DEVICE_SUM=2 # cartesian P2P ring: no 8 MB device-allreduce cliff (NRHS=12 abort)
|
||||
export GRID_ALLOC_NCACHE_LARGE=64
|
||||
export NRHS=4
|
||||
export PowerIterations=0
|
||||
export SmootherCoeffLog=0
|
||||
|
||||
# frozen-polynomial controls
|
||||
export PolyRecordIters=4 # outer steps recorded before the switch
|
||||
export FineChebLo=3.0 # harvested |R|<0.1 edge / PowerIteration edge x1.05
|
||||
export FineChebHi=137.0
|
||||
export CoarseChebLo=8.0
|
||||
export CoarseChebHi=47.0 # shift 2.0: edge 43.3 x1.08
|
||||
|
||||
run_mode () {
|
||||
name=$1; export FineSmootherMode=$2; export CoarseSmootherMode=$3
|
||||
echo "----- $name : FineSmootherMode=$FineSmootherMode CoarseSmootherMode=$CoarseSmootherMode -----"
|
||||
fname=log.modes.$name
|
||||
srun -N36 -n288 ./select_gpu $root/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix \
|
||||
--mpi ${MPI_GEOM} --grid $vol $OPTS1 --comms-overlap > $fname 2>&1
|
||||
echo " $(grep -h 'V2 3-level solve Nrhs' $fname | tr '\n' ' ')"
|
||||
echo " $(grep -h 'Fouter MrhsPGCR: Converged' $fname | sed 's/.*Converged/Converged/' | tr '\n' ' ')"
|
||||
echo " $(grep -h 'FINAL Nrhs .: worst' $fname | tr '\n' ' ')"
|
||||
grep -h "SwitchableSmoother\|GCRCoefficients .*calls" $fname | head -6
|
||||
}
|
||||
|
||||
run_mode M1_gcr_gcr gcr gcr
|
||||
run_mode M2_replay_replay replay replay
|
||||
run_mode M3_replay_gcr replay gcr
|
||||
run_mode M4_cheb_gcr cheb gcr
|
||||
run_mode M5_cheb_cheb cheb cheb
|
||||
|
||||
echo "========================================================="
|
||||
echo "summary"
|
||||
grep -h "V2 3-level solve Nrhs 1" log.modes.* | sed 's/.*V2/V2/'
|
||||
echo "========================================================="
|
||||
@@ -31,9 +31,6 @@ using namespace std;
|
||||
using namespace Grid;
|
||||
;
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
RealD SqrtApproximation(RealD x){
|
||||
return std::sqrt(x);
|
||||
}
|
||||
|
||||
@@ -112,9 +112,6 @@ void LoadBasis(aggregation &Agg, std::string file)
|
||||
#endif
|
||||
}
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
// Want Op in CoarsenOp to call MatPcDagMatPc
|
||||
template<class Field>
|
||||
@@ -132,19 +129,6 @@ public:
|
||||
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); }
|
||||
};
|
||||
/*
|
||||
template<class Field> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _SmootherOperator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
|
||||
@@ -291,63 +291,8 @@ public:
|
||||
};
|
||||
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
template<class Field> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _SmootherOperator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
// Field r(out.Grid());
|
||||
Cheby(_SmootherOperator,in,out);
|
||||
// _SmootherOperator.HermOp(out,r);
|
||||
// r=r-in;
|
||||
// RealD rr=norm2(r);
|
||||
// RealD ss=norm2(in);
|
||||
// std::cout << GridLogMessage<<" Chebyshev smoother resid "<<::sqrt(rr/ss)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template<class Field> class ChebyshevInverter : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _Operator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) :
|
||||
_Operator(Operator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field r(in.Grid());
|
||||
Field AinvR(in.Grid());
|
||||
_Operator.HermOp(out,r);
|
||||
r = in - r; // b - A x
|
||||
Cheby(_Operator,r,AinvR); // A^{-1} ( b - A x ) ~ A^{-1} b - x
|
||||
out = out + AinvR;
|
||||
_Operator.HermOp(out,r);
|
||||
r = in - r; // b - A x
|
||||
RealD rr = norm2(r);
|
||||
RealD ss = norm2(in);
|
||||
std::cout << "ChebshevInverse resid " <<::sqrt(rr/ss)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -291,62 +291,7 @@ public:
|
||||
};
|
||||
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
template<class Field> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _SmootherOperator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
// Field r(out.Grid());
|
||||
Cheby(_SmootherOperator,in,out);
|
||||
// _SmootherOperator.HermOp(out,r);
|
||||
// r=r-in;
|
||||
// RealD rr=norm2(r);
|
||||
// RealD ss=norm2(in);
|
||||
// std::cout << GridLogMessage<<" Chebyshev smoother resid "<<::sqrt(rr/ss)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field> class ChebyshevInverter : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _Operator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) :
|
||||
_Operator(Operator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field r(in.Grid());
|
||||
Field AinvR(in.Grid());
|
||||
_Operator.HermOp(out,r);
|
||||
r = in - r; // b - A x
|
||||
Cheby(_Operator,r,AinvR); // A^{-1} ( b - A x ) ~ A^{-1} b - x
|
||||
out = out + AinvR;
|
||||
_Operator.HermOp(out,r);
|
||||
r = in - r; // b - A x
|
||||
RealD rr = norm2(r);
|
||||
RealD ss = norm2(in);
|
||||
std::cout << "ChebshevInverse resid " <<::sqrt(rr/ss)<<std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -226,6 +226,7 @@ int main(int argc, char **argv)
|
||||
Stage("Grid Invert");
|
||||
{ GRID_TRACE("GridInvert"); RSI2.Invert(A); }
|
||||
double t3=usecond();
|
||||
RSI2.ReportTelemetry(grid); // SUMMA ring/gemm/alloc breakdown (outside the timed window)
|
||||
BlockCyclicRedistribute::CyclicToRows(grid,rowStart,A,&rows1d[0],myrows);
|
||||
double t4=usecond();
|
||||
double cert = Certify(grid,A0,A,nb,Pr,Pc);
|
||||
|
||||
@@ -87,9 +87,6 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
RealD SqrtApproximation(RealD x){
|
||||
return std::sqrt(x);
|
||||
}
|
||||
|
||||
@@ -48,64 +48,7 @@ using namespace Grid;
|
||||
* Lanczos:
|
||||
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
|
||||
*/
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
|
||||
@@ -48,9 +48,6 @@ using namespace Grid;
|
||||
* Lanczos:
|
||||
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
|
||||
*/
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field> class SolverWrapper : public LinearFunction<Field> {
|
||||
private:
|
||||
@@ -72,60 +69,6 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
|
||||
@@ -48,64 +48,7 @@ using namespace Grid;
|
||||
* Lanczos:
|
||||
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
|
||||
*/
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
|
||||
@@ -49,64 +49,7 @@ using namespace Grid;
|
||||
* Lanczos:
|
||||
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
|
||||
*/
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class RedBlackSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
|
||||
@@ -48,9 +48,6 @@ using namespace Grid;
|
||||
* Lanczos:
|
||||
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
|
||||
*/
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field> class SolverWrapper : public LinearFunction<Field> {
|
||||
private:
|
||||
@@ -72,60 +69,6 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
|
||||
@@ -48,64 +48,7 @@ using namespace Grid;
|
||||
* Lanczos:
|
||||
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
|
||||
*/
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
|
||||
@@ -112,64 +112,8 @@ public:
|
||||
};
|
||||
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
|
||||
|
||||
|
||||
@@ -109,63 +109,8 @@ public:
|
||||
};
|
||||
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & SmootherMatrix;
|
||||
FineOperator & SmootherOperator;
|
||||
RealD tol;
|
||||
RealD shift;
|
||||
int maxit;
|
||||
|
||||
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
|
||||
shift(_shift),tol(_tol),maxit(_maxit),
|
||||
SmootherOperator(_SmootherOperator),
|
||||
SmootherMatrix(_SmootherMatrix)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
ZeroGuesser<Field> Guess;
|
||||
ConjugateGradient<Field> CG(tol,maxit,false);
|
||||
|
||||
Field src(in.Grid());
|
||||
|
||||
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
|
||||
SmootherOperator.AdjOp(in,src);
|
||||
Guess(src,out);
|
||||
CG(MdagMOp,src,out);
|
||||
}
|
||||
};
|
||||
|
||||
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
|
||||
|
||||
|
||||
@@ -705,34 +705,7 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Fobj,class CComplex,int nbasis, class CoarseSolver>
|
||||
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
public:
|
||||
|
||||
@@ -725,34 +725,7 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
Matrix & _SmootherMatrix;
|
||||
FineOperator & _SmootherOperator;
|
||||
|
||||
Chebyshev<Field> Cheby;
|
||||
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
_SmootherMatrix(SmootherMatrix),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{};
|
||||
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
|
||||
_SmootherOperator.AdjOp(in,tmp);
|
||||
Cheby(MdagMOp,tmp,out);
|
||||
}
|
||||
};
|
||||
template<class Fobj,class CComplex,int nbasis, class CoarseSolver>
|
||||
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
public:
|
||||
|
||||
Reference in new issue
Block a user