Revision 795 branches/replacement/src/linad99/gaussher.cpp
gaussher.cpp (revision 795)  

4  4 
* Author: David Fournier 
5  5 
* Copyright (c) 2009 ADMB Foundation 
6  6 
*/ 
7 
#include <fvar.hpp> 

7 
/** 

8 
* \file 

9 
* This file has routines used for 

10 
* GaussHermite Quadrature and 

11 
* GaussLegendre Quadrature. 

12 
*/ 

8  13  
14 
#include <fvar.hpp> 

9  15 
static double eps=3.0e14; 
10  16 
static double pim=0.7511255444649427; 
11  17 
static int maxit=50; 
12  18  
13  
14  19 
void imtqlx ( const dvector& _d, const dvector& _e, const dvector& _z ); 
15  20  
16  21 
double sign ( double x ) 
...  ...  
28  33 
return value; 
29  34 
} 
30  35  
31  
32 
/** GaussHermite quadature. 

33  
34 
\n\n The implementation of this algorithm was inspired by 

35 
"Numerical Recipes in C", 2nd edition, 

36 
Press, Teukolsky, Vetterling, Flannery, chapter 4 

37  
38 
\deprecated Scheduled for replacement by 2010. 

36 
/** 

37 
* GaussHermite quadature. 

38 
* Computes a GaussHermite quadrature formula with simple knots. 

39 
* \param _t array of abscissa 

40 
* \param _wts array of corresponding wights 

39  41 
*/ 
40  
41  42 
void gauss_hermite (const dvector& _t,const dvector& _wts) 
42  
43 
//****************************************************************************80 

44  43 
// 
45  44 
// Purpose: 
46  45 
// 
47 
// CDGQF computes a Gauss quadrature formula with default A, B and simple knots.


46 
// computes a Gauss quadrature formula with simple knots.


48  47 
// 
49  48 
// Discussion: 
50  49 
// 
...  ...  
113  112 
bj(i) = sqrt((ilb+1)/2.0); 
114  113 
} 
115  114  
116  
117  115 
// Compute the knots and weights. 
118  
119  116 
if ( zemu <= 0.0 ) // Exit if the zeroth moment is not positive. 
120  117 
{ 
121  118 
cout << "\n"; 
...  ...  
146  143 
return; 
147  144 
} 
148  145  
149  
150 
/*void gauss_hermite(BOR_CONST dvector& _x,BOR_CONST dvector& _w) 

151 
{ 

152 
cout << "gauss_hermite" << endl; 

153 
dvector x=(dvector&) _x; 

154 
dvector w=(dvector&) _w; 

155 
int ximin=x.indexmin(); 

156 
int ximax=x.indexmax(); 

157 
int wimin=w.indexmin(); 

158 
int wimax=w.indexmax(); 

159 
if (ximin != wimin  ximax != wimax) 

160 
{ 

161 
cerr << " Vector size mismatch in Gauss_hermite routine" << endl; 

162 
ad_exit(1); 

163 
} 

164 


165 
x.shift(1); 

166 
w.shift(1); 

167 
int n=x.size(); 

168 
int i,its,j,m; 

169 
double p1,p2,p3,pp,z,z1; 

170  
171 
m=(n+1)/2; 

172 
for (i=1;i<=m;i++) 

173 
{ 

174 
if (i == 1) 

175 
{ 

176 
z=sqrt((double)(2*n+1))1.85575*pow((double)(2*n+1),0.16667); 

177 
} 

178 
else if (i == 2) 

179 
{ 

180 
z = 1.14*pow((double)n,0.426)/z; 
Also available in: Unified diff