(JEC) 23/11/15 avoid recomputation of the jlnp factors as much as possible

Makefile

@@ -52,6 +52,7 @@ all : confcheck makedir lib lagsht_testsuite
 install:
 	mkdir -p $(INCL)
 	cp $(SRC)/*.h $(INCL)
+ 
 
 
 .PHONY: check

setup.csh

+setenv LAGSHTROOTDIR /Users/campagne/Travail/Software/Laguerre/LagSHT

setup.sh

+LAGSHTROOTDIR=/Users/campagne/Travail/Software/Laguerre/LagSHT; export LAGSHTROOTDIR;

src/lagsht_testsuite.cc

@@ -25,11 +25,14 @@ struct PARAM {
   int Lmax;
   int N;
   r_8 R;
-  int Pmax; //JEC 7/10/15
+  int Pmax;
   int alpha;
   string geometry;
   int ntheta;
   int nphi;
+  string clenshawDir; //JEC 23/11/15
+  string jlnpDir;     //       "
+  bool recomputeJlnp; //       "
 } param ; 
 
 //simple random generator (taken from sharp_testsuite
@@ -293,6 +296,10 @@ void TestLaguerre2Bessel() {
   string geometry = param.geometry;
   int ntheta = param.ntheta;
   int nphi = param.nphi;
+  string clenshawDir = param.clenshawDir;   //JEC 23/11/15
+  string jlnpDir = param.jlnpDir;           //       "
+  bool recomputeJlnp = param.recomputeJlnp; //       "
+
 
   cout << " ______________ TestLaguerre2Bessel TEST ___________ " << endl;
   
@@ -374,7 +381,8 @@ void TestLaguerre2Bessel() {
   tstack_push("processing");
   cout << "tstack 2a Start" <<endl;
   tstack_push("Init Laguerre 2 Bessel");
-  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax);
+  Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax,
+			   clenshawDir,jlnpDir,recomputeJlnp);
   cout << "tstack 2a End" <<endl;
   tstack_pop("Init Laguerre 2 Bessel");
 
@@ -591,7 +599,7 @@ int main(int narg, char *arg[]) {
   int N = 128; 
   r_8 R = 1.;
   int Lmax = 128;
-  int Pmax = 0; //JEC 7/10/15
+  int Pmax = 0;
 
   int test=1;
 
@@ -599,6 +607,13 @@ int main(int narg, char *arg[]) {
   int ntheta = -1;
   int nphi = -1;
 
+  //JEC 23/11/15
+  char* pLAGSHTROOTDIR; pLAGSHTROOTDIR = getenv("LAGSHTROOTDIR");
+  string clenshawDir = string(pLAGSHTROOTDIR)+"/data";
+  string jlnpDir = string(pLAGSHTROOTDIR)+"/data";
+  bool recomputeJlnp = false;
+
+
   int ka=1;
   while (ka<narg) {
     if (strcmp(arg[ka],"-h")==0) {
@@ -644,6 +659,18 @@ int main(int narg, char *arg[]) {
       Pmax = atof(arg[ka+1]);
       ka+=2;      
     }    
+    else if (strcmp(arg[ka],"-clenPath")==0) {
+      clenshawDir = arg[ka+1];
+      ka+=2;
+    }
+    else if (strcmp(arg[ka],"-jlnpPath")==0) {
+      jlnpDir = arg[ka+1];
+      ka+=2;
+    }
+    else if(strcmp(arg[ka],"-jlnp")==0) {
+      recomputeJlnp = true;
+      ka+=1;
+    }
     else ka++;
   }//eo while
 
@@ -656,6 +683,11 @@ int main(int narg, char *arg[]) {
   param.ntheta = ntheta;
   param.nphi = nphi;
 
+  param.clenshawDir = clenshawDir;
+  param.jlnpDir     = jlnpDir;
+  param.recomputeJlnp = recomputeJlnp;
+
+
   cout << "Configuration parameters are set to: "
        << " Test number : " << test << "\n"
        << " Lmax, Nmax, Pmax, Rmax , alpha: " 
@@ -665,9 +697,14 @@ int main(int narg, char *arg[]) {
        << param.R << ", "
        << param.alpha << "\n"
        << "Geometry: " << param.geometry << " Ntheta, Nphi: " 
-       << param.ntheta << ", " << param.nphi
-       << endl;
-
+       << param.ntheta << ", " << param.nphi<<"\n";
+  if(recomputeJlnp){
+    cout << "rebuild Jlnp: \n Clenshaw Path<"<<param.clenshawDir<<">\n and jlnp Path<"
+	 << param.jlnpDir << ">"
+	 << endl;
+  } else {
+    cout << "Jlnp path <" << param.jlnpDir << ">" << endl;
+  }
 
   int rc=0;
   try {

src/laguerre2bessel.cc

@@ -7,6 +7,10 @@
 
 #include "walltimer.h"    //timing
 #include <fstream>       //dump
+#include <sstream>
+
+
+#define DEBUG 1
 
 namespace LagSHT {
 
@@ -16,25 +20,60 @@ namespace LagSHT {
      recursive
      Version 2: Jlnp calcule via une Discrete Fourier-Bessel Transform and a  Clenshaw-Curtis quadrature and global integration strategy. 
   */
