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
 * Gauss-Hermite Quadrature and
11
 * Gauss-Legendre Quadrature.
12
 */
8 13

  
14
#include <fvar.hpp>
9 15
static double eps=3.0e-14;
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
/** Gauss-Hermite 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
 * Gauss-Hermite quadature.
38
 * Computes a Gauss-Hermite 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((i-lb+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 zero-th 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;
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff