(JEC/MR) 21/4/15 correct acces to libsharp map2alm and alm2map for compliance...

(JEC/MR) 21/4/15 correct acces to libsharp map2alm and alm2map for compliance to C++11; isolate ANNServer to a specific new file. Modify Makefiles accordingly

Darwin_g++_omp_make.inc

@@ -11,7 +11,7 @@ CFLAGS = -fno-common -g -O -fPIC
 #  Def compilateur C++ et flags
 CXX =  g++
 #JEC add -fopenmp -O3 -Wa,-q
-CXXFLAGS = -fno-common -g -O3 -fPIC -fopenmp  -Wa,-q
+CXXFLAGS = -fno-common -g -O3 -fPIC -fopenmp  -Wa,-q -Wall -std=c++11
 # flags specifiques pour templates repository...
 CXXTEMPFLG =
 #  Compilo fortran

Linux_g++_make.inc

@@ -13,7 +13,7 @@ CFLAGS = -Wall -Wpointer-arith -Wmissing-prototypes -O -g -fPIC
 CXX =  g++
 #  Flag de warning -Wsynth NON inclus par defaut (04/2007)
 ###CXXFLAGS = -Wall -Wpointer-arith -O -g -fPIC  CNFPHFLF 
-CXXFLAGS = -Wall -Wpointer-arith -O3 -g -fPIC -fopenmp  -fno-common 
+CXXFLAGS = -Wall -Wpointer-arith -O3 -g -fPIC -fopenmp  -fno-common  -Wall -std=c++11
 # flags specifiques pour templates repository...
 CXXTEMPFLG =
 #  Compilo fortran

Makefile

@@ -36,10 +36,10 @@ EXE = ./Objs/
 
 #  Define our target list
 .PHONY: default
-default: makedir lagsht_testsuite
+default: makedir lib lagsht_testsuite
 
 .PHONY: all 
-all : makedir lagsht_testsuite
+all : makedir lib lagsht_testsuite
 
 .PHONY: check
 check : makedir lagsht_testsuite
@@ -76,6 +76,8 @@ CXXSHOBJ = laguerreBuilder.o \
 CXXHDR = lagsht_exceptions.h \
 	lagsht_numbers.h \
 	lagsht_utils.h \
+	lagsht_geom.h \
+	lagsht_annserver.h \
 	laguerreBuilder.h \
 	laguerreTransform.h \
 	lagSphericTransform.h \

lagSphericTransform.cc

@@ -83,7 +83,7 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
 }//SetThetaPhiMap
 
 void BaseLagSphTransform::SetPixelPositions() {
-  sharp_geom_info *ginfo; 
+  const sharp_geom_info *ginfo; 
   ginfo = sharpjob_.get_geom_info();
   for (int i=0; i<ginfo->npairs; ++i) {
       double theta1 = ginfo->pair[i].r1.theta;
@@ -118,7 +118,7 @@ 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");
+  if(flmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
     
   fijk.resize(Ntot3DPix);
 
@@ -151,7 +151,8 @@ void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
     cout << "Synthesis... at shell num: " << n << endl;
 #endif
       
-    sharpjob_.alm2map(&(flmk.data()[off7n].real()),&fijk.data()[off7k],false);
+    //    sharpjob_.alm2map(&(flmk.data()[off7n].real()),&fijk.data()[off7k],false); //the false is to forbidd accumulation
+   sharpjob_.alm2map(reinterpret_cast<const r_8*>(&(flmk.data()[off7n])),&fijk.data()[off7k],false); //the false is to forbidd accumulation
 
   }//loop on shell
 
@@ -168,7 +169,7 @@ 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");
+  if(fijk.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error");
   flmn.resize(NtotCoeff);
 
 #if DEBUG >= 1
@@ -191,7 +192,11 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
     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
+    //    sharpjob_.map2alm(&(fijk.data()[off7k]),&(flmk.data()[off7n].real()),false); //the false is to forbidd accumulation
+
+    sharpjob_.map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(flmk.data()[off7n])),false); //the false is to forbidd accumulation
+
+
   }//loop on shell
 
 #if DEBUG >= 1

lagSphericTransform.h

@@ -10,10 +10,9 @@ using namespace std;
 #include "sharp_almhelpers.h"
 #include "sharp_cxx.h"
 
-//ANN
-#include <ANN/ANN.h>
 
 #include "laguerreTransform.h"
