ADMB Documentation  11.1.2402
 All Classes Files Functions Variables Typedefs Friends Defines
gaussher.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: gaussher.cpp 2259 2014-08-27 23:16:51Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00014 #include <fvar.hpp>
00015 //static double eps=3.0e-14;
00016 //static double pim=0.7511255444649427;
00017 //static int maxit=50;
00018 
00019 void imtqlx ( const dvector& _d, const dvector& _e, const dvector& _z );
00020 
00025 double sign(const double x)
00026 {
00027   return x < 0.0 ? -1.0 : 1.0;
00028 }
00029 
00036 void gauss_hermite (const dvector& _t,const dvector& _wts)
00037 //
00038 //  Purpose:
00039 //
00040 //    computes a Gauss quadrature formula with simple knots.
00041 //
00042 //  Discussion:
00043 //
00044 //    This routine computes all the knots and weights of a Gauss quadrature
00045 //    formula with a classical weight function with default values for A and B,
00046 //    and only simple knots.
00047 //
00048 //    There are no moments checks and no printing is done.
00049 //
00050 //    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
00051 //
00052 //  Licensing:
00053 //
00054 //    This code is distributed under the GNU LGPL license.
00055 //
00056 //  Modified:
00057 //
00058 //    September 2010 by Derek Seiple
00059 //
00060 //  Author:
00061 //
00062 //    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
00063 //    C++ version by John Burkardt.
00064 //
00065 //  Reference:
00066 //
00067 //    Sylvan Elhay, Jaroslav Kautsky,
00068 //    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
00069 //    Interpolatory Quadrature,
00070 //    ACM Transactions on Mathematical Software,
00071 //    Volume 13, Number 4, December 1987, pages 399-415.
00072 //
00073 //  Parameters:
00074 //
00075 //    Output, double T[NT], the knots.
00076 //
00077 //    Output, double WTS[NT], the weights.
00078 //
00079 {
00080   dvector t=(dvector&) _t;
00081   dvector wts=(dvector&) _wts;
00082 
00083   if( t.indexmax()!=wts.indexmax() )
00084   {
00085     cerr << "Incompatible sizes in void "
00086         << "void gauss_hermite (const dvector& _t,const dvector& _wts)" << endl;
00087     ad_exit(-1);
00088   }
00089 
00090   int lb = t.indexmin();
00091   int ub = t.indexmax();
00092 
00093   dvector aj(lb,ub);
00094   dvector bj(lb,ub);
00095   double zemu;
00096   int i;
00097 
00098   //  Get the Jacobi matrix and zero-th moment.
00099   zemu = 1.772453850905516;
00100   for ( i = lb; i <= ub; i++ )
00101   {
00102     aj(i) = 0.0;
00103   }
00104   for ( i = lb; i <= ub; i++ )
00105   {
00106     bj(i) = sqrt((i-lb+1)/2.0);
00107   }
00108 
00109   //  Compute the knots and weights.
00110   if ( zemu <= 0.0 ) //  Exit if the zero-th moment is not positive.
00111   {
00112     cout << "\n";
00113     cout << "SGQF - Fatal error!\n";
00114     cout << "  ZEMU <= 0.\n";
00115     exit ( 1 );
00116   }
00117 
00118   //  Set up vectors for IMTQLX.
00119   for ( i = lb; i <= ub; i++ )
00120   {
00121     t(i) = aj(i);
00122   }
00123   wts(lb) = sqrt ( zemu );
00124   for ( i = lb+1; i <= ub; i++ )
00125   {
00126     wts(i) = 0.0;
00127   }
00128 
00129   //  Diagonalize the Jacobi matrix.
00130   imtqlx ( t, bj, wts );
00131 
00132   for ( i = lb; i <= ub; i++ )
00133   {
00134     wts(i) = wts(i) * wts(i);
00135   }
00136 
00137   return;
00138 }
00139 
00148 void gauss_legendre( double a, double b, const dvector& _t,
00149   const dvector& _wts )
00150 //
00151 //  Purpose:
00152 //
00153 //    computes knots and weights of a Gauss-Legendre quadrature formula.
00154 //
00155 //  Discussion:
00156 //
00157 //    The user may specify the interval (A,B).
00158 //
00159 //    Only simple knots are produced.
00160 //
00161 //    Use routine EIQFS to evaluate this quadrature formula.
00162 //
00163 //  Licensing:
00164 //
00165 //    This code is distributed under the GNU LGPL license.
00166 //
00167 //  Modified:
00168 //
00169 //    September 2010 by Derek Seiple
00170 //
00171 //  Author:
00172 //
00173 //    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
00174 //    C++ version by John Burkardt.
00175 //
00176 //  Reference:
00177 //
00178 //    Sylvan Elhay, Jaroslav Kautsky,
00179 //    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
00180 //    Interpolatory Quadrature,
00181 //    ACM Transactions on Mathematical Software,
00182 //    Volume 13, Number 4, December 1987, pages 399-415.
00183 //
00184 //  Parameters:
00185 //
00186 //    Input, double A, B, the interval endpoints, or
00187 //    other parameters.
00188 //
00189 //    Output, double T[NT], the knots.
00190 //
00191 //    Output, double WTS[NT], the weights.
00192 //
00193 {
00194   dvector t=(dvector&) _t;
00195   dvector wts=(dvector&) _wts;
00196 
00197   if( t.indexmax()!=wts.indexmax() )
00198   {
00199     cerr << "Incompatible sizes in void "
00200 "mygauss_legendre(double a, double b, const dvector& _t, const dvector& _wts)"
00201     << endl;
00202     ad_exit(-1);
00203   }
00204 
00205   t.shift(0);
00206   wts.shift(0);
00207   int nt = t.indexmax() + 1;
00208   int ub = nt-1;
00209 
00210   int i;
00211   int k;
00212   int l;
00213   double al;
00214   double ab;
00215   double abi;
00216   double abj;
00217   double be;
00218   double p;
00219   double shft;
00220   double slp;
00221   double temp;
00222   double tmp;
00223   double zemu;
00224 
00225   //  Compute the Gauss quadrature formula for default values of A and B.
00226   dvector aj(0,ub);
00227   dvector bj(0,ub);
00228 
00229   ab = 0.0;
00230   zemu = 2.0 / ( ab + 1.0 );
00231   for ( i = 0; i < nt; i++ )
00232   {
00233     aj[i] = 0.0;
00234   }
00235   for ( i = 1; i <= nt; i++ )
00236   {
00237     abi = i + ab * ( i % 2 );
00238     abj = 2 * i + ab;
00239     bj[i-1] = sqrt ( abi * abi / ( abj * abj - 1.0 ) );
00240   }
00241 
00242   //  Compute the knots and weights.
00243   if ( zemu <= 0.0 )  //  Exit if the zero-th moment is not positive.
00244   {
00245     cout << "\n";
00246     cout << "Fatal error!\n";
00247     cout << "  ZEMU <= 0.\n";
00248     exit ( 1 );
00249   }
00250 
00251   //  Set up vectors for IMTQLX.
00252   for ( i = 0; i < nt; i++ )
00253   {
00254     t[i] = aj[i];
00255   }
00256   wts[0] = sqrt ( zemu );
00257   for ( i = 1; i < nt; i++ )
00258   {
00259     wts[i] = 0.0;
00260   }
00261 
00262   //  Diagonalize the Jacobi matrix.
00263   imtqlx (t, bj, wts );
00264 
00265   for ( i = 0; i < nt; i++ )
00266   {
00267     wts[i] = wts[i] * wts[i];
00268   }
00269 
00270   //  Prepare to scale the quadrature formula to other weight function with
00271   //  valid A and B.
00272   ivector mlt(0,ub);
00273   for ( i = 0; i < nt; i++ )
00274   {
00275     mlt[i] = 1;
00276   }
00277   ivector ndx(0,ub);
00278   for ( i = 0; i < nt; i++ )
00279   {
00280     ndx[i] = i + 1;
00281   }
00282 
00283   dvector st(0,ub);
00284   dvector swts(0,ub);
00285 
00286   temp = 3.0e-14;
00287   al = 0.0;
00288   be = 0.0;
00289   if ( fabs ( b - a ) <= temp )
00290   {
00291     cout << "\n";
00292     cout << "Fatal error!\n";
00293     cout << "  |B - A| too small.\n";
00294     exit ( 1 );
00295   }
00296   shft = ( a + b ) / 2.0;
00297   slp = ( b - a ) / 2.0;
00298 
00299   p = pow ( slp, al + be + 1.0 );
00300 
00301   for ( k = 0; k < nt; k++ )
00302   {
00303     st[k] = shft + slp * t[k];
00304     l = abs ( ndx[k] );
00305 
00306     if ( l != 0 )
00307     {
00308       tmp = p;
00309       for ( i = l - 1; i <= l - 1 + mlt[k] - 1; i++ )
00310       {
00311         swts[i] = wts[i] * tmp;
00312         tmp = tmp * slp;
00313       }
00314     }
00315   }
00316 
00317   for(i=0;i<nt;i++)
00318   {
00319     t(i) = st(ub-i);
00320     wts(i) = swts(ub-i);
00321   }
00322 
00323   return;
00324 }
00325 
00331 void gauss_legendre(const dvector& _x, const dvector& _w)
00332 {
00333   gauss_legendre(0,1,_x,_w);
00334 }
00335 
00343 void normalized_gauss_hermite(const dvector& _x, const dvector& _w)
00344 {
00345   dvector& x=(dvector&) _x;
00346   dvector& w=(dvector&) _w;
00347   gauss_hermite(x,w);
00348   w=elem_prod(w,exp(square(x)));
00349   x*=sqrt(2.0);
00350   w*=sqrt(2.0);
00351 }
00352 
00359 void imtqlx( const dvector& _d, const dvector& _e, const dvector& _z )
00360 //
00361 //  Purpose:
00362 //
00363 //    IMTQLX diagonalizes a symmetric tridiagonal matrix.
00364 //
00365 //  Discussion:
00366 //
00367 //    This routine is a slightly modified version of the EISPACK routine to
00368 //    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
00369 //
00370 //    The authors thank the authors of EISPACK for permission to use this
00371 //    routine.
00372 //
00373 //    It has been modified to produce the product Q' * Z, where Z is an input
00374 //    vector and Q is the orthogonal matrix diagonalizing the input matrix.
00375 //   The changes consist (essentialy) of applying the orthogonal transformations
00376 //    directly to Z as they are generated.
00377 //
00378 //  Licensing:
00379 //
00380 //    This code is distributed under the GNU LGPL license.
00381 //
00382 //  Modified:
00383 //
00384 //    September 2010 by Derek Seiple
00385 //
00386 //  Author:
00387 //
00388 //    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
00389 //    C++ version by John Burkardt.
00390 //
00391 //  Reference:
00392 //
00393 //    Sylvan Elhay, Jaroslav Kautsky,
00394 //    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
00395 //    Interpolatory Quadrature,
00396 //    ACM Transactions on Mathematical Software,
00397 //    Volume 13, Number 4, December 1987, pages 399-415.
00398 //
00399 //    Roger Martin, James Wilkinson,
00400 //    The Implicit QL Algorithm,
00401 //    Numerische Mathematik,
00402 //    Volume 12, Number 5, December 1968, pages 377-383.
00403 //
00404 //  Parameters:
00405 //
00406 //    Input/output, double D(N), the diagonal entries of the matrix.
00407 //    On output, the information in D has been overwritten.
00408 //
00409 //    Input/output, double E(N), the subdiagonal entries of the
00410 //    matrix, in entries E(1) through E(N-1).  On output, the information in
00411 //    E has been overwritten.
00412 //
00413 //    Input/output, double Z(N).  On input, a vector.  On output,
00414 //    the value of Q' * Z, where Q is the matrix that diagonalizes the
00415 //    input symmetric tridiagonal matrix.
00416 //
00417 {
00418   dvector d=(dvector&) _d;
00419   dvector e=(dvector&) _e;
00420   dvector z=(dvector&) _z;
00421   if (d.indexmax()!=e.indexmax() || d.indexmax()!=z.indexmax() ||
00422     z.indexmax()!=e.indexmax() )
00423   {
00424     cerr << "Incompatible sizes in void "
00425          << "imtqlx ( const dvector& _d, const dvector& _e, const dvector& _z )"
00426          << endl;
00427     ad_exit(-1);
00428   }
00429 
00430   int lb = d.indexmin();
00431   d.shift(0);
00432   e.shift(0);
00433   z.shift(0);
00434   int n = d.indexmax() + 1;
00435 
00436   double b;
00437   double c;
00438   double f;
00439   double g;
00440   int i;
00441   int ii;
00442   int itn = 30;
00443   int j;
00444   int k;
00445   int l;
00446   int m;
00447   int mml;
00448   double p;
00449   double prec;
00450   double r;
00451   double s;
00452 
00453   prec = 3.0e-14;
00454 
00455   if ( n == 1 )
00456   {
00457     return;
00458   }
00459 
00460   e[n-1] = 0.0;
00461 
00462   for ( l = 1; l <= n; l++ )
00463   {
00464     j = 0;
00465     for ( ; ; )
00466     {
00467       for ( m = l; m <= n; m++ )
00468       {
00469         if ( m == n )
00470         {
00471           break;
00472         }
00473 
00474         if ( fabs ( e[m-1] ) <= prec * ( fabs ( d[m-1] ) + fabs ( d[m] ) ) )
00475         {
00476           break;
00477         }
00478       }
00479       p = d[l-1];
00480       if ( m == l )
00481       {
00482         break;
00483       }
00484       if ( itn <= j )
00485       {
00486         cout << "\n";
00487         cout << "IMTQLX - Fatal error!\n";
00488         cout << "  Iteration limit exceeded\n";
00489         exit ( 1 );
00490       }
00491       j = j + 1;
00492       g = ( d[l] - p ) / ( 2.0 * e[l-1] );
00493       r =  sqrt ( g * g + 1.0 );
00494       g = d[m-1] - p + e[l-1] / ( g + fabs ( r ) * sign ( g ) );
00495       s = 1.0;
00496       c = 1.0;
00497       p = 0.0;
00498       mml = m - l;
00499 
00500       for ( ii = 1; ii <= mml; ii++ )
00501       {
00502         i = m - ii;
00503         f = s * e[i-1];
00504         b = c * e[i-1];
00505 
00506         if ( fabs ( g ) <= fabs ( f ) )
00507         {
00508           c = g / f;
00509           r =  sqrt ( c * c + 1.0 );
00510           e[i] = f * r;
00511           s = 1.0 / r;
00512           c = c * s;
00513         }
00514         else
00515         {
00516           s = f / g;
00517           r =  sqrt ( s * s + 1.0 );
00518           e[i] = g * r;
00519           c = 1.0 / r;
00520           s = s * c;
00521         }
00522         g = d[i] - p;
00523         r = ( d[i-1] - g ) * s + 2.0 * c * b;
00524         p = s * r;
00525         d[i] = g + p;
00526         g = c * r - b;
00527         f = z[i];
00528         z[i] = s * z[i-1] + c * f;
00529         z[i-1] = c * z[i-1] - s * f;
00530       }
00531       d[l-1] = d[l-1] - p;
00532       e[l-1] = g;
00533       e[m-1] = 0.0;
00534     }
00535   }
00536 
00537   m = n;
00538 
00539   //  Sorting.
00540   for ( ii = 2; ii <= m; ii++ )
00541   {
00542     i = ii - 1;
00543     k = i;
00544     p = d[i-1];
00545 
00546     for ( j = ii; j <= n; j++ )
00547     {
00548       if ( d[j-1] > p )
00549       {
00550          k = j;
00551          p = d[j-1];
00552       }
00553     }
00554 
00555     if ( k != i )
00556     {
00557       d[k-1] = d[i-1];
00558       d[i-1] = p;
00559       p = z[i-1];
00560       z[i-1] = z[k-1];
00561       z[k-1] = p;
00562     }
00563   }
00564 
00565   d.shift(lb);
00566   e.shift(lb);
00567   z.shift(lb);
00568   return;
00569 }