(JEC) 25/3/15 first imprt

Makefile

+include $(SOPHYABASE)/include/sophyamake.inc
+
+# ===== Libsharp =====
+SHARPDIR = /Users/campagne/Travail/kits/libsharp-code/auto/
+SHARPLIB = $(SHARPDIR)/lib
+SHARPINC = $(SHARPDIR)/include
+SHARPLIBN = -lc_utils -lfftpack -lsharp
+
+CPPFLAGS  += -I$(SHARPINC) -I/opt/local/include
+LDFLAGS += -L$(SHARPLIB) $(SHARPLIBN) /opt/local/lib/libntl.a -lm
+
+OBJ = ./Objs/
+EXE = ./Objs/
+
+#  Define our target list
+.PHONY: default
+default: tstlag
+
+.PHONY: all 
+all : tstlag
+
+.PHONY: tidy
+tidy : 
+	rm *~
+
+.PHONY: clean
+clean :
+	rm -f ./Objs/*
+
+#C++ common Objects
+CXXOBJ = $(OBJ)laguerreBuilder.o \
+	$(OBJ)laguerreTransform.o \
+	$(OBJ)lagSphericTransform.o
+
+#C++ common Headers
+CXXHDR = laguerreBuilder.h \
+	laguerreTransform.h \
+	lagSphericTransform.h \
+	extrafunctions.h
+
+CXXFLAGS = -Wa,-q -march=native  -fopenmp -fno-common -O3
+CXXCOMPILE = $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c 
+
+#LDFLAGS +=  $(OPENMP) -lm
+
+
+#C++ rule for compiling
+$(OBJ)%.o: %.cc $(CXXHDR)
+	echo "compile... $@ depend $^"
+	$(CXXCOMPILE) -c $< -o $@
+
+######################
+.PHONY:  tstlag
+tstlag: $(EXE)tstlag
+	echo '---- tstlag made'
+
+$(EXE)tstlag : $(OBJ)tstlag.o $(CXXOBJ)
+	echo "Link..."
+	$(CXXLINK) -o $@ $(OBJ)tstlag.o  $(CXXOBJ) $(LDFLAGS) $(SOPHYAEXTSLBLIST)
+
+$(OBJ)tstlag.o: tstlag.cc $(CXXHDR)
+	echo "compile... $@"
+	$(CXXCOMPILE) $(FFLAGS) -c $< -o $@
+
+
+
+
+
+

extrafunctions.h

+#ifndef EXTRAFUNCTIONS_SEEN
+#define EXTRAFUNCTIONS_SEEN
+
+#ifdef SO_NTLXDBLE
+#include "NTL/xdouble.h" /* double : mantisse, long: exponent (NTL-library) */
+using namespace NTL;
+#endif
+
+
+
+namespace SOPHYA {
+//Utility fonctions to unify the code writing between "long double" (r_16) and NTL/xdouble
+#ifdef SO_NTLXDBLE  
+inline xdouble exp(const xdouble& x) {
+  return xexp(to_double(x));
+}
+#endif
+
+inline long double exp(long double x) {
+  return expl(x);
+}
+inline long double fabs(long double x) {
+  return fabsl(x);
+}
+
+
+
+inline double to_double(long double x) {
+  return double(x);
+}
+}//End namespace
+
+#endif

lagSphericTransform.cc

+#include <stdlib.h>
+#include <iostream>
+//Sophya
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include "pexceptions.h"
+
+#define LAGSPHERICTRANSFORM_CC_BFILE // avoid extern template declarations 
+#include "lagSphericTransform.h"
+
+#define DEBUG 0
+
+
+namespace SOPHYA {
+
+template <class T>
+void LaguerreSphericalTransform<T>::SetThetaPhiMap(string choice){
+  //Todo: on pourrait choisir ici son type de geoemetrie 2D: ici ECP = Fejer's 1 sum rule
+  
+  //libsharp defini Lmax = L-1 avec une somme l:0...Lmax inclus
+  //dans sharp_testsuite.c Fejer's 1st rule
+  //nrings = 2Lmax+1
+  //npixperring = 2Mmax+1
+  //et pour les Alm
+  //set_trinagular_alm (Lmax,Mmax).
+  
+  Nrings_ = 2*L_-1;
+  Nphi_   = 2*M_-1;
+  
+  sharpjob_.set_ECP_geometry(Nrings_,Nphi_);
+  sharpjob_.set_triangular_alm_info(L_-1,M_-1); 
+  
+  Npix_ = get_npix(sharpjob_.get_geom_info());    
+  Nalm_ = get_nalms(sharpjob_.get_alm_info());
+  
+  geomDone_ = true;
+  
+  cout << "SetThetaPhiMap: Nring, Nphi, Npix, Nalm " 
+       << Nrings_ << ", " 
+       << Nphi_ << ", "
+       << Npix_ << ", "
+       << Nalm_ 
+       << endl;
+}//SetThetaPhiMap
+
+template <class T>
+void LaguerreSphericalTransform<T>::Synthesis(const NDataBlock< complex<T> >& flmn, 
+					      NDataBlock<T>& fijk) {
+  int Ntot3DPix = Npix_*N_;
+  int NtotCoeff = Nalm_*N_;
+  if(flmn.Size() != NtotCoeff) throw SzMismatchError("SphLagTransform::Synthesis size error");
+    
+  fijk.ReSize(Ntot3DPix);
+
+  //Perform the Inverse Laguerre Transform first
+  NDataBlock< complex<T> > xflmk(NtotCoeff); //T version of the flmk
+    
+#if DEBUG >= 1
+  cout << "SphLagTransform::Synthesis type <T> I/O .... Go..." << endl;
+  cout << "Start   Inverse Laguerre...." << endl;
+#endif
+
+  lagTrans_.MultiSynthesis(flmn, xflmk, Nalm_);
+
+
+#if DEBUG >= 1
+  cout << "End   Inverse Laguerre...." << endl;
+  cout << "Start Inverse H.Spherical Transform " << endl;
+#endif
+
+
+  int off7n, off7k;
+  //perform the synthesis per shell
+  NDataBlock< complex<r_8> > flmk(Nalm_);    // r_8 version of the flmk for libsharp
+  NDataBlock<T> nodes(lagTrans_.GetNodes()); 
+
+  for(int n =0; n<N_; n++){
+    off7n = n*Nalm_;
+    off7k = n*Npix_; 
+
+    T scale = exp(-nodes(n)/2.); //<-------- use a scale to call sharp code
+
+#if DEBUG >= 1
+    cout << "Synthesis... at shell num: " << n << " scale = " << scale << endl;
+#endif
+      
+    //#pragma omp parallel for schedule(dynamic,1)
+    for(int i=0;i<Nalm_;i++){
+      xflmk(off7n+i) *= scale;
+      flmk(i) =  complex<r_8>(to_double(xflmk(off7n+i).real()),
+			      to_double(xflmk(off7n+i).imag()));
+#if DEBUG >= 2
+      cout << " Syntheis... flmk("<<i<<"): " << xflmk(off7n+i) << ", " << flmk(i) << endl;
+#endif
+    }
+
+    NDataBlock<r_8> rfijk(Npix_); //version for libsharp : todo see if declaration outside the loops is better
+    sharpjob_.alm2map(&(flmk.Data()[0].real()),&rfijk.Data()[0],false); //the false is to forbidd accumulation
+
+    //#pragma omp parallel for schedule(dynamic,1)
+    for(int i=0;i<Npix_;i++){
+      fijk(off7k+i) = rfijk(i); //store & switch r_8 -> T
+      fijk(off7k+i) /= scale;   //rescaling 
+#if DEBUG >= 2
+      cout << " Syntheis... rescaled fijk("<<i<<"): " << fijk(off7k+i) << " rfijk: "<< rfijk(i)  << endl;
+#endif      
+    }
+  }//loop on shell
+#if DEBUG >= 1
+  cout << "SphLagTransform::Synthesis.... Done" << endl;
+#endif
+}//Synthesis
+
+template <class T>
+void LaguerreSphericalTransform<T>::Analysis(const NDataBlock<T>& fijk, 
+					     NDataBlock< complex<T> >& flmn) {
+  int Ntot3DPix = Npix_*N_;
+  int NtotCoeff = Nalm_*N_;
+  if(fijk.Size() != Ntot3DPix) throw SzMismatchError("SphLagTransform::Analysis size error");
+  flmn.ReSize(NtotCoeff);
+
+#if DEBUG >= 1
+  cout << "SphLagTransform::Analysis xdbl.... Go..." << endl;
+#endif
+  //Perform the Spherical Hamonic analysis first
+    
+  NDataBlock< complex<T> > xflmk(NtotCoeff);
+  NDataBlock<T> nodes(lagTrans_.GetNodes());
+
+  NDataBlock< complex<r_8> > flmk(Nalm_); //used by lisharp
+  NDataBlock<r_8> rfijk(Npix_);           //used by libsharp
+
+  int off7n, off7k;
+  //Ylm analysis per shell
+  for(int n =0; n<N_; n++){
+    off7n = n*Nalm_;
+    off7k = n*Npix_; 
+
+#if DEBUG >= 1
+    cout << "Start Harmonic Spheric... at n = "<< n << endl;
+#endif
+    //rescalling des pixels values
+    T scale = exp(nodes(n)/2.);
+
+    //#pragma omp parallel for schedule(dynamic,1)
+    for(int i=0;i<Npix_;i++){
+      rfijk(i) = to_double(fijk(off7k+i)*scale); //store & switch T -> r_8 with a scale
+#if DEBUG >= 2
+      cout << "rfijk("<<i<<"): " << rfijk(i) << ", " << fijk(off7k+i) << endl;
+#endif
+    }
+     
+    sharpjob_.map2alm(&(rfijk.Data()[0]),&(flmk.Data()[0].real()),false); //the false is to forbidd accumulation
+    //#pragma omp parallel for schedule(dynamic,1)
+    for(int i=0; i<Nalm_; i++){
+      xflmk(i+off7n) = flmk(i);
+      xflmk(i+off7n) /= scale;  // rescaling
+#if DEBUG >= 2
+      cout << "flmk("<<i<<"): " << xflmk(i+off7n) << ", " << flmk(i+off7n) << endl;
+#endif
+    }
+  }//loop on shell
+
+#if DEBUG >= 1
+  cout << "End   Harmonic Spheric..." << endl;
+  cout << "Start Laguerre..." << endl;
+#endif
+
+
+  //Laguerre Transform
+  lagTrans_.MultiAnalysis(xflmk, flmn, Nalm_);
+
+#if DEBUG >= 1
+  cout << "End Laguerre..." << endl;
+#endif
+
+#if DEBUG >= 2
+  for (int i=0;i<NtotCoeff; i++) {
+    cout << "flmn("<<i<<"): " <<flmn(i) <<  endl;
+  }
+#endif
+#if DEBUG >= 1
+  cout << "SphLagTransform::Analysis.... Done" << endl;
+#endif
+}//Analysis
+
+///////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#ifdef SO_LDBLE128
+#pragma define_template LaguerreSphericalTransform<r_16>
+#endif
+#ifdef SO_NTLXDBLE
+#pragma define_template LaguerreSphericalTransform<xdouble>
+#endif
+#endif
+
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+#ifdef SO_LDBLE128
+template class  LaguerreSphericalTransform<r_16>;
+#endif
+#ifdef SO_NTLXDBLE
+template class LaguerreSphericalTransform<xdouble>;
+#endif
+#endif
+
+
+}// Fin de namespace

