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

(JEC/MR) 11/4/15 new code to compute nodes and weights

parent 12a3ffc3
...@@ -7,6 +7,10 @@ git clone https://gitlab.in2p3.fr/campagne/LagSHT.git ...@@ -7,6 +7,10 @@ git clone https://gitlab.in2p3.fr/campagne/LagSHT.git
Download en Personne enregistree sur gitlab Download en Personne enregistree sur gitlab
git clone git@gitlab.in2p3.fr:campagne/LagSHT.git git clone git@gitlab.in2p3.fr:campagne/LagSHT.git
# 11/4/15
* update of the code by M.R to compute the nodes/weights included in the LaguerreBuilder file.
* NB: use of `<cmath>` mandatory
# 10/4/15 # 10/4/15
* new code from M. Reinecke to compute the nodes/weights of the generalized Laguerre function quadrature. It extends greatly the validity of the computation (N ~ 20000 tested) wo the needs of special floating-point representation. It is done by proprely scaling factor do on the fly. * new code from M. Reinecke to compute the nodes/weights of the generalized Laguerre function quadrature. It extends greatly the validity of the computation (N ~ 20000 tested) wo the needs of special floating-point representation. It is done by proprely scaling factor do on the fly.
......
...@@ -142,6 +142,9 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk, ...@@ -142,6 +142,9 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
vector< complex<r_8> > flmk(NtotCoeff); vector< complex<r_8> > flmk(NtotCoeff);
int off7n, off7k; int off7n, off7k;
tstack_push("SHT Analysis");
//Ylm analysis per shell //Ylm analysis per shell
for(int n =0; n<N_; n++){ for(int n =0; n<N_; n++){
off7n = n*Nalm_; off7n = n*Nalm_;
...@@ -159,6 +162,8 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk, ...@@ -159,6 +162,8 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
cout << "Start Laguerre..." << endl; cout << "Start Laguerre..." << endl;
#endif #endif
tstack_pop("SHT Analysis");
tstack_push("Laguerre MultiAnalysis");
//Laguerre Transform //Laguerre Transform
lagTrans_.MultiAnalysis(flmk, flmn, Nalm_); lagTrans_.MultiAnalysis(flmk, flmn, Nalm_);
...@@ -167,6 +172,9 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk, ...@@ -167,6 +172,9 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
cout << "End Laguerre..." << endl; cout << "End Laguerre..." << endl;
#endif #endif
tstack_pop("Laguerre MultiAnalysis");
#if DEBUG >= 2 #if DEBUG >= 2
for (int i=0;i<NtotCoeff; i++) { for (int i=0;i<NtotCoeff; i++) {
cout << "flmn("<<i<<"): " <<flmn[i] << endl; cout << "flmn("<<i<<"): " <<flmn[i] << endl;
......
#include <stdlib.h>
#include <iostream>
#include <math.h> //fabsl, expl
#include "lagsht_exceptions.h" #include "lagsht_exceptions.h"
#include "laguerreBuilder.h" #include "laguerreBuilder.h"
...@@ -9,17 +5,17 @@ ...@@ -9,17 +5,17 @@
namespace LagSHT { namespace LagSHT {
void LaguerreFuncQuad::QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) { void LaguerreFuncQuad::QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) const {
mod_laguerre_nodes_weights(N_,alpha_,nodes,weights); nodes_and_weights(N_,nodes,weights);
}//LaguerreFuncQuad::QuadWeightNodes }//LaguerreFuncQuad::QuadWeightNodes
r_8 LaguerreFuncQuad::Value(r_8 x, r_8 scale){ r_8 LaguerreFuncQuad::Value(r_8 x, r_8 scale) const{
return mod_laguerre (x,N_,alpha_); return value(x,N_);
}//LaguerreFuncQuad::Value }//LaguerreFuncQuad::Value
void LaguerreFuncQuad::Values(r_8 x, vector<r_8>& vals, r_8 scale) { 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() != (N_+1) ) throw LagSHTError("Values with vals.size != (N+1)");
mod_laguerre_vector(x,N_,alpha_,vals); values(x,N_,vals);
}//LaguerreFuncQuad::Values }//LaguerreFuncQuad::Values
......
#ifndef LAGUERREBUILDER_SEEN #ifndef LAGUERREBUILDER_SEEN
#define LAGUERREBUILDER_SEEN #define LAGUERREBUILDER_SEEN
#include <cmath> //very important !!! #include <cmath> //very important for abs definition
#include <vector> #include <vector>
#include <iostream>
#include "lagsht_numbers.h" #include "lagsht_numbers.h"
using namespace std; using namespace std;
...@@ -20,167 +22,177 @@ namespace LagSHT { ...@@ -20,167 +22,177 @@ namespace LagSHT {
class LaguerreFuncQuad { class LaguerreFuncQuad {
public: public:
//! Creator //! Creator
LaguerreFuncQuad(int n, int alpha=2): LaguerreFuncQuad(int N, int alpha=2): N_(N), alpha_(alpha),
N_(n), ln2(0.693147180559945309417232121458176),
alpha_(alpha) inv_ln2(1.4426950408889634073599246810018921),
{} fsmall(std::ldexp(1.,-large_exponent2)),
fbig(std::ldexp(1.,large_exponent2)),
xi(N+2), cf(maxscale+1-minscale) {
using namespace std;
xi[0]=0.;
for (size_t i=1; i<xi.size(); ++i)
xi[i]=1./i;
for (size_t i=0; i<cf.size(); ++i)
cf[i] = ldexp(1.,(int(i)+minscale)*large_exponent2);
}
//! Destructor //! Destructor
virtual ~LaguerreFuncQuad() {} virtual ~LaguerreFuncQuad() {}
//! Getters //! Getters
int GetOrder() { return N_; } int GetOrder() const { return N_; }
int GetAlpha() { return alpha_; } int GetAlpha() const { return alpha_; }
//! Compute Weights and Nodes of the Gauss-Laguerre quadrature //! Compute Weights and Nodes of the Gauss-Laguerre quadrature
// void QuadWeightNodes(vector<r_16>& nodes, vector<r_16>& weights); void QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) const;
void QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights);
//! Compute L_n^alpha(x) * exp(-x/2) * scale with n = N_ //! Compute L_n^alpha(x) * exp(-x/2) * scale with n = N_
// r_16 Value(r_8 x, r_16 scale=1.0); r_8 Value(r_8 x, r_8 scale=1.0) const;
r_8 Value(r_8 x, r_8 scale=1.0);
//! Compute L_n^alpha(x) * exp(-x/2) * scale with n=0,...,N_ (included) //! Compute L_n^alpha(x) * exp(-x/2) * scale with n=0,...,N_ (included)
// void Values(r_8 x, vector<r_16>& vals, r_16 scale=1.0);
void Values(r_8 x, vector<r_8>& vals, r_8 scale=1.0); void Values(r_8 x, vector<r_8>& vals, r_8 scale=1.0);
protected: protected:
// returns the amount by which x must be shifted to go closer to a root of
// L_N^alpha(x)
double laguerre_newtoncorr (double x, int N, int alpha)
{
using namespace std;
if(N<1) return 0;
const int lim_exp=60; /*! returns the amount by which x must be shifted to go closer to a root of
const double scalefac=ldexp(1.,lim_exp); L_N^alpha(x)
const double xscalefac=ldexp(1.,-lim_exp); */
double newtoncorr (double x, int N) const {
double pnm1=0, pn=1; using namespace std;
if(N<1) return 0.;
double pnm1=0., pn=1.;
double c1=alpha_-1., c2=alpha_-1.;
for(int i=1;i<=N;i++) for(int i=1;i<=N;i++)
{ {
double pnm2 = pnm1; double pnm2 = pnm1;
pnm1 = pn; pnm1 = pn;
pn = ((2*i+alpha-1-x)*pnm1 - (i+alpha-1)*pnm2)/i; c1+=2.;
// if necessary, scale down c2+=1.;
// the function result is independent of the absolute scale pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
while (abs(pn)>=scalefac) // if necessary, scale down
{ pn*=xscalefac; pnm1*=xscalefac; } // the function result is independent of the absolute scale
while (abs(pn)>=fbig)
{ pn*=fsmall; pnm1*=fsmall; }
} }
double ppn=(N*pn-(N+alpha)*pnm1)/x; double ppn=(N*pn-(N+alpha_)*pnm1)/x;
return -pn/ppn; return -pn/ppn;
} }//newtoncorr
// returns L_N^alpha(x)*exp(-x/2) /*! returns L_N^alpha(x)*exp(-x/2)
double mod_laguerre (double x, int N, int alpha) */
{ double value(double x, int N) const {
using namespace std; using namespace std;
if (N<0) return 0; if (N<0) return 0;
const int lim_exp=60; double log2expfac = -0.5*x*inv_ln2;
const double scalefac=ldexp(1.,lim_exp); int scale=int(log2expfac/large_exponent2);
const double xscalefac=ldexp(1.,-lim_exp); log2expfac-=scale*large_exponent2;
const double log2=log(2.);
double pnm1=0, pn=1; double pnm1=0., pn=exp(ln2*log2expfac);
int scale=0; double c1=alpha_-1., c2=alpha_-1.;
for(int i=1;i<=N;i++) for(int i=1;i<=N;i++)
{ {
double pnm2 = pnm1; double pnm2 = pnm1;
pnm1 = pn; pnm1 = pn;
pn = ((2*i+alpha-1-x)*pnm1 - (i+alpha-1)*pnm2)/i; c1+=2.;
while (abs(pn)>=scalefac) // if necessary, scale down and count rescalings c2+=1.;
{ pn*=xscalefac; pnm1*=xscalefac; ++scale; } pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
while (abs(pn)>=fbig) // if necessary, scale down and count rescalings
{ pn*=fsmall; pnm1*=fsmall; ++scale; }
} }
// compute natural logarithm of the rescalings and add it to the Laguerre return (scale<minscale) ? 0 : pn*cf[scale-minscale];
// normalisation factor. }//value
double natlog = -0.5*x + log2*(scale*lim_exp);
// apply all corrections
return pn*exp(natlog);
}
// returns L_n^alpha(x)*exp(-x/2) [n=0..N] in res. /*! returns L_n^alpha(x)*exp(-x/2) [n=0..N] in res.
void mod_laguerre_vector (double x, int N, int alpha, std::vector<double> &res) */
{ void values (double x, int N, std::vector<double> &res) {
using namespace std; using namespace std;
if (N<0) { res.resize(0); return; } if (N<0) { res.resize(0); return; }
res.resize(N+1); res.resize(N+1);
res[0]=exp(-0.5*x); res[0]=exp(-0.5*x);
const int lim_exp=60; double log2expfac = -0.5*x*inv_ln2;
const double scalefac=ldexp(1.,lim_exp); int scale=int(log2expfac/large_exponent2);
const double xscalefac=ldexp(1.,-lim_exp); log2expfac-=scale*large_exponent2;
const double log2=log(2.);
double pnm1=0, pn=1; double pnm1=0., pn=exp(ln2*log2expfac);
int scale=0; double c1=alpha_-1., c2=alpha_-1.;
for(int i=1;i<=N;i++) for(int i=1;i<=N;i++)
{ {
double pnm2 = pnm1; double pnm2 = pnm1;
pnm1 = pn; pnm1 = pn;
pn = ((2*i+alpha-1-x)*pnm1 - (i+alpha-1)*pnm2)/i; c1+=2.;
while (abs(pn)>=scalefac) // if necessary, scale down and count rescalings c2+=1.;
{ pn*=xscalefac; pnm1*=xscalefac; ++scale; } pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
// compute natural logarithm of the rescalings and add it to the Laguerre while (abs(pn)>=fbig) // if necessary, scale down and count rescalings
// normalisation factor. { pn*=fsmall; pnm1*=fsmall; ++scale; }
double natlog = -0.5*x + log2*(scale*lim_exp); res[i] = (scale<minscale) ? 0 : pn*cf[scale-minscale];
// apply all corrections
res[i] = pn*exp(natlog);
} }
} }//values
// computes nodes and weights for Gauss-Laguerre integration of order N /*! computes nodes and weights for Gauss-Laguerre integration of order N
// NOTE: the weights are correct for modified Laguerre polynomials, i.e. NOTE: the weights are correct for modified Laguerre polynomials, i.e.
// L_N^alpha(x)*exp(-x/2) L_N^alpha(x)*exp(-x/2)
void mod_laguerre_nodes_weights(int N, int alpha, */
std::vector<double> &nodes, std::vector<double> &weights) void nodes_and_weights(int N, std::vector<double> &nodes,
{ std::vector<double> &weights) const {
using namespace std; using namespace std;
nodes.resize(N); nodes.resize(N);
weights.resize(N); weights.resize(N);
double tmp1 = N+alpha; double tmp1 = N+alpha_;
for (int i=1;i<alpha-1;i++) tmp1 *= (tmp1-1.); for (int i=1;i<alpha_-1;i++) tmp1 *= (tmp1-1.);
tmp1 /= (N+1.); tmp1 /= (N+1.);
double x=0.; double x=0.;
//Approximation by Sould & Secrest // Approximation by Stroud & Secrest, Gaussian quadrature formulas, 1966
for(int k=0; k<N; k++) for(int k=0; k<N; k++)
{ {
if (k==0) if (k==0)
x = (1.0+alpha)*(3.+0.92*alpha)/(1.+2.4*N+1.8*alpha); x = (1.0+alpha_)*(3.+0.92*alpha_)/(1.+2.4*N+1.8*alpha_);
else if (k==1) else if (k==1)
x += (15.+6.25*alpha)/(1.+0.9*alpha+2.5*N); x += (15.+6.25*alpha_)/(1.+0.9*alpha_+2.5*N);
else else
{ {
int km1 = k-1; int km1 = k-1;
double c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha/(1.+3.5*km1); double c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1);
double c1 = 1.+0.3*alpha; double c1 = 1.+0.3*alpha_;
x += c0*(x-nodes[k-2])/c1; x += c0*(x-nodes[k-2])/c1;
} }
// Refinement using Newton's method // Refinement using Newton's method
// Once we notice that we are close, do one more step and then exit // Once we notice that we are close, do one more step and then exit
bool done=false; bool done=false;
while (true) while (true)
{ {
double xold=x; double xold=x;
x+=laguerre_newtoncorr (x, N, alpha); x+=newtoncorr (x, N);
if (done) break; if (done) break;
if (abs(x-xold)<max(1e-12,1e-12*abs(x))) done=true; if (abs(x-xold)<max(1e-12,1e-12*abs(x))) done=true;
} }
nodes[k] = x;
double tmp0 = mod_laguerre(x, N+1, alpha); nodes[k] = x;
weights[k] = x*tmp1/(tmp0*tmp0);
}
}
double tmp0 = value(x, N+1);
weights[k] = x*tmp1/(tmp0*tmp0);
}
}//nodes_and_weights
protected: protected:
int N_; //!< Order of the polynomial int N_; //!< Order of the polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial int alpha_; //!< second parameter of the generalized Laguerre polynomial
private:
enum { large_exponent2=90, minscale=-14, maxscale=14 };
double ln2, inv_ln2, fsmall, fbig;
std::vector<double> xi;
std::vector<double> cf;
}; //LaguerreFuncQuad }; //LaguerreFuncQuad
......
#include <math.h> #include <cmath>
#include <iostream> #include <iostream>
#include <algorithm> // for copy #include <algorithm> // for copy
#include <iterator> // for ostream_iterator #include <iterator> // for ostream_iterator
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment