Commit bd7c2742 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 3/4/15 add specialization of template<r_8> that explicitly include a...

(JEC) 3/4/15 add specialization of template<r_8> that explicitly include a mixed r_16/r_8 algorithm and use generalized Laguerre Functions.
parent a472536c
......@@ -7,6 +7,9 @@ git clone https://gitlab.in2p3.fr/campagne/LagSHT.git
Download en Personne enregistree sur gitlab
git clone git@gitlab.in2p3.fr:campagne/LagSHT.git
# Version du 03/04/15
* refonte des algorithmes de Synthese/Analyse de la transformation de Laguerre en utilisant des functions generalisees de Laguerre et non les polynomes generalises. Ainsi, le facteur exponentiel permet de manipuler au mieux des nombres 64-bits sauf a quelques endroits ou on a besoin de 128-bits. L'ineterface avec `libsharp` se fait donc mieux.
* L'algorithme d'Analyse la transformation Harmonqiue Spherique + Laguerre a ete revu pour tirer profit des calculs iteratifs des polynomes de Laguerre comme dans le cza de la Synthese. Cela induit une acceleration significative.
# Version du 25/3/15
......
......@@ -13,8 +13,7 @@
namespace SOPHYA {
template <class T>
void LaguerreSphericalTransform<T>::SetThetaPhiMap(string choice){
void BaseLagSphTransform::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
......@@ -111,6 +110,9 @@ void LaguerreSphericalTransform<T>::Synthesis(const NDataBlock< complex<T> >& fl
#endif
}//Synthesis
template <class T>
void LaguerreSphericalTransform<T>::Analysis(const NDataBlock<T>& fijk,
NDataBlock< complex<T> >& flmn) {
......@@ -184,6 +186,8 @@ void LaguerreSphericalTransform<T>::Analysis(const NDataBlock<T>& fijk,
#endif
}//Analysis
///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
......@@ -205,5 +209,119 @@ template class LaguerreSphericalTransform<xdouble>;
#endif
#endif
template <>
LaguerreSphericalTransform<r_8>::LaguerreSphericalTransform(int Lmax, int Mmax, int Nmax, r_8 R, int alpha):
BaseLagSphTransform(Lmax,Mmax,Nmax,alpha), R_(R), lagTrans_(Nmax,R,alpha) {
}//Ctor
template <>
LaguerreSphericalTransform<r_8>::~LaguerreSphericalTransform() {}
template <>
void LaguerreSphericalTransform<r_8>::Synthesis(const NDataBlock< complex<r_8> >& flmn,
NDataBlock<r_8>& 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<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 >= 1
cout << "Synthesis... at shell num: " << n << " scale = " << scale << 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
template <>
void LaguerreSphericalTransform<r_8>::Analysis(const NDataBlock<r_8>& fijk,
NDataBlock< complex<r_8> >& 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<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()[0]),&(flmk.Data()[0].real()),false); //the false is to forbidd accumulation
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
......@@ -24,21 +24,17 @@ namespace SOPHYA {
//// --------- Class Spherical Harmonic Laguerre Transform -- //
////////////////////////////////////////////////////////////////
template<class T>
class LaguerreSphericalTransform {
class BaseLagSphTransform {
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),
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;
}//Ctor
if (Mmax<0) M_=Lmax;
}
//! Destructor
virtual ~LaguerreSphericalTransform() {}
virtual ~BaseLagSphTransform() {}
//! Number of colatitude rings of the 2D-sphere representation
int Nrings() { return Nrings_; }
//! Number of colatitude rings of the 2D-sphere representation
......@@ -50,26 +46,9 @@ class LaguerreSphericalTransform {
//! 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
......@@ -94,11 +73,8 @@ 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
......@@ -106,6 +82,93 @@ protected:
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 !!!
};
template<class T>
class LaguerreSphericalTransform : public BaseLagSphTransform {
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
LaguerreSphericalTransform(int Lmax, int Mmax, int Nmax, r_8 R, int alpha=2):
BaseLagSphTransform(Lmax,Mmax,Nmax,alpha), R_(R), lagTrans_(Nmax,(T)R,alpha) {
}//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
......@@ -113,6 +176,8 @@ protected:
//--------- extern template declarations (if needed) -----------
#if defined ( NEED_EXT_DECL_TEMP ) && !defined( LAGSPHERICTRANSFORM_CC_BFILE )
extern template class LaguerreSphericalTransform<r_8>;
#ifdef SO_LDBLE128
extern template class LaguerreSphericalTransform<r_16>;
#endif
......
......@@ -187,4 +187,29 @@ template class LaguerrePolyQuad<xdouble>;
#endif
#endif
void LaguerreFuncQuad::QuadWeightNodes(NDataBlock<r_8>& nodes, NDataBlock<r_8>& weights) {
NDataBlock<r_16> xnodes(N_);
NDataBlock<r_16> xweights(N_);
lagPoly_.QuadWeightNodes(xnodes,xweights);
nodes.ReSize(N_);
weights.ReSize(N_);
for(int i=0;i<N_;i++){
nodes(i) = (r_8)(xnodes(i));
weights(i) = (r_8)(expl(xnodes(i))*xweights(i));
}
}//LaguerreFuncQuad::QuadWeightNodes
r_16 LaguerreFuncQuad::Value(r_8 x, r_16 scale){
return lagPoly_.Value((r_16)x,scale*expl(-(r_16)x/2.0));
}//LaguerreFuncQuad::Value
void LaguerreFuncQuad::Values(r_8 x, NDataBlock<r_16>& vals, r_16 scale) {
if(N_<0) throw MathExc("Values call with N<=0");
if(vals.Size() != (N_+1) ) throw MathExc("Values with vals.size != (N+1)");
lagPoly_.Values((r_16)x,vals,scale*expl(-(r_16)x/2.0));
}//LaguerreFuncQuad::Values
}// Fin de namespace
......@@ -15,9 +15,10 @@ using namespace NTL;
namespace SOPHYA {
////////////////////////////////////////////////////////////////
//// ------------------- Class LaguerrePolyQuad -------------- //
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////
//// --Class generalized Laguerre Polynomial and Gauss-Laguerre Quadrature -------------- //
////////////////////////////////////////////////////////////////////////////////////////////
// T: NTLxdouble: tested
// r_16 : tested
......@@ -71,6 +72,47 @@ extern template class LaguerrePolyQuad<xdouble>;
#endif // Fin de if defined ( NEED_EXT_DECL_TEMP )
/////////////////////////////////////////////////////////////////////////////////////////////
//// - -Class generalized Laguerre Function and Gauss-Laguerre Quadrature -------------- //
//// The Laguerre Function is scaled by exp(-r/2) compared to the Laguerre Polynomials //
//// the exp(r_i) scales the weights compared to the Polynomial Gauss-Laguerre Quadrature//
//// (nb. an public heritance from LaguerrePolyQuad would made confusion) //
/////////////////////////////////////////////////////////////////////////////////////////////
class LaguerreFuncQuad {
public:
//! Creator
LaguerreFuncQuad(int n, int alpha=2):
N_(n),
alpha_(alpha),
lagPoly_(n,alpha)
{}
//! Destructor
virtual ~LaguerreFuncQuad() {}
//! Getters
int GetOrder() { return N_; }
int GetAlpha() { return alpha_; }
//! Compute Weights and Nodes of the Gauss-Laguerre quadrature
void QuadWeightNodes(NDataBlock<r_8>& nodes, NDataBlock<r_8>& weights);
//! Compute L_n^alpha(x) * scale with n = N_
r_16 Value(r_8 x, r_16 scale=1.0);
//! Compute L_n^alpha(x) * scale with n=0,...,N_ (included)
void Values(r_8 x, NDataBlock<r_16>& vals, r_16 scale=1.0);
protected:
int N_; //!< Order of the polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial
LaguerrePolyQuad<r_16> lagPoly_;
}; //LaguerreFuncQuad
} // Fin du namespace
......
......@@ -8,13 +8,13 @@
#define LAGUERRETRANSFORM_CC_BFILE // avoid extern template declarations
#include "laguerreTransform.h"
#define DEBUG 0
#define DEBUG 1
namespace SOPHYA {
template <class T>
LaguerreTransform<T>::LaguerreTransform(int N, T R, int alpha) :
N_(N), alpha_(alpha), R_(R) {
BaseLaguerreTransform(N,alpha), R_(R) {
if(alpha<0) throw MathExc("LagTransform call with alpha<0");
LaguerrePolyQuad<T> lag(N_,alpha_);
......@@ -35,7 +35,7 @@ LaguerreTransform<T>::LaguerreTransform(int N, T R, int alpha) :
#if DEBUG >= 1
cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
cout << "Xdble Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
#endif
if (fabs(sumR-sumRTh)>(0.000000001*sumRTh)) {
......@@ -65,37 +65,71 @@ void LaguerreTransform<T>::MultiAnalysis(const NDataBlock< complex<T> >& fi, NDa
fn.ReSize(N_*stride);
#if DEBUG >= 1
cout << "Analysis...." << endl;
cout << "Multi Analysis...." << endl;
fi.Print(cout);
#endif
T invalphaFact = ((T)1.0)/alphaFact_; //1/alpha!
// T factn = ((T)1.0)/alphaFact_;
NDataBlock<T> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
facts(0) = invalphaFact;
for(int n=1;n<N_;n++) facts(n) = facts(n-1)*sqrt((T)( ((T)n)/((T)(n+alpha_)) ) );
// a voir si on peut interchanger la loop sur n et i mais alors la loop sur les "l" changerait aussi...
// au pire c'est une matrice NxN des facteurs Ln(r_i)= Matrice(n,i) qui serait commune a la synthese et a l'analyse...
LaguerrePolyQuad<T> lag(N_-1);
NDataBlock<T> LnAll(N_); //all the values Lag_n(r) n:0...N-1
int off7i, off7n;
for (int n=0; n<N_;n++){
LaguerrePolyQuad<T> lag(n);
off7n = n*stride;
for (int i = 0; i<N_; i++ ){
T r = nodes_(i);
for (int i = 0; i<N_; i++ ){
off7i = i*stride;
T r = nodes_(i);
lag.Values(r,LnAll);
for (int n=0; n<N_;n++){
off7n = n*stride;
//Lni = sqrt(n!/(n+alpha)!) L_n^(alpha)(r_i)
T Lni = lag.Value(r,facts(n));
off7i = i*stride;
//#pragma omp parallel for schedule(dynamic,1)
T Lni = LnAll(n)*facts(n);
//#pragma omp parallel for
for (int l=0; l<stride; l++) {
fn(l+off7n) += weights_(i) * fi(l+off7i) * Lni;
}
}//loop on i
}//loop on n
}//MultiAnalysis
//Version non optimisee au niveau des loop i et n: on peut les inverser
// template <class T>
// void LaguerreTransform<T>::MultiAnalysis(const NDataBlock< complex<T> >& fi, NDataBlock< complex<T> >& fn, int stride) {
// if(fi.Size() != N_*stride) throw SzMismatchError("LagTransform::Analysis size error");
// fn.ReSize(N_*stride);
// #if DEBUG >= 1
// cout << "Multi Analysis...." << endl;
// fi.Print(cout);
// #endif
// T invalphaFact = ((T)1.0)/alphaFact_; //1/alpha!
// NDataBlock<T> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
// facts(0) = invalphaFact;
// for(int n=1;n<N_;n++) facts(n) = facts(n-1)*sqrt((T)( ((T)n)/((T)(n+alpha_)) ) );
// int off7i, off7n;
// for (int n=0; n<N_;n++){
// LaguerrePolyQuad<T> lag(n);
// off7n = n*stride;
// for (int i = 0; i<N_; i++ ){
// T r = nodes_(i);
// //Lni = sqrt(n!/(n+alpha)!) L_n^(alpha)(r_i)
// T Lni = lag.Value(r,facts(n));
// off7i = i*stride;
// //#pragma omp parallel for
// for (int l=0; l<stride; l++) {
// fn(l+off7n) += weights_(i) * fi(l+off7i) * Lni;
// }
// }//loop on i
// }//loop on n
// }//MultiAnalysis
template <class T>
void LaguerreTransform<T>::MultiSynthesis(const NDataBlock< complex<T> >& fn, NDataBlock< complex<T> >& fi, int stride) {
if(fn.Size() != N_*stride) throw SzMismatchError("LagTransform::Synthesis size error");
......@@ -128,7 +162,7 @@ void LaguerreTransform<T>::MultiSynthesis(const NDataBlock< complex<T> >& fn, ND
//flm(ri) = Sum_{n=0}^{N-1} flmn * sqrt(n!/(n+alpha)!) L_n^(alpha)(r_i)
//(lm) is represented by a single index l
off7n = n*stride;
//#pragma omp parallel for schedule(dynamic,1)
//#pragma omp parallel for
for (int l=0; l<stride; l++) {
fi(l+off7i) += Lni * fn(l+off7n);
}
......@@ -159,6 +193,146 @@ template class LaguerreTransform<xdouble>;
#endif
template<>
LaguerreTransform<r_8>::LaguerreTransform(int N, r_8 R, int alpha) :
BaseLaguerreTransform(N,alpha), R_(R) {
cout << "LaguerreTransform<r_8> start...." <<endl;
if(alpha<0) throw MathExc("LagTransform call with alpha<0");
LaguerreFuncQuad lag(N_,alpha_);
lag.QuadWeightNodes(nodes_,weights_);
//verification
r_8 sumR = 0.;
r_8 sumW = 0.;
r_8 tau = nodes_(N_-1)/708.0; //double number 708 ~ ln(minimum positive double)
for (int i=0;i<N_;i++){
sumR += nodes_(i);
sumW += weights_(i)*exp(-nodes_(i)/tau);
}
int alphaFact = 1;
for(int i=1;i<=alpha_; i++) alphaFact *= i;
r_8 sumRTh = (r_8)(N_*(N_+alpha_));
r_8 sumWTh = (r_8)(alphaFact)*pow(tau,(r_8)(alpha_+1));
#if DEBUG >= 1
cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
#endif
if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< "WARNING!!!! verify LaguerrePolyQuad sum of nodes not accurate: "
<< " current value " << sumR << " theoretical value : " << sumRTh
<< " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< endl;
}
if (fabs(sumW-sumWTh)>(0.000001*sumWTh)) {
cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< "WARNING!!!! verify LaguerrePolyQuad sum of weights not accurate: "
<< " current value " << sumW << " theoretical value : " << sumRTh
<< " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"
<< endl;
}
alphaFact_ = sqrt((r_8)alphaFact); //sqrt(alpha!)
}//Ctor
template <>
void LaguerreTransform<r_8>::MultiAnalysis(const NDataBlock< complex<r_8> >& fi, NDataBlock< complex<r_8> >& fn, int stride) {
if(fi.Size() != N_*stride) throw SzMismatchError("LagTransform::Analysis size error");
fn.ReSize(N_*stride);
#if DEBUG >= 1
cout << "Analysis with Function...." << endl;
fi.Print(cout);
#endif
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
NDataBlock<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);
NDataBlock<r_16> LnAll(N_); //all the values Lag_n(r) n:0...N-1
int off7i, off7n;
//Try i-loop before n-loop (the reverse is the natural order)
for (int i = 0; i<N_; i++ ){
off7i = i*stride;
r_8 r = nodes_(i);
lag.Values(r,LnAll);
for (int n=0; n<N_;n++){
// LaguerreFuncQuad lag(n);