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