(JEC) 9/4/15 new program decoupled from SOphya

Makefile

+# ===== Libsharp =====
+SHARPDIR = /Users/campagne/Travail/kits/libsharp-code/auto/
+SHARPLIB = $(SHARPDIR)/lib
+SHARPINC = $(SHARPDIR)/include
+SHARPLIBN = -lsharp -lc_utils -lfftpack 
+
+
+OBJ = ./Objs/
+EXE = ./Objs/
+
+#  Define our target list
+.PHONY: default
+default: lagsht_testsuite
+
+.PHONY: all 
+all : lagsht_testsuite
+
+.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 = lagsht_exceptions.h \
+	lagsht_numbers.h \
+	lagsht_utils.h \
+	laguerreBuilder.h \
+	laguerreTransform.h \
+	lagSphericTransform.h
+
+
+CXX = g++
+CPPFLAGS = -DDarwin
+CXXFLAGS = -Wa,-q -march=native -fno-common -O3 -fopenmp
+
+CPPFLAGS  += -I$(SHARPINC) -I/opt/local/include
+LDFLAGS += -L$(SHARPLIB) $(SHARPLIBN) -lm
+
+
+CXXCOMPILE = $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c 
+CXXLINK = $(CXX) $(CXXFLAGS) -bind_at_load
+
+
+#C++ rule for compiling
+$(OBJ)%.o: %.cc $(CXXHDR)
+	echo "compile... $@ depend $^"
+	$(CXXCOMPILE) -c $< -o $@
+
+######################
+.PHONY:  lagsht_testsuite
+lagsht_testsuite: $(EXE)lagsht_testsuite
+	echo '---- lagsht_testsuite made'
+
+$(EXE)lagsht_testsuite : $(OBJ)lagsht_testsuite.o $(CXXOBJ)
+	echo "Link..."
+	$(CXXLINK) -o $@ $(OBJ)lagsht_testsuite.o  $(CXXOBJ) $(LDFLAGS)
+
+$(OBJ)lagsht_testsuite.o: lagsht_testsuite.cc $(CXXHDR)
+	echo "compile... $@"
+	$(CXXCOMPILE) -c $< -o $@
+
+
+
+
+
+

lagSphericTransform.cc