-  Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
-				   LaguerreTransform* lagTrans,
-				   int Nmax, int Pmax, r_8 Rmax
-				   ):
-    sphere_(sphere),lagTrans_(lagTrans),
-    Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
-
-    if(sphere_) {
-      //store localy usefull parameters for the algorithms Analysis/Synthesis
-      Npix_ = sphere_->Npix();
-      Nalm_ = sphere_->Nalm();
-      Lmax_ = sphere_->Lmax();
-      Mmax_ = sphere_->Mmax();
-    } else {
-      throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
-    }
- 
-    ComputeJInteg();
-
+Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
+				 LaguerreTransform* lagTrans,
+				 int Nmax, int Pmax, r_8 Rmax,
+				 string clenshawDir,
+				 string jlnpDir,
+				 bool recomputeJlnp
+				 ):
+  sphere_(sphere),lagTrans_(lagTrans),
+  Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax), PPmax_(Pmax) {
+  
+  if(sphere_) {
+    //store localy usefull parameters for the algorithms Analysis/Synthesis
+    Npix_ = sphere_->Npix();
+    Nalm_ = sphere_->Nalm();
+    Lmax_ = sphere_->Lmax();
+    Mmax_ = sphere_->Mmax();
+  } else {
+    throw LagSHTError("LaguerreSphericalTransform::Ctor sphere_ ptr NULL");
+  }
+  
+  jnu_ = new Bessel(Lmax_,max(Pmax_,PPmax_));  //the max is here just in case PPmax =/= Pmax
+  
+  if(!recomputeJlnp){
+    
+    //read the Rmax-independant values of the jlnp and rescale
+    stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
+    
+    std::ifstream ifs(jlnpDir+"/"+ss.str(), std::ifstream::in);
+    if(!ifs.is_open())
+      throw LagSHTError("ERROR: jlnp file open failed.:"+ss.str());
+    
+    r_8 tau = lagTrans_->GetTau();
+    r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
+    int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
+    Jlnp_.resize(Nmax_*Jstride);
+    
+    int pcur, lcur, ncur;
+    r_8 val;
+    for(int p=0; p<Pmax_; p++){
+      for(int l=0; l<Lmax_; l++){
+	int ploff7 = p + l*Pmax_;
+	for (int n=0; n<Nmax_; n++){
+	  ifs >> lcur >> ncur >> pcur >> val;         //read
+// 	  if( (l!=lcur) || (n!=ncur) || (p!=pcur))                         //check
+// 	    throw LagSHTError("ERROR: jlnp file read failed.:"+ss.str());
+	  Jlnp_[ploff7 + n*Jstride] =  tau3sqrt * val;                     //rescaling
+	}//n
+      }//l
+    }//p
+    
+  } else {
+    ComputeJInteg(clenshawDir,jlnpDir);
+  }//jlnp omputation/load
+  
 }//Ctor
 
 