+#include "lagsht_annserver.h"
 #include "lagsht_exceptions.h"
 
 
@@ -24,85 +23,6 @@ namespace LagSHT {
 //// --------- Class Spherical Harmonic Laguerre Transform -- //
 ////////////////////////////////////////////////////////////////
 
-  typedef pair<r_8,r_8>    pixpos_t;
-  typedef vector<pixpos_t> geom_t;
-
-
-  class ANNServer {
-    
-  public:
-    //! Constructor
-    ANNServer(): nNear_(1), dim_(2), eps_(0), 
-		 npts_(0), dataPts_(NULL), kdTree_(NULL), nnIdx_(NULL), dists_(NULL) {
-      
-      nnIdx_ = new ANNidx[nNear_];  // allocate near neigh indices
-      dists_ = new ANNdist[nNear_]; // allocate near neighbor dists      
-    }
-
-    //! Constructor
-    ANNServer(geom_t vec):
-      nNear_(1), dim_(2), eps_(0),npts_(0), dataPts_(NULL), kdTree_(NULL) {
-      nnIdx_ = new ANNidx[nNear_];  // allocate near neigh indices
-      dists_ = new ANNdist[nNear_]; // allocate near neighbor dists      
-    }//Ctor
-
-    //! Set kdtree
-    void SetkdTree(geom_t vec) {
-
-      npts_ = vec.size();
-      dataPts_ = annAllocPts(npts_, dim_);
-      if(dataPts_ == NULL) throw LagSHTError("ANN Failure dataPts alloc");
-
-      for(int i=0;i<npts_;i++){
-	ANNpoint p = dataPts_[i];
-	p[0] = vec[i].first;
-	p[1] = vec[i].second;
-      }
-      kdTree_ = new ANNkd_tree(dataPts_,npts_,dim_);
-      if(kdTree_ == NULL) throw LagSHTError("ANN Failure kdTree alloc");
-      
-    }//SetkdTree
-
-    //! Destructor
-    virtual ~ANNServer() {
-      if (dataPts_) annDeallocPts(dataPts_);
-      if (nnIdx_) delete [] nnIdx_;	
-      if (dists_) delete [] dists_;
-      if (kdTree_) delete kdTree_;      
-      annClose();
-    }//Dtor
-
-    int GetClosestId(r_8 theta, r_8 phi) const {
-
-      ANNpoint queryPt;
-      queryPt = annAllocPt(dim_);
-      if(queryPt == NULL) throw LagSHTError("Get Failure queryPt alloc");
-      queryPt[0] = theta;
-      queryPt[1] = phi;
-      
-      if (kdTree_ == NULL) throw LagSHTError("Get Failure kdTree alloc");
-      if (nnIdx_ == NULL) throw LagSHTError("Get Failure nnIdx alloc");
-      if (dists_ == NULL) throw LagSHTError("Get Failure dists alloc");
-      
-      kdTree_->annkSearch(queryPt,nNear_,nnIdx_,dists_,eps_);
-      annDeallocPt(queryPt);
-      
-      return nnIdx_[0]; //first closest
-    }//GetClosestId
-
-
-  protected:
-    int nNear_; //nomber of solution to the closest point
-    int dim_;   //dimsion of the point space
-    r_8 eps_;   // error bound
-    int npts_;  //number of points to fill the kdTree
-
-    ANNkd_tree* kdTree_;  //the tree
-    ANNidxArray  nnIdx_; // near neighbor indices
-    ANNdistArray dists_; // near neighbor distances
-    ANNpointArray dataPts_; //the initial grid points
-
-  };//ANNServer
 
 class BaseLagSphTransform {
  public:

lagsht_annserver.h

+#ifndef LAGSHTANNSERVER_SEEN
+#define LAGSHTANNSERVER_SEEN
+
+//ANN
+#include <ANN/ANN.h>
+
+#include "lagsht_geom.h"
+#include "lagsht_exceptions.h"
+
+#define DEBUB 0
+
+namespace LagSHT {
+
+////////////////////////////////////////////////////////////////////////////
+//// --------- Class gateway to Approximate Nearest Neighbor Searching -- //
+////////////////////////////////////////////////////////////////////////////
+
+class ANNServer {
+    
+public:
+  //! Constructor default
+  ANNServer(): nNear_(1), dim_(2), eps_(0), 
+	       npts_(0), kdTree_(NULL), nnIdx_(NULL), 
+	       dists_(NULL), dataPts_(NULL) {
+      
+    nnIdx_ = new ANNidx[nNear_];  // near neighboors indices
+    dists_ = new ANNdist[nNear_]; //  --------------- dists      
+  }
+
+  //! Constructor
+  ANNServer(geom_t vec):
+    nNear_(1), dim_(2), eps_(0),npts_(0), kdTree_(NULL),  dataPts_(NULL) {
+    nnIdx_ = new ANNidx[nNear_];
+    dists_ = new ANNdist[nNear_];
+  }//Ctor
+
+  //! Set kdtree
+  void SetkdTree(geom_t vec) {
+
+    npts_ = vec.size();
+    dataPts_ = annAllocPts(npts_, dim_);
+    if(dataPts_ == NULL) throw LagSHTError("ANN Failure dataPts alloc");
+
+    for(int i=0;i<npts_;i++){
+      ANNpoint p = dataPts_[i];
+      p[0] = vec[i].first;
+      p[1] = vec[i].second;
+    }
+    kdTree_ = new ANNkd_tree(dataPts_,npts_,dim_);
+    if(kdTree_ == NULL) throw LagSHTError("ANN Failure kdTree alloc");
+      
+  }//SetkdTree
+
+  //! Destructor
+  virtual ~ANNServer() {
+    if (dataPts_) annDeallocPts(dataPts_);
+    if (nnIdx_) delete [] nnIdx_;	
+    if (dists_) delete [] dists_;
+    if (kdTree_) delete kdTree_;      
+    annClose();
+  }//Dtor
+
+  //! Get the index of the closet point (theta, phi)
+  int GetClosestId(r_8 theta, r_8 phi) const {
+
+    ANNpoint queryPt;
+    queryPt = annAllocPt(dim_);
+
+#if DEBUG >=1
+    if(queryPt == NULL) throw LagSHTError("Get Failure queryPt alloc");
+#endif
+
+    queryPt[0] = theta;
+    queryPt[1] = phi;
+      
+#if DEBUG >=1
+    if (kdTree_ == NULL) throw LagSHTError("Get Failure kdTree alloc");
+    if (nnIdx_ == NULL) throw LagSHTError("Get Failure nnIdx alloc");
+    if (dists_ == NULL) throw LagSHTError("Get Failure dists alloc");
+#endif
+      
+    kdTree_->annkSearch(queryPt,nNear_,nnIdx_,dists_,eps_);
+    annDeallocPt(queryPt);
+      
+    return nnIdx_[0]; //first closest
+  }//GetClosestId
+
+
+protected:
+  int nNear_; //<! nomber of solution to the closest point
+  int dim_;   //<! dimsion of the point space
+  r_8 eps_;   //<! error bound
+  int npts_;  //<! number of points to fill the kdTree
+
+  ANNkd_tree* kdTree_;  //<! the tree
+  ANNidxArray  nnIdx_; //<! near neighbor indices
+  ANNdistArray dists_; //<! near neighbor distances
+  ANNpointArray dataPts_; //<! the initial grid points
+
+};//ANNServer
+
+}//end namespace
+
+#endif //LAGSHTANNSERVER

lagsht_geom.h

+#ifndef LAGSHTGEOM_SEEN
+#define LAGSHTGEOM_SEEN
+
+namespace LagSHT {
+
+typedef pair<r_8,r_8>    pixpos_t;
+typedef vector<pixpos_t> geom_t;
+
+}//end namespace
+
+#endif //LAGSHTGEOM

lagsht_testsuite.cc

@@ -40,9 +40,10 @@ static double drand (double min, double max, int *state){
 //----------------------------------------------------
 void TestBasic() {
 
-  int N = param.N;
-  int alpha = param.alpha;
-  int Lmax = param.Lmax;
+  // to access some general parameters... use param.
+  // int N = param.N;
+  // int alpha = param.alpha;
+  // int Lmax = param.Lmax;
 
 
 
@@ -71,7 +72,6 @@ void LaguerreQuadrature() {
 
   int N = param.N;
   int alpha = param.alpha;
-  int Lmax = param.Lmax;
 
   cout << " ___________  LaguerreQuadrature  TEST _____________ " << endl;
 
@@ -100,8 +100,6 @@ void LaguerreQuadrature() {
 //--------------------------------------
 void MultiLaguerreTransformSdtAlone() {
   int N = param.N;
-  int alpha = param.alpha;
-  int Lmax = param.Lmax;
   r_8 R  = param.R;
 
   cout << " ___________   MultiLaguerreTransformSdtAlone  TEST _____________ " << endl;
@@ -143,7 +141,7 @@ void MultiLaguerreTransformSdtAlone() {
 
   r_8 err_abs(0.);
   r_8 err_rel(0.);
-  int imax;
+  int imax=-1;
   for(int i=0;i<Ntot;i++){
     if(i<10)cout << "("<<i<<"): " << fnOrig[i] << " <-> " << fi[i] << " <-> "  
 		 << fn[i] << endl;
@@ -170,7 +168,6 @@ void MultiLaguerreTransformSdtAlone() {
 void MultiSphericalLaguerreTransform() {
 
   int N = param.N;
-  int alpha = param.alpha;
   int Lmax = param.Lmax;
   r_8 R = param.R;
   string geometry = param.geometry;
@@ -245,7 +242,7 @@ void MultiSphericalLaguerreTransform() {
 
   r_8 err_abs(0.);
   r_8 err_rel(0.);
-  int imax;
+  int imax = -1;
   for(int i=0;i<Ntot;i++){
     if(i>(Ntot-9)) cout << "(" << i << ") : flnm Orig " <<  flmnOrig[i] << " <-> flmn Rec " <<  flmn[i] << endl;
 
@@ -276,9 +273,7 @@ void MultiSphericalLaguerreTransform() {
 
 int main(int narg, char *arg[]) {
 
-  int debug=0; 
-
-   size_t maxmemsize =  getMemorySize()/1e6;
+   unsigned int maxmemsize =  getMemorySize()/1e6;
    cout << "Max Memory size: " <<  maxmemsize << " MBytes" <<  endl;
 
   int N = 5; 
@@ -302,10 +297,6 @@ int main(int narg, char *arg[]) {
 	   << endl; 
       return 0;
     }
-    else if (strcmp(arg[ka],"-debug")==0) {
-      debug=atoi(arg[ka+1]);
-      ka+=2;
-    }
     else if (strcmp(arg[ka],"-n")==0) {
       N=atoi(arg[ka+1]);
       ka+=2;
@@ -347,7 +338,7 @@ int main(int narg, char *arg[]) {
     if(N>5500) throw LagSHTError("WARNING: N>5500 the algorithm may fail due to weight estimates"); 
 
 
-    size_t estimateMem = 8*N*(uint_8)((2*Lmax+1)*(2*Lmax+1))/1e6; //Npix * 8Bytes and then in MBytes
+    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;

laguerreBuilder.cc

@@ -14,7 +14,7 @@ r_8 LaguerreFuncQuad::Value(r_8 x, r_8 scale) const{
 }//LaguerreFuncQuad::Value
 
 void LaguerreFuncQuad::Values(r_8 x, vector<r_8>& vals, r_8 scale) {
-  if(vals.size() != (N_+1) ) throw LagSHTError("Values with vals.size != (N+1)");
+  if(vals.size() != (size_t)(N_+1) ) throw LagSHTError("Values with vals.size != (N+1)");
   values(x,N_,vals);
 }//LaguerreFuncQuad::Values
 

laguerreTransform.cc

@@ -85,7 +85,7 @@ int LaguerreTransform::R2Index(r_8 r) const {
 void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi, 
 				      vector< complex<r_8> >& fn, int stride) {
 
-    if(fi.size() != N_*stride) throw LagSHTError("LagTransform::Analysis size error");
+  if(fi.size() != (size_t)(N_*stride) ) throw LagSHTError("LagTransform::Analysis size error");
     fn.resize(N_*stride);
 
 #if DEBUG >= 1
@@ -147,7 +147,7 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
 void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn, 
 				       vector< complex<r_8> >& fi, int stride) {
 
-    if(fn.size() != N_*stride) throw LagSHTError("LagTransform::Synthesis size error");
+  if(fn.size() != (size_t)(N_*stride) ) throw LagSHTError("LagTransform::Synthesis size error");
     fi.resize(N_*stride);
 
 #if DEBUG >= 1
@@ -156,8 +156,6 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
 #if DEBUG >= 2
     std::copy(fn.begin(), fn.end(), std::ostream_iterator<char>(std::cout, " "));
 #endif
-    int off7k, off7n; 
-    int threadId;
     r_8 invalphaFact = 1.0/alphaFact_; //1/alpha! 
 
     vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif