Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 2c2d72af authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 25/4/15 transfert the test of nodes and weights in laguerreBuilder.h...

(JEC) 25/4/15 transfert the test of nodes and weights in laguerreBuilder.h nodes_and_weights; use the test on memory in the case:4 of lagsht_testsuite.cc to allow laguerreBuilder to be called with large values of n (tested with n:20000)
parent 4ffbdb5f
......@@ -411,16 +411,6 @@ int main(int narg, char *arg[]) {
int rc=0;
try {
if(N>5500) throw LagSHTError("WARNING: N>5500 the algorithm may fail due to weight estimates");
unsigned int estimateMem = 8*N*(uint_8)((2*Lmax+1)*(2*Lmax+1))/1e6; //Npix * 8Bytes and then in MBytes
if ( estimateMem > (0.9*maxmemsize) ) {
char buff[80];
cout << "estimate Mem usage " << estimateMem << " MB" << endl;
sprintf(buff,"Main FATAL: estimate Mem %u MB > (90%%) Max Mem %u MB ",estimateMem,maxmemsize);
throw LagSHTError(buff);
}//test memory
switch(test) {
......@@ -437,11 +427,25 @@ int main(int narg, char *arg[]) {
TestPixelization();
break;
case 4:
{
if(N>5500) throw LagSHTError("WARNING: N>5500 the algorithm may fail due to weight estimates");
unsigned int estimateMem = 8*N*(uint_8)((2*Lmax+1)*(2*Lmax+1))/1e6; //Npix * 8Bytes and then in MBytes
if ( estimateMem > (0.9*maxmemsize) ) {
char buff[80];
cout << "estimate Mem usage " << estimateMem << " MB" << endl;
sprintf(buff,"Main FATAL: estimate Mem %u MB > (90%%) Max Mem %u MB ",estimateMem,maxmemsize);
throw LagSHTError(buff);
}//test memory
tstack_push("MultiSphericalLaguerreTransform");
MultiSphericalLaguerreTransform();
tstack_pop("MultiSphericalLaguerreTransform");
tstack_report("MultiSphericalLaguerreTransform");
}
break;
default:
throw LagSHTError("EROOR Test type unkwown");
}//end of switch
......
......@@ -9,6 +9,9 @@
using namespace std;
#define DEBUG 1
namespace LagSHT {
......@@ -149,8 +152,7 @@ protected:
double x=0.;
// Approximation by Stroud & Secrest, Gaussian quadrature formulas, 1966
for(int k=0; k<N; k++)
{
for(int k=0; k<N; k++) {
if (k==0)
x = (1.0+alpha_)*(3.+0.92*alpha_)/(1.+2.4*N+1.8*alpha_);
else if (k==1)
......@@ -179,7 +181,45 @@ protected:
double tmp0 = value(x, N+1);
weights[k] = x*tmp1/(tmp0*tmp0);
}//end for
//verification
r_16 sumR = 0.;
r_16 sumW = 0.;
r_16 wip;
for (int i=0;i<N;i++){
sumR += nodes[i];
wip = exp(log(abs(weights[i]))-nodes[i]); if(weights[i]<0) wip *= -1.0;
sumW += wip;
// sumW += weights[i]*exp(-nodes[i]/tau);
}
int alphaFact = 1;
for(int i=1;i<=alpha_; i++) alphaFact *= i;
r_16 sumRTh = (r_16)(N*(N+alpha_));
r_16 sumWTh = (r_16)(alphaFact);
#if DEBUG >= 1
cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
#endif
if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< "WARNING!!!! verify LaguerrePolyQuad sum of nodes not accurate: "
<< " current value " << sumR << " theoretical value : " << sumRTh
<< " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< endl;
}
if (fabs(sumW-sumWTh)>(0.000001*sumWTh)) {
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< "WARNING!!!! verify LaguerrePolyQuad sum of weights not accurate: "
<< " current value " << sumW << " theoretical value : " << sumWTh
<< " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< endl;
}
}//nodes_and_weights
protected:
......@@ -199,5 +239,6 @@ private:
} // Fin du namespace
#undef DEBUG
#endif //LAGUERREBUILDER_SEEN
......@@ -32,41 +32,44 @@ LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha) :
LaguerreFuncQuad lag(N_,alpha_);
lag.QuadWeightNodes(nodes_,weights_);
//verification
r_16 sumR = 0.;
r_16 sumW = 0.;
r_16 tau = nodes_[N_-1]/708.0; //double number 708 ~ ln(minimum positive double)
for (int i=0;i<N_;i++){
sumR += nodes_[i];
sumW += weights_[i]*exp(-nodes_[i]/tau);
}
int alphaFact = 1;
for(int i=1;i<=alpha_; i++) alphaFact *= i;
r_16 sumRTh = (r_16)(N_*(N_+alpha_));
r_16 sumWTh = (r_16)(alphaFact)*pow(tau,(r_16)(alpha_+1));
#if DEBUG >= 1
cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
#endif
if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< "WARNING!!!! verify LaguerrePolyQuad sum of nodes not accurate: "
<< " current value " << sumR << " theoretical value : " << sumRTh
<< " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< endl;
}
if (fabs(sumW-sumWTh)>(0.000001*sumWTh)) {
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< "WARNING!!!! verify LaguerrePolyQuad sum of weights not accurate: "
<< " current value " << sumW << " theoretical value : " << sumWTh
<< " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< endl;
}
// //verification
// r_16 sumR = 0.;
// r_16 sumW = 0.;
// r_16 tau = nodes_[N_-1]/708.0; //double number 708 ~ ln(minimum positive double)
// for (int i=0;i<N_;i++){
// sumR += nodes_[i];
// sumW += weights_[i]*exp(-nodes_[i]/tau);
// }
// r_16 sumRTh = (r_16)(N_*(N_+alpha_));
// r_16 sumWTh = (r_16)(alphaFact)*pow(tau,(r_16)(alpha_+1));
// #if DEBUG >= 1
// cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
// cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
// #endif
// if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
// cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
// << "WARNING!!!! verify LaguerrePolyQuad sum of nodes not accurate: "
// << " current value " << sumR << " theoretical value : " << sumRTh
// << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
// << endl;
// }
// if (fabs(sumW-sumWTh)>(0.000001*sumWTh)) {
// cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
// << "WARNING!!!! verify LaguerrePolyQuad sum of weights not accurate: "
// << " current value " << sumW << " theoretical value : " << sumWTh
// << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
// << endl;
// }
alphaFact_ = sqrt((r_8)alphaFact); //sqrt(alpha!)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment