(JEC) 09/07/2015 Laguerre->Bessel transformation first attempt

Makefile

@@ -18,6 +18,10 @@ else
 # ===== nNative BLAS (Darwin)
 	BLASYES = -DCBLAS -flax-vector-conversions 
 	BLASLIBN = -framework Accelerate
+	ifeq ($(BLAS), 0)
+		BLASYES =
+		BLASLIBN =
+	endif
 endif
 
 
@@ -92,13 +96,17 @@ CXXOBJ = $(OBJ)laguerreBuilder.o \
 	$(OBJ)laguerreTransform.o \
 	$(OBJ)lagsht_spheregeom.o \
 	$(OBJ)lagSphericTransform.o \
-	$(OBJ)walltimer.o
+	$(OBJ)walltimer.o \
+	$(OBJ)lagsht_bessel.o \
+	$(OBJ)laguerre2bessel.o
 
 CXXSHOBJ = laguerreBuilder.o \
 	laguerreTransform.o \
 	lagsht_spheregeom.o \
 	lagSphericTransform.o \
-	walltimer.o
+	walltimer.o \
+	lagsht_bessel.o \
+	laguerre2bessel.o	
 
 
 #C++ common Headers
@@ -111,16 +119,14 @@ CXXHDR = $(SRC)lagsht_exceptions.h \
 	 $(SRC)laguerreBuilder.h \
 	 $(SRC)laguerreTransform.h \
 	 $(SRC)lagSphericTransform.h \
-	 $(SRC)walltimer.h
-
-
-
+	 $(SRC)walltimer.h \
+	 $(SRC)lagsht_bessel.h \
+	 $(SRC)laguerre2bessel.h
 
 CPPFLAGS  += $(BLASYES) $(SHARPINC) $(BLASINC)
 LDFLAGS   += $(SHARPLIBN) $(BLASLIBN)  -lm
 
 
-
 #C++ rule for compiling
 $(OBJ)%.o: $(SRC)%.cc $(CXXHDR)
 	echo "compile... $<"

src/lagSphericTransform.cc

@@ -5,7 +5,6 @@
 
 #include "lagSphericTransform.h"
 #include "lagsht_exceptions.h"
-//#include "lagsht_healpixhelper.h"
 
 #include "walltimer.h"    //timing
 

src/lagSphericTransform.h

 #ifndef LAGSPHERICTRANSFORM_SEEN
 #define LAGSPHERICTRANSFORM_SEEN
-
 // #include <utility>      // std::pair, std::make_pair
 // using namespace std;
 
@@ -10,6 +9,7 @@
 // #include "sharp_almhelpers.h"
 // #include "sharp_cxx.h"
 
+
 #include "lagsht_spheregeom.h"
 #include "laguerreTransform.h"
 //#include "lagsht_geom.h"
@@ -22,14 +22,24 @@ namespace LagSHT {
 ////////////////////////////////////////////////////////////////
 //// --------- Class Spherical Harmonic Laguerre Transform -- //
 ////////////////////////////////////////////////////////////////
-
+/* Example to access the Laguerre Coeff flmn and the intermediate alm(r_k) Coeff.
+  for(int n=0;n<Nmax;n++){
+    for(int l=0;l<Lmax;l++){
+      for (int m=0;m<=l;m++) {
+	int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering 
+	r_8 flmn2 = flmn[id].real()*flmn[id].real() + flmn[id].imag()*flmn[id].imag();
+	r_8 almk2 = almk[id].real()*almk[id].real() + almk[id].imag()*almk[id].imag();
+      }
+    }
+  }
+ */
 class LaguerreSphericalTransform {
  public:
 
 
   LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax, 
 			     r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2): 
-    R_(R), Lmax_(Lmax), Mmax_(Mmax) {
+    Lmax_(Lmax), Mmax_(Mmax), R_(R) {
     
     //Factory
     if (spheregeom == "ECP" || spheregeom == "Fejer1") {
@@ -80,9 +90,13 @@ class LaguerreSphericalTransform {
   //! Laguerre Transform helper
   LaguerreTransform* GetLagTransform() const { return lagTrans_;}
 
+
+
    //! Getters
-  int GetOrder() { return N_; }
-  int GetAlpha() { return alpha_;}
+  int GetOrder()   const { return N_; }
+  int GetAlpha()   const { return alpha_; }
+  int GetNPixels() const { return Npix_; }
+  int GetNAlm()    const { return Nalm_; }
 
 
   //! Get the f(l,m,n) coefficient in the vector collection
@@ -121,7 +135,7 @@ class LaguerreSphericalTransform {
 
   //! Analysis
   /*! \brief Pixels -> Coeffs with Input/Output using T floating representation
-    \input fijk: 3D real function values
+     \input fijk: 3D real function values
     \output flmn: 3D complex spherical-laguerre coefficients
     \output almk: 2D alm complex spherical coefficients for each subshell k
   */
@@ -130,9 +144,7 @@ class LaguerreSphericalTransform {
 		vector< complex<r_8> >& almk
 		);
 
-
 protected:
-  r_8 R_; //!< Radial maximal value (not yet used)
 
   BaseGeometry* sphere_; //<! the 2D pixelization of the shells 
 
@@ -148,6 +160,9 @@ protected:
   int N_; //!< Order of the Generalized Lagurerre polynomial
   int alpha_;  //!< second parameter of the generalized Laguerre polynomial 
 
+  r_8 R_; //!< Radial maximal value (not yet used)
+
+
 }; //class
 
 

src/lagsht_bessel.cc

+#include "lagsht_bessel.h"
+#include "lagsht_exceptions.h"
+#include "laguerreTransform.h"
+
+namespace LagSHT {
+  
+void Bessel::BesselRoots() {
+  //  using namespace std;
+
+  string msg;
+  try {
+    unsigned int n_zeros = nmax_;
+    
+    for (int lcur =0; lcur<lmax_; lcur++){
+      //j(l,x) ~ J(l+1/2,x)/sqrt(x). 
+      r_8 order = (r_8)lcur + 0.5;
+      vector<r_8> zeros;
+      boost::math::cyl_bessel_j_zero(order,1,n_zeros, back_inserter(zeros));
+      copy(zeros.begin(),zeros.end(),qln_.begin() + lcur*nmax_);
+      
+    }//l loop   
+    
+  }  catch(LagSHTError& ex) {
+    msg = "BesselRoots ERROR: ";
+    msg += ex.what(); 
+    throw(msg);
+  }
+  catch (...) {
+    msg = "Bessel BesselRoots catched unknown (...) exception  ";
+    throw(msg);      
+  } 
+}//BesselRoots
+  
+}//namespace 

src/lagsht_bessel.h

+#ifndef LAGSHTBESSEL_SEEN
+#define LAGSHTBESSEL_SEEN
+//STD
+#include <vector>
+#include <algorithm>
+//Boost
+#include "boost/math/special_functions/bessel.hpp"
+
+#include "lagsht_numbers.h"
+using namespace std;
+
+namespace LagSHT {
+
+
+////////////////////////////////////////////////////////////////
+//// --------- Class Spherical Bessel utils                -- //
+////////////////////////////////////////////////////////////////
+
+
+class  Bessel {
+  
+ public:
+  //Ctor
+  Bessel(int lmax, int nmax): lmax_(lmax),nmax_(lmax) { 
+    Init(lmax,nmax); 
+  }
+  
+
+  //Dtor
+  virtual ~Bessel() {}
+
+  //Initializor
+  void Init(int lmax, int nmax) {
+    
+    lmax_= lmax;  
+    nmax_ = nmax;
+ 
+    //Set proper dimension of qln zeros matrix
+    qln_.resize(lmax*nmax);
+
+    //compute the Bessel roots
+    BesselRoots();
+
+  }//end Init
+
+  // Sperical Bessel functions
+  static inline r_8 jn(int n, r_8 x) { return boost::math::sph_bessel(n, x); }
+
+  //Helper
+  r_8 operator()(int l, int n) const { return qln_[l*nmax_+n];} //qln(l,n)
+  void GetVecRoots(std::vector<r_8>& roots, int l) const {
+    roots.resize(nmax_);
+    copy(qln_.begin()+l*nmax_,qln_.begin()+(l+1)*nmax_,roots.begin());
+  }
+
+  //Find Bessel root (qln) values such that j_l(qln) = 0 for l=O,...,Lmax and n=1,..,Nmax
+  void BesselRoots();
+
+ 
+ protected:
+
+  int lmax_;  //!< l=0,...,lmax-1  order of Spherical Bessel indice
+  int nmax_;  //!< n=0,...,nmax-1  index of the roots of Spherical Bessel func
+  
+  vector<r_8> qln_;  //!< roots: Care qln(l,n) l=0,...,lmax-1; n=0,...,nmax-1 but stands for n=1,..,nmax
+
+private:
+  
+  /*
+    Explicit prohibit Copy Ctor
+  */
+  Bessel(const Bessel& );
+  
+  /*
+    Explicit prohibit  Assignment Operator
+  */
+  Bessel& operator=(const Bessel&);
+}; //class
+
+}//end namespace
+#endif //LAGSHTBESSEL_SEEN

src/laguerre2bessel.cc

+#include "laguerre2bessel.h"
+#include "lagsht_bessel.h"
+#include "laguerreTransform.h"
+
+namespace LagSHT {
+
+Laguerre2Bessel::Laguerre2Bessel(int Lmax, int Mmax, int Nmax, int Pmax, r_8 Rmax):
+  Lmax_(Lmax), Mmax_(Mmax), Nmax_(Nmax), Pmax_(Pmax), Rmax_(Rmax) {
+
+    /*!compute J_{lpn} by Gauss-Laguerre quadrature. 
+      J_{lpn} = J_{ln}(k_{lp}) = \int_0^\infty dr r^2 {\cal L}^{(2)}_n (r) j_l(k_{lp}r)
+      J_{lpn} approximated by Sum_{k=0}^K0 w_k {\cal L}^{(2)}(r_k) j_l(k_{lp} r_k)
+      Use laguerreTransform class even if it computes more coeff as needed.
+    */
+
+    K0_ = 4*Nmax; //Heuristic Quadrature Order; should be > Nmax  
+    LaguerreTransform lagT(K0_,Rmax);
+    int stride = Lmax*Pmax; // nbre of (l,p) couples
+    vector<r_8> nodes = lagT.GetNodes(); //r_k of Gauss-Laguerre Quadrature
+    Bessel jnu(Lmax,Pmax); //compute  qlp: l:0...,Lmax-1 et p:0,...,Pmax-1
+    
+    vector<r_8> jlpk(stride*K0_); //values of j_l(k_{lp} r_k)
+    int lp=0;
+    for (int lc = 0; lc < Lmax; lc++) {
+      for (int pc = 0; pc < Pmax; pc++) {
+	r_8 klp = jnu(lc,pc)/Rmax;
+	for (int k=0; k<K0_ ; k++) {
+	  jlpk[lp+k*stride] = Bessel::jn(lc, klp * nodes[k]);
+	}//end k 
+	lp++;
+      }//end pc
+    }//end lc
+  
+    //J_{lpn} approximated by Sum_{k=0}^K0 w_k {\cal L}^{(2)}(r_k) j_l(k_{lp} r_k)
+    //    vector<r_8> Jlpn; //will be resize to stride*K0 elements while we need n:0...,Nmax-1 
+    lagT.MultiAnalysisR(jlpk,Jlpn_,stride);    
+}//Ctor
+
+/*!
+  FBlmp = FB_{lm}(k_{lp}) = Sqrt[2/Pi] Sum_n FL_{lmn} J_{lpn}
+*/
+void Laguerre2Bessel::Transform(const vector< complex<r_8> >& FLlmn,
+				int Nalm,
+				vector< complex<r_8> >& FBlmp){
+  
+  int Jstride = Lmax_*Pmax_;
+  r_8 fact = sqrt(2.0/M_PI);
+  //try here a pedestrian computation
+  for(int p=0; p<Pmax_; p++){
+    for(int l=0; l<Lmax_; l++){
+      int ploff7 = p + l*Pmax_;
+      for (int m=0; m<=l; m++) {
+	int almoff7 = l+m*Lmax_-m*(m+1)/2;
+	int idBess = p*Nalm + almoff7;
+	FBlmp[idBess] = 0.0;
+	for (int n=0; n<Nmax_; n++){
+	  int idLag =  n*Nalm + almoff7;
+	  FBlmp[idBess] += FLlmn[idLag]*Jlpn_[ploff7 + n*Jstride];
+	}//n
+	FBlmp[idBess] *= fact;
+      }//m
+    }//l
+  }//p
+  
+}//Transform
+
+
+void Laguerre2Bessel::Print(ostream& os, int what){
+  //debug of the Bessel-Laguerre scalar product
+  if(what>0){
+
+    int stride = Lmax_*Pmax_;
+    int lp=0;
+    for (int lc = 0; lc < Lmax_; lc++) {
+      for (int pc = 0; pc < Pmax_; pc++) {
+	for (int nc = 0; nc < Nmax_; nc++) {
+	  os << lc<<" "<<pc<<" "<<nc<<" " << setprecision(12)<< Jlpn_[lp + stride*nc] << endl;
+	}//loop nc
+	lp++;
+      }//loop pc
+    }//loop lc
+  }//what
+}//Print
+
+}//namespace

src/laguerre2bessel.h

+#ifndef LAGUERRE2BESSEL_SEEN
+#define LAGUERRE2BESSEL_SEEN
+
+#include <iostream>
+#include <vector>
+#include <complex>
+
+#include "lagsht_numbers.h"
+
+using namespace std;
+
+namespace LagSHT {
+
+////////////////////////////////////////////////////////////////////////////////
+//// --------- Class to transform Spherical Laguerre into Spherical Bessel -- //
+////////////////////////////////////////////////////////////////////////////////
+
+class Laguerre2Bessel {
+ public:
+
+  //!Ctor: TODO see if Mmax is finally used or not
+  Laguerre2Bessel(int Lmax, int Mmax, int Nmax, int Pmax, r_8 Rmax);
+
+  //! Transform Fourier-Laguerre into Fourier-Bessel
+  /*! \brief Fourier-Laguerre Coeffs -> Fourier-Bessel Coeffs
+    \input flmn: 3D complex spherical-laguerre coefficients (see LaguerreSphericalTransform)
+    \input Nalm: number of (l,m) couples to access the data
+    \output flmp: 3D complex spherical-bessel coefficients
+   */
+  void Transform(const vector< complex<r_8> >& FLlmn,
+		 int Nalm,
+		 vector< complex<r_8> >& FBlmp);
+
+  //!Debug
+  void Print(ostream& os, int what = -1);
+
+ 
+ protected:
+  int Lmax_; //!< Spherical harmonic limit such that l:0,...,Lmax-1   m:0,...,l
+  int Mmax_; //!< maximum value of m (NOT USED)
+  int Nmax_; //!< Max. Order of the Generalized Lagurerre polynomial  (alpha=2 by default)
+
+  int Pmax_; //!< Max. order of the Bessel coeff.
+  r_8 Rmax_; //!< Radial maximal value  
+
+  
+  int K0_; //!< Order of the LaguerreTransform to compute the Bessel-Laguerre scalar products Jlpn
+  vector<r_8> Jlpn_; //!< the  Bessel-Laguerre scalar product
+
+}; //class
+
+
+}//end namespace
+#endif

src/laguerreTransform.cc

@@ -91,7 +91,7 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
       vtmp.assign(N_,0.);
       for (int i = 0; i<N_; i++ ){
 	complex<r_8> fli = fi[l+ i*stride];
-	r_8 wi = weights_[i];
+	//////////JEC TO BE CHECK IF UNUSED	r_8 wi = weights_[i];
 	for (int n=0; n<N_; n++){
 	  //	  vtmp[n] +=  fli * wi * LnkMtx[n*N_+i];
 	  vtmp[n] +=  fli * LnkMtx[n+N_*i];
@@ -105,6 +105,65 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
 
 }//MultiAnalysis
 
+  //JEC 9/6/15 real case
+void LaguerreTransform::MultiAnalysisR(const vector<r_8>& fi, 
+				      vector<r_8>& fn, int stride) {
+
+  if(fi.size() != (size_t)(N_*stride) ) throw LagSHTError("LagTransform::Analysis size error");
+    fn.resize(N_*stride);
+
+#if DEBUG >= 1
+    cout << "Analysis with Function...." << endl;
+#endif
+#if DEBUG >= 2
+    std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
+#endif                     
+    r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!                   
+
+    vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
+    facts[0] = invalphaFact;
+    for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
+    
+    LaguerreFuncQuad lag(N_-1);
+    vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
+
+    vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
+    for (int k=0; k<N_; k++){
+      r_8 rk = nodes_[k];
+      lag.Values(rk,LnAll);
+      for (int n = 0; n<N_; n++ ){
+	//	LnkMtx[n*N_+k] = LnAll[n]*facts[n];
+	LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k];
+      }
+    }
+
+
+
+    //#ifdef CBLAS  
+//     cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fn[0]), 2*stride);
+//    throw LagSHTError("CBLAS to be adapted for REAL vectors");
+    //#else
+    vector<r_8> vtmp(N_);
+    for (int l=0; l<stride; l++) {
+      vtmp.assign(N_,0.);
+      for (int i = 0; i<N_; i++ ){
+	r_8 fli = fi[l+ i*stride];
+	//////////JEC TO BE CHECK IF UNUSED r_8 wi = weights_[i];
+	for (int n=0; n<N_; n++){
+	  //	  vtmp[n] +=  fli * wi * LnkMtx[n*N_+i];
+	  vtmp[n] +=  fli * LnkMtx[n+N_*i];
+	}//loop on k
+      }//loop on n
+      for (int n=0; n<N_; n++) {
+	fn[l+n*stride] += vtmp[n];
+      }
+    }//loop on l
+    //#endif
+
+}//MultiAnalysis
+
+
+
 
 void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn, 
 				       vector< complex<r_8> >& fi, int stride) {

src/laguerreTransform.h

@@ -61,6 +61,17 @@ class LaguerreTransform : public BaseLaguerreTransform {
   void MultiAnalysis(const vector< complex<r_8> >& fi, vector< complex<r_8> >& fn, int stride);
 
 
+
+  //! Multi Analysis at once (Values -> Coefficients)
+  /*!
+    \param fi: values f(r_i) with r_i zero L^(alpha)_N (Real case)
+    \param fn : result
+    \param stride : step between 2 conscecutive bunches of fi & fn 
+  */
+  void MultiAnalysisR(const vector<r_8>& fi, vector<r_8>& fn, int stride);
+
+
+
   //! Multiple Synthesis at once (Coefficientss -> Values)
   /*!
   \param fn: (input) complex coefficients of Laguerre Transform