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

(JEC/MR) 10/4/15 new code for Nodes/Weights computation which extend the...

(JEC/MR) 10/4/15 new code for Nodes/Weights computation which extend the validity wo needs of more 64-bits floating-points
parent f0d1bbd2
......@@ -69,6 +69,17 @@ $(EXE)lagsht_testsuite : $(OBJ)lagsht_testsuite.o $(CXXOBJ)
$(OBJ)lagsht_testsuite.o: lagsht_testsuite.cc $(CXXHDR)
echo "compile... $@"
$(CXXCOMPILE) -c $< -o $@
######################
.PHONY: laguerre
laguerre: $(EXE)laguerre
echo '---- laguerre made'
$(EXE)laguerre : $(OBJ)laguerre.o
echo "Link..."
$(CXXLINK) -o $@ $(OBJ)laguerre.o
$(OBJ)laguerre.o: laguerre.cc
echo "compile... $@"
$(CXXCOMPILE) -c $< -o $@
......
......@@ -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
# 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.
# 9/4/15
* new set of files for a Stand-a-lone version wo `Sophya`
* extension of the computation of the weights (Laguerre function) up to N ~ 5500 by transfert of the exponential scale a the very beginning (see laguerreBuilder.cc)
......
......@@ -3,7 +3,7 @@
#include <fstream>
#include <string.h>
#include <limits> // std::numeric_limits
#include <math.h> //log
//#include <math.h> //log
#include <string>
#include <typeinfo>
using namespace std;
......@@ -73,21 +73,24 @@ void LaguerreQuadrature() {
LaguerreFuncQuad lagFunc(N,alpha);
cout << "N = " << lagFunc.GetOrder() << " alpha = " << lagFunc.GetAlpha() << endl;
vector<r_16> nodes;
vector<r_16> weights;
vector<r_8> nodes;
vector<r_8> 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;
{
char buff1[80], buff2[80];
ofstream ofs1, ofs2;
sprintf(buff1,"lagNodes-%d-Func.txt",N);
sprintf(buff2,"lagWeights-%d-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 << weights[i] << endl;
}
}
}//LaguerreQuadrature
//--------------------------------------
......
......@@ -9,185 +9,17 @@
namespace LagSHT {
void LaguerrePolyQuad::QuadWeightNodes(vector<r_16>& nodes, vector<r_16>& weights, bool rescale){
nodes.resize(N_);
weights.resize(N_);
r_16 x;
r_16 tN = (r_16)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
r_16 tmp1 = tN + (r_16)alpha_;
int alphaM1 = alpha_-1;
for (int i=1;i<alphaM1;i++) tmp1 *= (tmp1-(r_16)1.);
tmp1 /= (tN+1.);
//Approximation by Sould & Secrest
for(int k=0;k<N_;k++){
if (k == 0) {
x = (r_16)((1.0+alpha_)*(3.+0.92*alpha_)/(1.+2.4*N_+1.8*alpha_));
} else if (k == 1) {
x += (r_16)((15.+6.25*alpha_)/(1.+0.9*alpha_+2.5*N_));
} else {//k>=2
int km1 = k-1;
r_16 c0 = (r_16)((1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1));
r_16 c1 = (r_16)(1.+0.3*alpha_);
x += c0*(x-nodes[k-2])/c1;
}
r_16 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
r_16 tmp0 = (!rescale) ? Value(x) : Value(x,expl(-x/2.)); //L_{N+1}(x_i) or L_{N+1}(x_i) Exp(-x_i/2)
tmp0 *= tmp0; //[L_{N+1}(x_i)]^2 or [L_{N+1}(x_i)]^2 Exp[-x_i]
N_--; //N+1-1 = N
tmp0 = 1./tmp0; // 1./[L_{N+1}(x_i)]^2 or Exp[x_i]/[L_{N+1}(x_i)]^2
weights[k] = tmp1 * tmp0 * x;
if(isnan( weights[k] )) cout << "LaguerrePolyQuad::QuadWeightNodes Pb weights["<<k<<"] nan: tmp0, tmp1, x " << tmp0 << ", " << tmp1 << ", " << x <<endl;
#if DEBUG >= 2
cout << " weight: " << weights(k) << endl;
#endif
}//loop k
}//QuadWeightNodes
r_16 LaguerrePolyQuad::Value(r_16 x, r_16 scale) {
r_16 Zero = (r_16)0.0;
r_16 One = (r_16)1.0;
if(N_<0) return Zero;
if(N_==0) return One*scale;
r_16 Two = (r_16)2.0;
r_16 alphaM1 = (r_16)alpha_ - One;
r_16 pnm1 = Zero;
r_16 pn = One*scale;
for(int i=1;i<=N_;i++){
r_16 ti = (r_16)i;
r_16 pnm2 = pnm1;
pnm1 = pn;
pn = ( (Two*ti + alphaM1 -x)*pnm1 - (ti + alphaM1)*pnm2 )/ti;
}
return pn;
}//Value
void LaguerrePolyQuad::Values(r_16 x, vector<r_16>& vals, r_16 scale) {
if(N_<0) throw LagSHTError("Values call with N<=0");
if(vals.size() != (N_+1) ) throw LagSHTError("Values with vals.size != (N+1)");
r_16 One = (r_16)1.0;
vals[0] = One*scale;
if(N_==0) return;
r_16 Zero = (r_16)0.0;
r_16 Two = (r_16)2.0;
r_16 alphaM1 = (r_16)alpha_ - One;
r_16 pnm1 = Zero;
r_16 pn = One*scale;
for(int i=1;i<=N_;i++){
r_16 ti = (r_16)i;
r_16 pnm2 = pnm1;
pnm1 = pn;
pn = ( (Two*ti + alphaM1 -x)*pnm1 - (ti + alphaM1)*pnm2 )/ti;
vals[i] = pn;
}
}//Values
r_16 LaguerrePolyQuad::LLpRatio(r_16 x) {
if (N_<=0) throw LagSHTError("LLpRatio call with N<=0");
r_16 pnm1 = (r_16)0.0;
r_16 pn = (r_16)1.0; //one can use here a scaling factor for numerical purpose
r_16 One = (r_16)1.0;
r_16 Two = (r_16)2.0;
r_16 talpha = (r_16)alpha_;
r_16 alphaM1 = talpha - One;
for(int i=1;i<=N_;i++){
r_16 ti = (r_16)i;
r_16 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)
r_16 tN = (r_16)N_;
r_16 ppn = (tN * pn - (tN+talpha)*pnm1)/x;
return pn/ppn;
}//LLpRatio
void LaguerrePolyQuad::RootUpdate(r_16& x) {
const int NLoops=100;
const r_16 newtonEnd= (r_16)(1.e-10);
r_16 xold = (r_16)0;
r_16 One = (r_16)1.0;
r_16 talpha = (r_16)(alpha_);
r_16 alphaP1 = One + talpha;
r_16 tN = (r_16)N_;
r_16 OneOv2 = (r_16)0.5;
//use simple Newton's method
for(int l=0; l<NLoops; l++){
r_16 LLpR = LLpRatio(x);
xold = x;
r_16 eps = LLpR;
x = xold - eps;
if(fabsl(x-xold)<newtonEnd) break;
}
}//RootUpdate
void LaguerreFuncQuad::QuadWeightNodes(vector<r_16>& nodes, vector<r_16>& weights) {
nodes.resize(N_);
weights.resize(N_);
lagPoly_.QuadWeightNodes(nodes,weights,true);
// vector<r_16> xnodes(N_);
// vector<r_16> xweights(N_);
// lagPoly_.QuadWeightNodes(xnodes,xweights,false);
for(int i=0;i<N_;i++){
// if (nodes[i] != xnodes[i]) cout << "LaguerreFuncQuad:nodes pb ("<<i<<")" << nodes[i] << ", " << xnodes[i] << endl;
// xweights[i] = (r_16)(expl(xnodes[i])*xweights[i]);
// if (fabsl(weights[i]-xweights[i])>1e-10*weights[i]) cout << "LaguerreFuncQuad:weights pb0 ("<<i<<")" << weights[i] << ", " << xweights[i] << ", " << weights[i]-xweights[i] << endl;
if(isnan(weights[i])) cout << "LaguerreFuncQuad:weights pb1: ("<<i<<")" << weights[i] << ", " << expl(nodes[i]/2.) << endl;
}
void LaguerreFuncQuad::QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) {
mod_laguerre_nodes_weights(N_,alpha_,nodes,weights);
}//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));
r_8 LaguerreFuncQuad::Value(r_8 x, r_8 scale){
return mod_laguerre (x,N_,alpha_);
}//LaguerreFuncQuad::Value
void LaguerreFuncQuad::Values(r_8 x, vector<r_16>& vals, r_16 scale) {
if(N_<0) throw LagSHTError("Values call with N<=0");
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)");
lagPoly_.Values((r_16)x,vals,scale*expl(-(r_16)x/2.0));
mod_laguerre_vector(x,N_,alpha_,vals);
}//LaguerreFuncQuad::Values
......
#ifndef LAGUERREBUILDER_SEEN
#define LAGUERREBUILDER_SEEN
#include <cmath> //very important !!!
#include <vector>
#include "lagsht_numbers.h"
......@@ -8,52 +9,8 @@ using namespace std;
namespace LagSHT {
////////////////////////////////////////////////////////////////////////////////////////////
//// --Class generalized Laguerre Polynomial and Gauss-Laguerre Quadrature -------------- //
////////////////////////////////////////////////////////////////////////////////////////////
//! Coumputation of the Zeros, Weights and values of Generalized Laguerre Polynomials (n, alpha)
class LaguerrePolyQuad {
public:
//! Creator
LaguerrePolyQuad(int n, int alpha=2): N_(n), alpha_(alpha) {}
//! Destructor
virtual ~LaguerrePolyQuad() {}
//! Getters
int GetOrder() { return N_; }
int GetAlpha() { return alpha_; }
//! Compute Weights and Nodes of the Gauss-Laguerre quadrature
/*!
\input scale : used to recale the weights
*/
void QuadWeightNodes(vector<r_16>& nodes, vector<r_16>& weights, bool rescale = false);
//! Compute L_n^alpha(x) * scale with n = N_
r_16 Value(r_16 x, r_16 scale= (r_16)1.0);
//! Compute L_n^alpha(x) * scale with n=0,...,N_ (included)
void Values(r_16 x, vector<r_16>& vals, r_16 scale= (r_16)1.0);
//! Compute L_n^alpha(x)/L'_n^alpha(x) with n = N_
r_16 LLpRatio(r_16 x);
protected:
//! Newton's Method to converge towards the true zero of the polynomail
void RootUpdate(r_16& x);
protected:
int N_; //!< Order of the polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial
}; //class LaguerrePolyQuad
/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
//// - -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//
......@@ -65,8 +22,7 @@ public:
//! Creator
LaguerreFuncQuad(int n, int alpha=2):
N_(n),
alpha_(alpha),
lagPoly_(n,alpha)
alpha_(alpha)
{}
//! Destructor
virtual ~LaguerreFuncQuad() {}
......@@ -76,20 +32,155 @@ public:
int GetAlpha() { return alpha_; }
//! Compute Weights and Nodes of the Gauss-Laguerre quadrature
void QuadWeightNodes(vector<r_16>& nodes, vector<r_16>& weights);
// void QuadWeightNodes(vector<r_16>& nodes, vector<r_16>& weights);
void QuadWeightNodes(vector<r_8>& nodes, vector<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) * 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);
//! Compute L_n^alpha(x) * scale with n=0,...,N_ (included)
void Values(r_8 x, vector<r_16>& vals, r_16 scale=1.0);
//! 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);
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;
const double scalefac=ldexp(1.,lim_exp);
const double xscalefac=ldexp(1.,-lim_exp);
double pnm1=0, pn=1;
for(int i=1;i<=N;i++)
{
double pnm2 = pnm1;
pnm1 = pn;
pn = ((2*i+alpha-1-x)*pnm1 - (i+alpha-1)*pnm2)/i;
// if necessary, scale down
// the function result is independent of the absolute scale
while (abs(pn)>=scalefac)
{ pn*=xscalefac; pnm1*=xscalefac; }
}
double ppn=(N*pn-(N+alpha)*pnm1)/x;
return -pn/ppn;
}
// returns L_N^alpha(x)*exp(-x/2)
double mod_laguerre (double x, int N, int alpha)
{
using namespace std;
if (N<0) return 0;
const int lim_exp=60;
const double scalefac=ldexp(1.,lim_exp);
const double xscalefac=ldexp(1.,-lim_exp);
const double log2=log(2.);
double pnm1=0, pn=1;
int scale=0;
for(int i=1;i<=N;i++)
{
double pnm2 = pnm1;
pnm1 = pn;
pn = ((2*i+alpha-1-x)*pnm1 - (i+alpha-1)*pnm2)/i;
while (abs(pn)>=scalefac) // if necessary, scale down and count rescalings
{ pn*=xscalefac; pnm1*=xscalefac; ++scale; }
}
// compute natural logarithm of the rescalings and add it to the Laguerre
// normalisation factor.
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.
void mod_laguerre_vector (double x, int N, int alpha, std::vector<double> &res)
{
using namespace std;
if (N<0) { res.resize(0); return; }
res.resize(N+1);
res[0]=exp(-0.5*x);
const int lim_exp=60;
const double scalefac=ldexp(1.,lim_exp);
const double xscalefac=ldexp(1.,-lim_exp);
const double log2=log(2.);
double pnm1=0, pn=1;
int scale=0;
for(int i=1;i<=N;i++)
{
double pnm2 = pnm1;
pnm1 = pn;
pn = ((2*i+alpha-1-x)*pnm1 - (i+alpha-1)*pnm2)/i;
while (abs(pn)>=scalefac) // if necessary, scale down and count rescalings
{ pn*=xscalefac; pnm1*=xscalefac; ++scale; }
// compute natural logarithm of the rescalings and add it to the Laguerre
// normalisation factor.
double natlog = -0.5*x + log2*(scale*lim_exp);
// apply all corrections
res[i] = pn*exp(natlog);
}
}
// computes nodes and weights for Gauss-Laguerre integration of order N
// NOTE: the weights are correct for modified Laguerre polynomials, i.e.
// L_N^alpha(x)*exp(-x/2)
void mod_laguerre_nodes_weights(int N, int alpha,
std::vector<double> &nodes, std::vector<double> &weights)
{
using namespace std;
nodes.resize(N);
weights.resize(N);
double tmp1 = N+alpha;
for (int i=1;i<alpha-1;i++) tmp1 *= (tmp1-1.);
tmp1 /= (N+1.);
double x=0.;
//Approximation by Sould & Secrest
for(int k=0; k<N; k++)
{
if (k==0)
x = (1.0+alpha)*(3.+0.92*alpha)/(1.+2.4*N+1.8*alpha);
else if (k==1)
x += (15.+6.25*alpha)/(1.+0.9*alpha+2.5*N);
else
{
int km1 = k-1;
double c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha/(1.+3.5*km1);
double c1 = 1.+0.3*alpha;
x += c0*(x-nodes[k-2])/c1;
}
// Refinement using Newton's method
// Once we notice that we are close, do one more step and then exit
bool done=false;
while (true)
{
double xold=x;
x+=laguerre_newtoncorr (x, N, alpha);
if (done) break;
if (abs(x-xold)<max(1e-12,1e-12*abs(x))) done=true;
}
nodes[k] = x;
double tmp0 = mod_laguerre(x, N+1, alpha);
weights[k] = x*tmp1/(tmp0*tmp0);
}
}
protected:
int N_; //!< Order of the polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial
LaguerrePolyQuad lagPoly_;
}; //LaguerreFuncQuad
......
......@@ -84,7 +84,8 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
LaguerreFuncQuad lag(N_-1);
vector<r_16> LnAll(N_); //all the values Lag_n(r) n:0...N-1
// vector<r_16> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
// Timer tm("lagTrans analysis Test");
......@@ -145,7 +146,8 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
//TODO : voir de meme si cette initialization sert
LaguerreFuncQuad lag(N_-1);
vector<r_16> LnAll(N_); //all the values Lag_n(r) n:0...N-1
// vector<r_16> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
......
......@@ -38,8 +38,10 @@ class LaguerreTransform : public BaseLaguerreTransform {
virtual ~LaguerreTransform() {}
//! Getters
vector<r_16> GetNodes() const { return nodes_; } //const ?
vector<r_16> GetWeights() const { return weights_; } //const ?
/* vector<r_16> GetNodes() const { return nodes_; } //const ? */
/* vector<r_16> GetWeights() const { return weights_; } //const ? */
vector<r_8> GetNodes() const { return nodes_; } //const ?
vector<r_8> GetWeights() const { return weights_; } //const ?
//! Multi Analysis at once (Values -> Coefficients)
......@@ -80,8 +82,10 @@ protected:
r_8 alphaFact_; //!< sqrt(alpha!)
r_8 R_; //!< maximal radius (not used yet)
vector<r_16> nodes_; //!< nodes of the Gauss-Laguerre quadrature used
vector<r_16> weights_; //!< weights of the Gauss-Laguerre quadrature used
/* vector<r_16> nodes_; //!< nodes of the Gauss-Laguerre quadrature used */
/* vector<r_16> weights_; //!< weights of the Gauss-Laguerre quadrature used */
vector<r_8> nodes_; //!< nodes of the Gauss-Laguerre quadrature used
vector<r_8> weights_; //!< weights of the Gauss-Laguerre quadrature used
};//LaguerreTransform
......
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