lagSphericTransform.h

+#ifndef LAGSPHERICTRANSFORM_SEEN
+#define LAGSPHERICTRANSFORM_SEEN
+
+//Sophya
+#include "machdefs.h"
+#include "ndatablock.h"
+
+//Sharp
+#include "sharp_lowlevel.h"
+#include "sharp_geomhelpers.h"
+#include "sharp_almhelpers.h"
+#include "sharp_cxx.h"
+
+// extern "C" {
+// #include <complex.h>
+// }
+
+
+#include "laguerreTransform.h"
+#include "extrafunctions.h"
+
+namespace SOPHYA {
+////////////////////////////////////////////////////////////////
+//// --------- Class Spherical Harmonic Laguerre Transform -- //
+////////////////////////////////////////////////////////////////
+
+
+template<class T>
+class LaguerreSphericalTransform {
+ public:
+  //! Creator
+  LaguerreSphericalTransform(int Lmax, int Mmax, int Nmax, r_8 R, int alpha=2): 
+    L_(Lmax), M_(Mmax), N_(Nmax), R_(R), alpha_(alpha), lagTrans_(Nmax,(T)R,alpha), 
+    Nrings_(0), Nphi_(0), Npix_(0), Nalm_(0),geomDone_(false) {
+    
+    if (Mmax<0) M_=Lmax;
+  }//Ctor
+    
+  //! Destructor
+  virtual ~LaguerreSphericalTransform() {}
+  
+  //! Number of colatitude rings of the 2D-sphere representation
+  int Nrings() { return Nrings_; }
+  //! Number of colatitude rings of the 2D-sphere representation
+  int NTheta() { return Nrings(); }
+  //! Number of phi values (pixels) per colatitude rings
+  int NPhi()   { return Nphi_; }
+  //! Total number of 2D pixels
+  int Npix()   { return Npix_; }
+  //! Total number of Alm coefficients
+  int Nalm()   { return Nalm_; }
+
+
+  //! Choose the 2D map representation & initialization of libsharp
+  void SetThetaPhiMap(string choice ="");
+
+
+ //! Synthesis 
+ /*! \brief Coeffs -> Pixels with Input/Output using T floating representation
+   \param flmn: 3D complex spherical-laguerre coefficients
+   \param fijk: 3D real function values
+ */
+  void Synthesis(const NDataBlock< complex<T> >& flmn, NDataBlock<T>& fijk);
+
+  //! Analysis
+  /*! \brief Pixels -> Coeffs with Input/Output using T floating representation
+    \param fijk: 3D real function values
+    \param flmn: 3D complex spherical-laguerre coefficients
+  */
+  void Analysis(const NDataBlock<T>& fijk, NDataBlock< complex<T> >& flmn);
+
+
+protected:
+
+  //from sharp test suite
+  static ptrdiff_t get_npix(const sharp_geom_info *ginfo) {
+    ptrdiff_t res=0;
+    for (int i=0; i<ginfo->npairs; ++i){
+      res += ginfo->pair[i].r1.nph;
+      if (ginfo->pair[i].r2.nph>0) res += ginfo->pair[i].r2.nph;
+    }
+    return res;
+  }
+  
+  //from sharp test suite
+  static ptrdiff_t get_nalms(const sharp_alm_info *ainfo) {
+    ptrdiff_t res=0;
+    for (int i=0; i<ainfo->nm; ++i)
+      res += ainfo->lmax-ainfo->mval[i]+1;
+    return res;
+  }
+
+protected:
+  int L_; //!< Spherical harmonic limit such that l:0,...,L-1   m:0,...,l
+  int M_; //!< maximum value of m
+  int N_; //!< Radial band limit such that n:0,...,N-1
+  r_8 R_; //!< Radial maximal value (not yet used)
+  int alpha_; //!< parameter of the generalized Laguerre Polynomials
+
+  LaguerreTransform<T> lagTrans_; //!< the Laguerre Transform
+
+  int Nrings_; //!< Number of colatitude rings of the 2D-sphere representation
+  int Nphi_; //!< Number of phi values (pixels) per colatitude rings
+  int Npix_; //!< Total number of 2D pixels
+  int Nalm_; //!< Total number of Alm coefficients
+  bool geomDone_; //!< Set to true once the geometry is properly set (See SetThetaPhiMap(string) )
+
+  sharp_cxxjob<r_8> sharpjob_; //!< libsharp job // voir si r_8 -> r_16 marcherait !!!
+  
+}; //class
+
+
+//--------- extern template declarations (if needed) -----------
+#if defined ( NEED_EXT_DECL_TEMP ) && !defined( LAGSPHERICTRANSFORM_CC_BFILE )
+
+#ifdef SO_LDBLE128
+extern template class LaguerreSphericalTransform<r_16>;
+#endif
+#ifdef SO_NTLXDBLE
+extern template class LaguerreSphericalTransform<xdouble>;
+#endif
+
+#endif // Fin de if defined ( NEED_EXT_DECL_TEMP )
+
+
+
+} //fin du namesape
+
+#endif // LAGSPHERICTRANSFORM