+#include <stdlib.h>
+#include <iostream>
+
+#include "lagSphericTransform.h"
+#include "lagsht_exceptions.h"
+
+#define DEBUG 0
+
+namespace LagSHT {
+
+void BaseLagSphTransform::SetThetaPhiMap(string choice){
+  //Todo: here to be developped the way to initialize libsharp
+  //the arguments can give the user the choice of the geometry
+  // this may be inspired from the "get_infos" routine in sharp_testsuite.c
+  
+  //libsharp defini Lmax = L-1 avec une somme l:0...Lmax inclus
+  //dans sharp_testsuite.c Fejer's 1st rule
+  //nrings = 2Lmax+1
+  //npixper ring = 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
+
+
+
+void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn, 
+					   vector<r_8>& fijk) {
+  int Ntot3DPix = Npix_*N_;
+  int NtotCoeff = Nalm_*N_;
+  if(flmn.size() != NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
+    
+  fijk.resize(Ntot3DPix);
+
+  //Perform the Inverse Laguerre Transform first
+  vector< complex<r_8> > flmk(NtotCoeff); 
+
+#if DEBUG >= 1
+  cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
+  cout << "Start   Inverse Laguerre...." << endl;
+#endif
+
+  lagTrans_.MultiSynthesis(flmn, flmk, 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
+  for(int n =0; n<N_; n++){
+    off7n = n*Nalm_;
+    off7k = n*Npix_; 
+
+#if DEBUG >= 2
+    cout << "Synthesis... at shell num: " << n << endl;
+#endif
+      
+    sharpjob_.alm2map(&(flmk.data()[off7n].real()),&fijk.data()[off7k],false);
+
+  }//loop on shell
+#if DEBUG >= 1
+  cout << "SphLagTransform::Synthesis.... Done" << endl;
+#endif
+}//Synthesis
+
+
+
+void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk, 
+					  vector< complex<r_8> >& flmn) {
+  int Ntot3DPix = Npix_*N_;
+  int NtotCoeff = Nalm_*N_;
+  if(fijk.size() != Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error");
+  flmn.resize(NtotCoeff);
+
+#if DEBUG >= 1
+  cout << "SphLagTransform::Analysis.... Go..." << endl;
+#endif
+  //Perform the Spherical Hamonic analysis first
+    
+  vector< complex<r_8> > flmk(NtotCoeff);
+
+  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
+
+    sharpjob_.map2alm(&(fijk.data()[off7k]),&(flmk.data()[off7n].real()),false); //the false is to forbidd accumulation
+  }//loop on shell
+
+#if DEBUG >= 1
+  cout << "End   Harmonic Spheric..." << endl;
+  cout << "Start Laguerre..." << endl;
+#endif
+
+
+  //Laguerre Transform
+  lagTrans_.MultiAnalysis(flmk, 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
+
+
+
+}// Fin de namespace

lagSphericTransform.h

+#ifndef LAGSPHERICTRANSFORM_SEEN
+#define LAGSPHERICTRANSFORM_SEEN
+//Sharp
+#include "sharp_lowlevel.h"
+#include "sharp_geomhelpers.h"
+#include "sharp_almhelpers.h"
+#include "sharp_cxx.h"
+
+
+#include "laguerreTransform.h"
+#include "lagsht_exceptions.h"
+
+
+namespace LagSHT {
+
+
+////////////////////////////////////////////////////////////////
+//// --------- Class Spherical Harmonic Laguerre Transform -- //
+////////////////////////////////////////////////////////////////
+
+class BaseLagSphTransform {
+ public:
+
+  typedef pair<r_8,r_8> pixpos_t;
+  typedef vector<pixpos_t> geom_t;
+
+  //! Creator
+  BaseLagSphTransform(int Lmax, int Mmax, int Nmax, int alpha=2):
+    L_(Lmax), M_(Mmax), N_(Nmax), alpha_(alpha),
+    Nrings_(0), Nphi_(0), Npix_(0), Nalm_(0),geomDone_(false) {
+     if (Mmax<0) M_=Lmax;
+  }
+  //! Destructor
+  virtual ~BaseLagSphTransform() {}
+
+  //! 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_; }
+  
+  //! List of pixel positions
+  geom_t GetPixelPositions() {
+    throw LagSHTError("BaseLagSphTransform::GetPixelPositions not yet implemented");
+  }
+  
+  
+  //! Choose the 2D map representation & initialization of libsharp
+  void SetThetaPhiMap(string choice ="");
+
+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
+  int alpha_; //!< parameter of the generalized Laguerre Polynomials
+
+  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 extend the validity of the computation
+
+};
+
+class LaguerreSphericalTransform : public BaseLagSphTransform {
+ public:
+
+  //! Constructor
+  LaguerreSphericalTransform(int Lmax, int Mmax, int Nmax, r_8 R, int alpha=2): 
+    BaseLagSphTransform(Lmax,Mmax,Nmax,alpha), R_(R), lagTrans_(Nmax,(r_8)R,alpha) {
+  }//Ctor
+
+    
+  //! Destructor
+  virtual ~LaguerreSphericalTransform() {}
+  
+
+ //! 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 vector< complex<r_8> >& flmn, vector<r_8>& 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 vector<r_8>& fijk, vector< complex<r_8> >& flmn);
+
+
+protected:
+  r_8 R_; //!< Radial maximal value (not yet used)
+  LaguerreTransform lagTrans_; //!< the Laguerre Transform
+  
+}; //class
+
+
+}//end namespace
+#endif
+ 

lagsht_exceptions.h

+#ifndef LAGSHTEXCEPTIONS_SEEN
+#define LAGSHTEXCEPTIONS_SEEN
+
+#include <exception>
+#include <string>
+
+namespace LagSHT {
+
+class LagSHTError : public std::exception {
+public:    
+  explicit LagSHTError(const std::string& m) :  msg_(m) { }
+  explicit LagSHTError(const char* m): msg_(m) { }
+  virtual ~LagSHTError() throw() {}
+  //! Implementation of std::exception what() method, returning the exception message
+  virtual const char* what() const throw() { return msg_.c_str(); }
+  
+private:
+  std::string msg_;
+  
+};//LagSHTError
+
+}//end namespace
+
+#endif //LAGSHTEXCEPTIONS_SEEN

lagsht_numbers.h

+#ifndef LAGSHTNUMBERS_SEEN
+#define LAGSHTNUMBERS_SEEN
+
+typedef signed char int_1;
+typedef unsigned char uint_1;
+typedef short int_2;
+typedef unsigned short uint_2;
+typedef int int_4;
+typedef unsigned int uint_4;
+typedef long long int_8;
+typedef unsigned long long uint_8;
+typedef float r_4;
+typedef double r_8;
+typedef long double r_16;
+
+
+#endif //LAGSHTNUMBERS_SEEN

lagsht_testsuite.cc

+#include <stdlib.h>
+#include <iostream>
+#include <fstream>
+#include <string.h>
+#include <limits>       // std::numeric_limits
+#include <math.h> //log
+#include <string>
+#include <typeinfo>
+using namespace std;
+
+#include "lagsht_exceptions.h"
+#include "lagsht_utils.h" //getMemorySize
+
+
+#include "lagSphericTransform.h"
+using namespace LagSHT;
+
+#define DEBUG 1
+
+//-------- Parameters set in the main and  used in the different test functions
+struct PARAM {
+  int Lmax;
+  int N;
+  r_8 R;
+  int alpha;
+} param ; 
+
+//simple random generator (taken from sharp_testsuite
+static double drand (double min, double max, int *state){
+  *state = (((*state) * 1103515245) + 12345) & 0x7fffffff;
+  return min + (max-min)*(*state)/(0x7fffffff+1.0);
+}//rand
+
+
+//----------------------------------------------------
+void TestBasic() {
+
+  int N = param.N;
+  int alpha = param.alpha;
+  int Lmax = param.Lmax;
+
+
+
+  cout << " ___________  BASIC TEST _____________ " << endl;
+  
+
+  cout << " sizeof(double): " << sizeof(r_8) << " bytes" 
+       << " min " << std::numeric_limits<r_8>::min() << "-ln(min): " << -log(std::numeric_limits<r_8>::min()) 
+       << " max " << std::numeric_limits<r_8>::max()  
+       << " min exp_10 " << std::numeric_limits<r_8>::min_exponent10
+       << " max exp_10 " << std::numeric_limits<r_8>::max_exponent10  
+       <<  endl;
+  
+
+
+  cout << " sizeof(long double): " << sizeof(r_16) << " bytes" 
+       << " min " << std::numeric_limits<r_16>::min() <<  "-ln(min): " << -logl(std::numeric_limits<r_16>::min()) 
+       << " max " << std::numeric_limits<r_16>::max()  
+       << " min exp_10 " << std::numeric_limits<r_16>::min_exponent10
+       << " max exp_10 " << std::numeric_limits<r_16>::max_exponent10  
+       <<  endl;
+}//TestBasic
+
+//----------------
+void LaguerreQuadrature() {
+
+  int N = param.N;
+  int alpha = param.alpha;
+  int Lmax = param.Lmax;
+
+  cout << " ___________  LaguerreQuadrature  TEST _____________ " << endl;
+
+  LaguerreFuncQuad lagFunc(N,alpha);
+  cout << "N = " << lagFunc.GetOrder() << " alpha = " << lagFunc.GetAlpha() << endl;
+  vector<r_16> nodes;
+  vector<r_16> weights;
+  lagFunc.QuadWeightNodes(nodes,weights);
+
+  char buff1[80], buff2[80];
+  ofstream ofs1, ofs2;
+  sprintf(buff1,"lagNodes-%d-r16-Func.txt",N);
+  sprintf(buff2,"lagWeights-%d-r16-Func.txt",N);
+  ofs1.open(buff1,ofstream::out);
+  ofs2.open(buff2,ofstream::out);
+  for (int i=0;i<N;i++){
+    ofs1 << nodes[i] << endl;
+    ofs2 << log10(weights[i]) << endl;
+  }
+
+}//LaguerreQuadrature
+
+//--------------------------------------
+void MultiLaguerreTransformSdtAlone() {
+  int N = param.N;
+  int alpha = param.alpha;
+  int Lmax = param.Lmax;
+  r_8 R  = param.R;
+
+  cout << " ___________   MultiLaguerreTransformSdtAlone  TEST _____________ " << endl;
+
+  //    Timer tm("Test 3 (r_16)");
+      
+  int mapSize=2; //this is the stride
+      
+  int Ntot = N*mapSize;
+
+#if DEBUG >= 1
+  cout << "Start random generation...." << endl;
+#endif
+  vector< complex<r_8> > fnOrig(Ntot);
+  int state=1234567 + 8912 ; //random seed
+  for(int i=0;i<Ntot;i++){
+    r_8 rv = drand(-1,1,&state);
+    r_8 iv = drand(-1,1,&state);
+    fnOrig[i] = complex<r_8>(rv,iv);
+  }
+#if DEBUG >=2 
+  std::copy(fnOrig.begin(), fnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
+#endif
+    
+  LaguerreTransform trans(N,R);
+
+  vector< complex<r_8> > fi(Ntot);
+#if DEBUG >= 1
+  cout << "Start MultiSynthesis...." << endl;
+#endif
+  trans.MultiSynthesis(fnOrig,fi,mapSize);
+    
+#if DEBUG >=2 
+  std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
+#endif
+
+  vector< complex<r_8> > fn(Ntot);
+  trans.MultiAnalysis(fi,fn,mapSize);
+
+  r_8 err_abs(0.);
+  r_8 err_rel(0.);
+  int imax;
+  for(int i=0;i<Ntot;i++){
+    if(i<10)cout << "("<<i<<"): " << fnOrig[i] << " <-> " << fi[i] << " <-> "  
+		 << fn[i] << endl;
+    complex<r_8> cdiff = (fnOrig[i] -  fn[i]) * conj(fnOrig[i] -  fn[i]);
+    r_16 diff = sqrt(cdiff.real());
+    if(diff>err_abs){ 
+      err_abs = diff;
+      imax = i;
+    }
+    complex<r_8> foriCAbs = fnOrig[i]*conj(fnOrig[i]);
+    r_8 foriAbs = sqrt(foriCAbs.real());
+    r_8 relatdiff = diff/foriAbs;
+	
+    if((relatdiff)>err_rel) err_rel = relatdiff;
+  }
+  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+  cout << "r_16 Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
+  cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl; 
+  
+
+}//MultiLaguerreTransformSdtAlone
+
+//---------------------------------------
+void MultiSphericalLaguerreTransform() {
+
+  int N = param.N;
+  int alpha = param.alpha;
+  int Lmax = param.Lmax;
+  r_8 R = param.R;
+
+
+  cout << " ___________  MultiSphericalLaguerreTransform   TEST _____________ " << endl;
+
+  int state = 1234567 + 8912 ; //random seed
+
+
+  LaguerreSphericalTransform sphlagtrans(Lmax,-1,N,R);
+      
+  sphlagtrans.SetThetaPhiMap();
+
+  int Nalm = sphlagtrans.Nalm();
+  int Ntot = N * Nalm; //total number of 3D-pixels & number of Coefficients of the Spherical Laguerre transform
+  int Npix =  sphlagtrans.Npix();
+  int NpTot= N * Npix; 
+      
+  vector< complex<r_8> > flmnOrig(Ntot);
+  int id=-1;
+  for(int n=0;n<N;n++){ //on each shell
+    for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
+      for(int l=m;l<Lmax;l++){
+	id++;
+	r_8 rv = drand(-1,1,&state);
+	r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
+	flmnOrig[id] = complex<r_8>(rv,iv);	    
+      }//end l-loop
+    }//end m-loop
+  }//end loop on shell
+
+#if DEBUG >=2 
+  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
+#endif
+      
+  //  Timer tm("Test 4 r_8 function");
+
+#if DEBUG >= 1
+  cout << "Main:Synthesis r_8 function..." << endl;
+#endif
+  vector<r_8> fijk(NpTot);
+  sphlagtrans.Synthesis(flmnOrig,fijk);
+#if DEBUG >= 2
+  for (int i=0; i<NpTot; i++){
+    cout << "fijk("<<i<<"): " << fijk(i) << endl;
+  }
+#endif
+  //  tm.Split("SHTL r_8/end synthesis");
+
+
+#if DEBUG >= 1
+  cout << "Main:Analysis r_8 function..." << endl;
+#endif
+  vector< complex<r_8> > flmn(Ntot);
+  sphlagtrans.Analysis(fijk,flmn);
+    
+  //  tm.Split("SHTL r_8/end analysis");
+
+
+  cout << "Error analysis..." << endl;
+
+  r_8 err_abs(0.);
+  r_8 err_rel(0.);
+  int imax;
+  for(int i=0;i<Ntot;i++){
+    if(i>(Ntot-9)) cout << "(" << i << ") : flnm Orig " <<  flmnOrig[i] << " <-> flmn Rec " <<  flmn[i] << endl;
+
+    complex<r_8> cdiff = (flmnOrig[i] -  flmn[i]) * conj(flmnOrig[i] -  flmn[i]);