@@ -54,10 +93,12 @@ namespace LagSHT {
  JEC: 12/11/15 Add the contribution of Integral_{r_{Nmax-1}}^{r_{Nmax-1} + 3(r_{Nmax-1}-r_{Nmax-2})} by  Clenshaw-Curtis quadrature
 
  */
-void Laguerre2Bessel::ComputeJInteg(){
+void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
 
   r_8 tau = lagTrans_->GetTau();
   r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
+
+//   cout << "(JEC) tau3sqrt= " << tau3sqrt << endl;
   
   //JEC 12/11/15 Start
   r_8 rNm1 = Rmax_ / tau; // R = r[N-1] voir acces direct  = (lagTrans_->GetNodes()).back()
@@ -67,14 +108,14 @@ void Laguerre2Bessel::ComputeJInteg(){
   r_8 integUpperBound = rNm1+deltaR; //for the integral 
   r_8 tol = 1e-10;
   typedef  ClenshawCurtisQuadrature<double> Integrator_t;
-  Integrator_t theQuad(20,"/Users/campagne/Travail/Software/Laguerre/LagSHT/data/ClenshawCurtisRuleData-20.txt",false); //false=do not recompute the weights and nodes of the quadrature
+  string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-20.txt";
+  Integrator_t theQuad(20,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
   
   Quadrature<r_8,Integrator_t>::values_t integ_val;
   //JEC 12/11/15 end
   
   int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
-  Jlnp_.resize(Nmax_*Jstride);
-  
+  Jlnp_.resize(Nmax_*Jstride);  
 
   //JEC 16/11/15 START: Todo all that might be in common somewhere (hint: laguerreTransform.h)
   int alpha = lagTrans_->GetAlpha();
@@ -91,18 +132,16 @@ void Laguerre2Bessel::ComputeJInteg(){
   //END
   
     
-  int PPmax = 1024; // Sum maximum To be put in parameter
-  int PPmaxm1 = PPmax -1;
+  int PPmaxm1 = PPmax_ -1;
 
-  jnu_ = new Bessel(Lmax_,max(Pmax_,PPmax));  //To do put it in constructor
-  
-  
   tstack_push("J Integ compute ");
 
+  stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
   std::ofstream ofs;
-  ofs.open ("jlnp.txt", std::ofstream::out);
+  ofs.open (jlnpDir+"/"+ss.str(), std::ofstream::out);
 
-   
+
+//   cout << "(JEC) Jlpn recompute..." << endl;
   for(int p=0; p<Pmax_; p++){
     for(int l=0; l<Lmax_; l++){
       
@@ -123,7 +162,7 @@ void Laguerre2Bessel::ComputeJInteg(){
 	//contribution of the integral: [0,rNm1] estimated
 	//using Discrete Fourier-Bessel Transform
 
-	for (int pp=0; pp<PPmax; pp++){
+	for (int pp=0; pp<PPmax_; pp++){
 	  r_8 qlpp = (*jnu_)(l,pp);           //ql,pp
 	  r_8 jlp1 = Bessel::jn(l+1,qlpp);       //j(l+1,ql,pp)
 	  
@@ -135,23 +174,26 @@ void Laguerre2Bessel::ComputeJInteg(){
 	}//pp-loop
 	
 
-	Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
-	
-
+	//	Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
+	//JEC 23/11/15 save the Rmax-independant part
+	Jlnp_[ploff7 + n*Jstride] = coeff * facts[n] * sum;
 	
 	//contribution [rNm1, rNm1+deltaR] estimated with Clenshaw-Curtis quadrature (JEC 12/11/15) Start
 	IntFunc1D f(l,qlp/rNm1,Ln);
 	theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
 	integ_val = theQuad.GlobalAdapStrat(tol,2,1);
-	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (tau3sqrt*facts[n]*sqrt(2.0/M_PI));
+	//	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (tau3sqrt*facts[n]*sqrt(2.0/M_PI));
+	//JEC 23/11/15 save the Rmax-independant part
+	Jlnp_[ploff7 + n*Jstride] += integ_val.first * (facts[n]*sqrt(2.0/M_PI));
 	
-	//End
-
 		
+	//JEC 23/11/15 save the Rmax-independant part
 	ofs <<l<<" "
 	    <<n<<" "
-	    <<p<< " " << setprecision(12) << Jlnp_[ploff7 + n*Jstride] << endl;
-	
+	    <<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
+		
+	Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;
+
 	delete Ln;
 	
       }//n
@@ -161,7 +203,6 @@ void Laguerre2Bessel::ComputeJInteg(){
   ofs.close();
 
   tstack_pop("J Integ compute ");
-  cout << "J Integ compute  END" << endl;
   
 }//ComputeJlpnInteg
 
@@ -179,9 +220,10 @@ void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
   
   int Jstride = Lmax_*Pmax_;
 
+#if DEBUG>=1
   std::ofstream ofs;
   ofs.open ("FBlmp.txt", std::ofstream::out);
-
+#endif
 
   //try here a pedestrian computation
   for(int p=0; p<Pmax_; p++){
@@ -200,17 +242,19 @@ void Laguerre2Bessel::Lag2BessCoeff(const vector< complex<r_8> >& FLlmn,
 	}//n
 	FBlmp[idBess] = sum;
 
+#if DEBUG>=1
        	ofs <<l<<" "
 	    <<m<<" "
 	    <<p<< " " << setprecision(12) << 	FBlmp[idBess] << endl;
-
+#endif
 
       }//m
     }//l
   }//p
 
+#if DEBUG>=1
   ofs.close();
-
+#endif
   
 
 }//Transform
@@ -255,6 +299,7 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
   //pedestrian computation
   r_8 fact = sqrt(2.*M_PI)/(Rmax_*Rmax_*Rmax_);
 
+  vector< complex<r_8> > sumElmt(Pmax_);
   for(int k=0;k<Nmax_;k++){
 
     r_8 rk_red = vrk[k]/rNm1;
@@ -265,24 +310,6 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
 	int lmOff7 = l+m*Lmax_-m*(m+1)/2;
 	int id0 = k*Nalm_ + lmOff7;
 
-
-	//JEC 18/11/15 simple summation algo
-// 	complex<r_8> sum = 0.;
-// 	for(int p=0;p<Pmax_;p++){
-
-// 	  r_8 qlp = (*jnu_)(l,p);
-// 	  r_8 jlp1 = Bessel::jn(l+1,qlp);	  
-// 	  r_8 kapalp = fact/(jlp1*jlp1);
-
-// 	  int id1 = p*Nalm_ + lmOff7;
-
-//  	  complex<r_8> tmp = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
-// 	  sum += tmp;
-// 	}//loop on p
-// 	FBalmk[id0] = sum;
-
-	//JEC 18/11/15 Kahan summation algo
-	vector< complex<r_8> > sumElmt(Pmax_);
 	for(int p=0;p<Pmax_;p++){ //Fill
  	  r_8 qlp = (*jnu_)(l,p);
  	  r_8 jlp1 = Bessel::jn(l+1,qlp);	  
@@ -291,7 +318,7 @@ void Laguerre2Bessel::Synthesis(const vector< complex<r_8> >& FBlmp,
 
   	  sumElmt[p] = kapalp*Bessel::jn(l,qlp*rk_red)*FBlmp[id1];
 	}
-	//	FBalmk[id0] = LagSHT::accumulate(sumElmt.begin(),sumElmt.end());
+	//kahan summation	FBalmk[id0] = LagSHT::accumulate(sumElmt.begin(),sumElmt.end());
 	FBalmk[id0] = std::accumulate(sumElmt.begin(),sumElmt.end(),complex<r_8>(0.,0.));
 
       }//loop on m

src/laguerre2bessel.h

@@ -33,6 +33,7 @@
 #include "laguerreTransform.h"
 #include "lagsht_spheregeom.h"
 #include "lagsht_bessel.h"
+#include "lagsht_exceptions.h"
 
 #include "quadinteg.h" //JEC 12/11/15
 
@@ -44,27 +45,27 @@ namespace LagSHT {
 //// --------- Class to transform Spherical Laguerre into Spherical Bessel -- //
 ////////////////////////////////////////////////////////////////////////////////
 
-//implementation of the Kahan summation algorithm
-template <class InIt>
-typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
-    typedef typename std::iterator_traits<InIt>::value_type real;
-    real sum = real();
-    real running_error = real();
-    real temp;
-    real difference;
-
-    for (; begin != end; ++begin) {
-        difference = *begin;
-        difference -= running_error;
-        temp = sum;
-        temp += difference;
-        running_error = temp;
-        running_error -= sum;
-        running_error -= difference;
-        sum = std::move(temp);
-    }
-    return sum;
-}
+/* //implementation of the Kahan summation algorithm */
+/* template <class InIt> */
+/* typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) { */
+/*     typedef typename std::iterator_traits<InIt>::value_type real; */
+/*     real sum = real(); */
+/*     real running_error = real(); */
+/*     real temp; */
+/*     real difference; */
+
+/*     for (; begin != end; ++begin) { */
+/*         difference = *begin; */
+/*         difference -= running_error; */
+/*         temp = sum; */
+/*         temp += difference; */
+/*         running_error = temp; */
+/*         running_error -= sum; */
+/*         running_error -= difference; */
+/*         sum = std::move(temp); */
+/*     } */
+/*     return sum; */
+/* } */
 
 class Laguerre2Bessel {
 public:
@@ -72,7 +73,11 @@ public:
   //! Ctor
   Laguerre2Bessel(BaseGeometry* sphere, 
 		  LaguerreTransform* lagTrans,
-		  int Nmax, int Pmax, r_8 Rmax);
+		  int Nmax, int Pmax, r_8 Rmax,
+		  string clenshawDir,
+		  string jlnpDir,
+		  bool recomputeJlnp=false
+		  );
 
   //! Destructor
   virtual ~Laguerre2Bessel() {
@@ -108,8 +113,10 @@ private:
   /*! \brief
       J_{lnp} = J_{ln}(k_{lp}) 
               = tau^(3/2) \int_0^\infty dx x^2 {\cal L}^{(2)}_n (x) j_l(tau k_{lp} x)
+      \input  clenshawDir directory to find the Clenshaw-Curtis quadrature weight & nodes
+      \input jlnpDir Directory where to save the result
   */
-  void ComputeJInteg();
+  void ComputeJInteg(string clenshawDir, string jlnpDir);
 
  
   /*!
@@ -146,6 +153,8 @@ protected:
   int Pmax_; //!< Max. order of the Bessel coeff.
   r_8 Rmax_; //!< Radial maximal value  
 
+  int PPmax_; //!< cut off in the alm computation
+
   
   Bessel*<