laguerreBuilder.cc

+#include <stdlib.h>
+#include <iostream>
+//Sophya
+#include "sopnamsp.h"
+#include "machdefs.h"
+#include "pexceptions.h"
+
+#define LAGUERREBUILDER_CC_BFILE // avoid extern template declarations 
+#include "laguerreBuilder.h"
+
+#define DEBUG 1
+
+namespace SOPHYA {
+
+template <class T>
+void LaguerrePolyQuad<T>::QuadWeightNodes(NDataBlock<T>& nodes, NDataBlock<T>& weights){
+  nodes.ReSize(N_);
+  weights.ReSize(N_);
+  T x;
+  T tN = (T)N_;
+
+  //tmp1 = Gamma(n+alpha+1)/(n! (n+1)^2) dans le cas ou alpha est un entier >0
+  //     = (n+alpha)!/(n! (n+1)^2) c'est le coeff du weight
+  T tmp1 = tN + (T)alpha_;
+  int alphaM1 = alpha_-1;
+  for (int i=1;i<alphaM1;i++) tmp1 *= (tmp1-(T)1.);
+  tmp1 /= (tN+1.);
+    
+  //Approximation by Sould & Secrest
+  for(int k=0;k<N_;k++){ 
+    if (k == 0) {
+      x = (T)((1.0+alpha_)*(3.+0.92*alpha_)/(1.+2.4*N_+1.8*alpha_));
+    } else if (k == 1) {
+      x += (T)((15.+6.25*alpha_)/(1.+0.9*alpha_+2.5*N_));
+    } else {//k>=2
+      int km1 = k-1;
+      T c0 = (T)((1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1));
+      T c1 = (T)(1.+0.3*alpha_);
+      x += c0*(x-nodes(k-2))/c1;
+    }
+
+    T xapp = x;
+      
+    //Refinement
+    RootUpdate(x);
+
+#if DEBUG >= 2
+    cout << " Update: " << x;
+#endif
+
+    nodes(k) = x;
+
+#if DEBUG >= 2
+    cout << " nodes: " << nodes(k);
+#endif	
+    N_++; //N+1
+    T tmp0 = Value(x); //L_{N+1}(x_i)
+    tmp0 *= tmp0; //[L_{N+1}(x_i)]^2
+    N_--; //N+1-1 = N
+    tmp0 = 1./tmp0; // 1./[L_{N+1}(x_i)]^2
+
+    weights(k) =  tmp1 * tmp0 * x;
+
+#if DEBUG >= 2
+    cout << " weight: " << weights(k) << endl;
+#endif      
+  }//loop k
+
+}//QuadWeightNodes
+
+
+template <class T>
+T LaguerrePolyQuad<T>::Value(T x, T scale) {
+  T Zero = (T)0.0;
+  T One  = (T)1.0;
+  
+  if(N_<0)  return Zero;
+  if(N_==0) return One*scale;
+  
+  T Two  = (T)2.0;
+  T alphaM1 = (T)alpha_ - One;
+  
+  T pnm1 = Zero;
+  T pn   = One*scale;
+  for(int i=1;i<=N_;i++){
+    T ti = (T)i;
+    T pnm2 = pnm1;
+    pnm1 = pn;
+    pn = ( (Two*ti + alphaM1 -x)*pnm1 - (ti + alphaM1)*pnm2 )/ti;
+  }
+  return pn;
+}//Value
+
+
+template <class T>
+void  LaguerrePolyQuad<T>::Values(T x, NDataBlock<T>& vals, T scale) {
+  if(N_<0) throw MathExc("Values call with N<=0");
+  if(vals.Size() != (N_+1) ) throw MathExc("Values with vals.size != (N+1)");
+  
+  T One  = (T)1.0;
+  
+  vals(0) =  One*scale;
+  if(N_==0) return;
+  
+  T Zero = (T)0.0;
+  T Two  = (T)2.0;
+  
+  T alphaM1 = (T)alpha_ - One;
+  
+  T pnm1 = Zero;
+  T pn   = One*scale;
+
+  for(int i=1;i<=N_;i++){
+    T ti = (T)i;
+    T pnm2 = pnm1;
+    pnm1 = pn;
+    pn = ( (Two*ti + alphaM1 -x)*pnm1 - (ti + alphaM1)*pnm2 )/ti;
+    vals(i) = pn;
+  }
+}//Values
+
+
+template <class T>
+T LaguerrePolyQuad<T>::LLpRatio(T x) {
+
+  if (N_<=0) throw MathExc("LLpRatio call with N<=0");
+  T pnm1 = (T)0.0;
+  T pn   = (T)1.0; //one can use here a scaling factor for numerical purpose
+  T One  = (T)1.0;
+  T Two  = (T)2.0;
+  T talpha = (T)alpha_;
+  T alphaM1 = talpha - One;
+  for(int i=1;i<=N_;i++){
+    T ti = (T)i;
+    T pnm2 = pnm1;
+    pnm1 = pn;
+    pn = ( (Two*ti + alphaM1 -x)*pnm1 - (ti + alphaM1)*pnm2 )/ti;
+  }
+  //pn = L_n^alpha(x), pnm1 = L_(n-1)^alpha(x)
+  T tN = (T)N_;
+  T ppn = (tN * pn - (tN+talpha)*pnm1)/x;
+
+  return pn/ppn;
+}//LLpRatio
+
+template <class T>
+void LaguerrePolyQuad<T>::RootUpdate(T& x) {
+  const int NLoops=100;
+  const T newtonEnd= (T)(1.e-10); 
+  T xold = (T)0;
+  T One = (T)1.0;
+  T talpha = (T)(alpha_);
+  T alphaP1 = One + talpha;
+  T tN = (T)N_;
+  T OneOv2 = (T)0.5;
+  //use simple Newton's method
+  int l;
+  for(l=0; l<NLoops; l++){
+    T LLpR = LLpRatio(x);
+    xold = x;
+    //      T eps = (One -OneOv2*LLpR*(One-alphaP1/xold))*LLpR;
+    T eps = LLpR;
+    x = xold - eps;
+    //      cout << "Update ("<<l<<"): " << xold << " " << x << endl;
+    if(fabs(x-xold)<newtonEnd) break;
+  } 
+}//RootUpdate
+
+///////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#ifdef SO_LDBLE128
+#pragma define_template LaguerrePolyQuad<r_16>
+#endif
+#ifdef SO_NTLXDBLE
+#pragma define_template LaguerrePolyQuad<xdouble>
+#endif
+#endif
+