portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Approx tests and lanczos improvement

Browse Source
This commit is contained in:
paboyle 2017-10-13 11:30:50 +01:00
parent 9fe6ac71ea
commit cb9ff20249
9 changed files with 823 additions and 712 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
lib/algorithms/approx/Chebyshev.h
View File

@ -83,8 +83,10 @@ namespace Grid {
public:
void csv(std::ostream &out){
RealD diff = hi-lo;
for (RealD x=lo-0.2*diff; x<hi+0.2*diff; x+=(hi-lo)/1000) {
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;
}

1399
lib/algorithms/iterative/BlockImplicitlyRestartedLanczos/BlockImplicitlyRestartedLanczos.h
View File

File diff suppressed because it is too large Load Diff

5
lib/algorithms/iterative/BlockImplicitlyRestartedLanczos/FieldBasisVector.h
View File

@ -48,7 +48,7 @@ class BasisFieldVector {
}
}
void rotate(std::vector<RealD>& Qt,int j0, int j1, int k0,int k1,int Nm) {
void rotate(Eigen::MatrixXd& Qt,int j0, int j1, int k0,int k1,int Nm) {
GridBase* grid = _v[0]._grid;
@ -62,7 +62,7 @@ class BasisFieldVector {
for(int j=j0; j<j1; ++j){
for(int k=k0; k<k1; ++k){
B[j] +=Qt[k+Nm*j] * _v[k]._odata[ss];
B[j] +=Qt(j,k) * _v[k]._odata[ss];
}
}
for(int j=j0; j<j1; ++j){
@ -70,7 +70,6 @@ class BasisFieldVector {
}
}
}
}
size_t size() const {

3
lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
View File

@ -244,6 +244,7 @@ until convergence
// Implicitly shifted QR transformations
Qt = Eigen::MatrixXd::Identity(Nm,Nm);
for(int ip=k2; ip<Nm; ++ip){
// Eigen replacement for qr_decomp ???
qr_decomp(eval,lme,Nm,Nm,Qt,eval2[ip],k1,Nm);
@ -319,7 +320,7 @@ until convergence
} // end of iter loop
std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
std::cout<< GridLogError <<" ImplicitlyRestartedLanczos::calc() NOT converged.";
std::cout << GridLogError <<" ImplicitlyRestartedLanczos::calc() NOT converged.";
std::cout << GridLogMessage <<"**************************************************************************"<< std::endl;
abort();

10
lib/log/Log.cc
View File

@ -59,13 +59,15 @@ void GridLogTimestamp(int on){
}
Colours GridLogColours(0);
GridLogger GridLogError(1, "Error", GridLogColours, "RED");
GridLogger GridLogIRL (1, "IRL" , GridLogColours, "NORMAL");
GridLogger GridLogSolver (1, "Solver", GridLogColours, "NORMAL");
GridLogger GridLogError (1, "Error" , GridLogColours, "RED");
GridLogger GridLogWarning(1, "Warning", GridLogColours, "YELLOW");
GridLogger GridLogMessage(1, "Message", GridLogColours, "NORMAL");
GridLogger GridLogDebug(1, "Debug", GridLogColours, "PURPLE");
GridLogger GridLogDebug (1, "Debug", GridLogColours, "PURPLE");
GridLogger GridLogPerformance(1, "Performance", GridLogColours, "GREEN");
GridLogger GridLogIterative(1, "Iterative", GridLogColours, "BLUE");
GridLogger GridLogIntegrator(1, "Integrator", GridLogColours, "BLUE");
GridLogger GridLogIterative (1, "Iterative", GridLogColours, "BLUE");
GridLogger GridLogIntegrator (1, "Integrator", GridLogColours, "BLUE");
void GridLogConfigure(std::vector<std::string> &logstreams) {
GridLogError.Active(0);

13
lib/log/Log.h
View File

@ -85,6 +85,7 @@ class Logger {
protected:
Colours &Painter;
int active;
int timing_mode;
static int timestamp;
std::string name, topName;
std::string COLOUR;
@ -101,20 +102,24 @@ public:
name(nm),
topName(topNm),
Painter(col_class),
timing_mode(0),
COLOUR(col) {} ;
void Active(int on) {active = on;};
int isActive(void) {return active;};
static void Timestamp(int on) {timestamp = on;};
void Reset(void) { StopWatch.Reset(); }
void TimingMode(int on) { timing_mode = on; if(on) Reset(); }
friend std::ostream& operator<< (std::ostream& stream, Logger& log){
if ( log.active ) {
stream << log.background()<< std::setw(8) << std::left << log.topName << log.background()<< " : ";
stream << log.colour() << std::setw(10) << std::left << log.name << log.background() << " : ";
stream << log.background()<< std::left << log.topName << log.background()<< " : ";
stream << log.colour() << std::left << log.name << log.background() << " : ";
if ( log.timestamp ) {
StopWatch.Stop();
GridTime now = StopWatch.Elapsed();
if ( log.timing_mode==1 ) StopWatch.Reset();
StopWatch.Start();
stream << log.evidence()<< now << log.background() << " : " ;
}
@ -135,6 +140,8 @@ public:
void GridLogConfigure(std::vector<std::string> &logstreams);
extern GridLogger GridLogIRL;
extern GridLogger GridLogSolver;
extern GridLogger GridLogError;
extern GridLogger GridLogWarning;
extern GridLogger GridLogMessage;

2
lib/threads/Threads.h
View File

@ -51,7 +51,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#define PARALLEL_CRITICAL
#endif
#define parallel_region PARALLEL_REGION
#define parallel_for PARALLEL_FOR_LOOP for
#define parallel_for_internal PARALLEL_FOR_LOOP_INTERN for
#define parallel_for_nest2 PARALLEL_NESTED_LOOP2 for
namespace Grid {

36
tests/debug/Test_cheby.cc
View File

@ -37,8 +37,15 @@ RealD InverseApproximation(RealD x){
RealD SqrtApproximation(RealD x){
return std::sqrt(x);
}
RealD Approximation32(RealD x){
return std::pow(x,-1.0/32.0);
}
RealD Approximation2(RealD x){
return std::pow(x,-1.0/2.0);
}
RealD StepFunction(RealD x){
if ( x<0.1 ) return 1.0;
if ( x<10.0 ) return 1.0;
else return 0.0;
}
@ -56,7 +63,6 @@ int main (int argc, char ** argv)
Chebyshev<LatticeFermion> ChebyInv(lo,hi,2000,InverseApproximation);
{
std::ofstream of("chebyinv");
ChebyInv.csv(of);
@ -78,7 +84,6 @@ int main (int argc, char ** argv)
ChebyStep.JacksonSmooth();
{
std::ofstream of("chebystepjack");
ChebyStep.csv(of);
@ -100,5 +105,30 @@ int main (int argc, char ** argv)
ChebyNE.csv(of);
}
lo=0.0;
hi=4.0;
Chebyshev<LatticeFermion> Cheby32(lo,hi,2000,Approximation32);
{
std::ofstream of("cheby32");
Cheby32.csv(of);
}
Cheby32.JacksonSmooth();
{
std::ofstream of("cheby32jack");
Cheby32.csv(of);
}
Chebyshev<LatticeFermion> ChebySqrt(lo,hi,2000,Approximation2);
{
std::ofstream of("chebysqrt");
ChebySqrt.csv(of);
}
ChebySqrt.JacksonSmooth();
{
std::ofstream of("chebysqrtjack");
ChebySqrt.csv(of);
}
Grid_finalize();
}

61
tests/hmc/Test_remez.cc
View File

@ -38,11 +38,11 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "Testing Remez"<<std::endl;
double lo=0.01;
double hi=1.0;
double lo=1.0e-3;
double hi=5.0;
int precision=64;
int degree=10;
AlgRemez remez(0.001,1.0,precision);
int degree=16;
AlgRemez remez(lo,hi,precision);
////////////////////////////////////////
// sqrt and inverse sqrt
@ -50,21 +50,50 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/2)"<<std::endl;
remez.generateApprox(degree,1,2);
MultiShiftFunction Sqrt(remez,1.0,false);
MultiShiftFunction InvSqrt(remez,1.0,true);
MultiShiftFunction Root2(remez,1.0,false);
MultiShiftFunction InvRoot2(remez,1.0,true);
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/4)"<<std::endl;
remez.generateApprox(degree,1,4);
MultiShiftFunction SqrtSqrt(remez,1.0,false);
MultiShiftFunction InvSqrtSqrt(remez,1.0,true);
MultiShiftFunction Root4(remez,1.0,false);
MultiShiftFunction InvRoot4(remez,1.0,true);
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/8)"<<std::endl;
remez.generateApprox(degree,1,8);
MultiShiftFunction Root8(remez,1.0,false);
MultiShiftFunction InvRoot8(remez,1.0,true);
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/16)"<<std::endl;
remez.generateApprox(degree,1,16);
MultiShiftFunction Root16(remez,1.0,false);
MultiShiftFunction InvRoot16(remez,1.0,true);
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/32)"<<std::endl;
remez.generateApprox(degree,1,32);
MultiShiftFunction Root32(remez,1.0,false);
MultiShiftFunction InvRoot32(remez,1.0,true);
ofstream gnuplot(std::string("Sqrt.gnu"),std::ios::out|std::ios::trunc);
Sqrt.gnuplot(gnuplot);
ofstream gnuplot(std::string("Root2.gnu"),std::ios::out|std::ios::trunc);
Root2.gnuplot(gnuplot);
ofstream gnuplot_i2(std::string("InvRoot2.gnu"),std::ios::out|std::ios::trunc);
InvRoot2.gnuplot(gnuplot_i2);
ofstream gnuplot_i4(std::string("InvRoot4.gnu"),std::ios::out|std::ios::trunc);
InvRoot4.gnuplot(gnuplot_i4);
ofstream gnuplot_i8(std::string("InvRoot8.gnu"),std::ios::out|std::ios::trunc);
InvRoot8.gnuplot(gnuplot_i8);
ofstream gnuplot_i16(std::string("InvRoot16.gnu"),std::ios::out|std::ios::trunc);
InvRoot16.gnuplot(gnuplot_i16);
ofstream gnuplot_i32(std::string("InvRoot32.gnu"),std::ios::out|std::ios::trunc);
InvRoot32.gnuplot(gnuplot_i32);
ofstream gnuplot_inv(std::string("InvSqrt.gnu"),std::ios::out|std::ios::trunc);
InvSqrt.gnuplot(gnuplot);
double x=0.6789;
double sx=std::sqrt(x);
@ -72,10 +101,10 @@ int main (int argc, char ** argv)
double isx=1.0/sx;
double issx=1.0/ssx;
double asx =Sqrt.approx(x);
double assx =SqrtSqrt.approx(x);
double aisx =InvSqrt.approx(x);
double aissx=InvSqrtSqrt.approx(x);
double asx =Root2.approx(x);
double assx =Root4.approx(x);
double aisx =InvRoot2.approx(x);
double aissx=InvRoot4.approx(x);
std::cout<<GridLogMessage << "x^(1/2) : "<<sx<<" "<<asx<<std::endl;
std::cout<<GridLogMessage << "x^(1/4) : "<<ssx<<" "<<assx<<std::endl;