Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

ali.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                        Numerical Math Package
00006  *                   Aitken-Lagrange interpolation
00007  *
00008  * This package allows one to interpolate a function value for a given
00009  * argument value using function values tabulated over either uniform or
00010  * non-uniform grid. The latter is specified by a vector of node point
00011  * abscissae. The uniform grid is specified by the size of a grid mesh
00012  * and the abscissa of the first grid point.
00013  *
00014  * Synopsis
00015  *      double ali(q,x0,s,y)
00016  *      double q                        Argument value specified
00017  *      double x0                       Abscissae for the 1. grid point
00018  *      double s                        Grid mesh, >0
00019  *      VECTOR y                        Vector of function values
00020  *                                      tabulated at points
00021  *                                      x0 + s*(i-y.q_lwb()))
00022  *                                      The vector must contain at 
00023  *                                      least 2 elements
00024  *
00025  *      double ali(q,x,y)
00026  *      const double q                  Argument value specified
00027  *      const VECTOR x                  Vector of grid node abscissae
00028  *      const VECTOR y                  Vector of function values
00029  *                                      tabulated at points x[i]
00030  *                                      The vector must contain at 
00031  *                                      least 2 elements
00032  * Output
00033  *      Both functions return the interpolated value of y(q)
00034  *      Interpolating process finishes either
00035  *              - if the difference between two successive interpolated
00036  *                values is absolutely less than EPSILON
00037  *              - if the absolute value of this difference stops
00038  *                diminishing
00039  *              - after (N-1) steps, N being the no. of elements in vector y
00040  *
00041  * Algorithm
00042  *      Aitken scheme of Lagrange interpolation
00043  *      Algorithm is described in the book
00044  *              "Mathematical software for computers", Institute of
00045  *               Mathematics of the Belorussian Academy of Sciences,
00046  *               Minsk, 1974, issue #4,
00047  *               p. 146 (description of ALI, DALI subroutines)
00048  *               p. 180 (description of ATSE, DATSE subroutines)
00049  *      The book essentially describes IBM's SSP package.
00050  *
00051  * $Id: ali.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00052  *
00053  ************************************************************************
00054  */
00055 
00056 
00057 #include "LAStreams.h"
00058 #include "math_num.h"
00059 #include "std.h"
00060 #include <float.h>
00061 
00062 
00063 namespace linalg 
00064 {
00065     using namespace linalg;
00066     
00067 //#define DEBUG
00068 
00069 /*
00070  *------------------------------------------------------------------------
00071  *                   Class that handles the interpolation
00072  */
00073 
00074 class ALInterp
00075 {
00076   Vector arg;                   // [1:n] Arranged table of arguments
00077   Vector val;                   // [1:n] Arranged table of function values
00078   double q;                     // Argument value the function is to be
00079                                 // interpolated at
00080 
00081   ALInterp(const ALInterp&);            // Deliberately unimplemented:
00082   void operator = (const ALInterp&);    // no copying/cloning allowed!
00083 
00084 public:
00085                                 // Construct the arranged tables for the
00086                                 // uniform grid
00087   ALInterp(const double q, const double x0, const double s, const Vector& y);
00088                                 // Construct the arranged tables for the
00089                                 // non-uniform grid
00090   ALInterp(const double q, const Vector& x, const Vector& y);
00091 
00092   double interpolate(void);     // Perform actual interpolation
00093 };
00094 
00095 /*
00096  *------------------------------------------------------------------------
00097  *
00098  *              Arranging data for the Aitken-Lagrange interpolation
00099  *
00100  * Abscissae (arg) and ordinates (val) of the grid points should be arranged
00101  * in such a way that the distance abs(q-arg[i]) would increase as i
00102  * increases. 
00103  * Here q is the point the function is to be interpolated at.
00104  *
00105  */
00106 
00107                                 // Construct the arranged tables for the
00108                                 // uniform grid
00109 ALInterp::ALInterp(const double _q, const double x0,
00110                    const double s, const Vector& y)
00111         : arg(y.q_no_elems()), val(y.q_no_elems()), q(_q)
00112 {
00113   const int n = y.q_no_elems();
00114   assure( n > 1, "Vector y (function values) must have at least 2 points");
00115   assure( s > 0, "The grid mesh has to be positive");
00116 
00117                         // First find the index of the grid node which
00118                         // is closest to q. Assign index 1 to this
00119                         // node. Then look at neighboring grid nodes
00120                         // and assign indices to them
00121                         // (kind of breadth-first search)
00122   int js = (int)( (q-x0)/s + 1.5 );     // Index j for the point x0+s*j
00123                                         // which is the closest to q
00124   if( js < 1 )                          // Check for the case of extrapolation
00125     js = 1;                             // to the left end
00126   else if( js > n )
00127     js = n;                             // or to the right end
00128 
00129                                         // Direction to the next closest
00130                                         // to q grid node
00131   bool right_pt_is_closer = q > x0 + (js-1)*s;
00132 
00133   register int jcurr = js, jleft = js, jright = js;
00134                                         // Pick up elements x0+s*i
00135                                         // in the neighborhood of q
00136   for(LAStreamOut args(arg), vals(val); !args.eof(); )
00137   {
00138     args.get() = x0 + (jcurr-1)*s;
00139     vals.get() = y(jcurr-1+y.q_lwb());  // Once the closest to q point js
00140     if( jright >= n )                   // is found, we pick up points
00141       right_pt_is_closer = false;       // alternatively to the right
00142     if( jleft <= 1 )                    // and to the left of the js
00143       right_pt_is_closer = true;        // further and further
00144     if( right_pt_is_closer )
00145       jcurr = ++jright, right_pt_is_closer = false;
00146     else
00147       jcurr = --jleft, right_pt_is_closer = true;
00148   }
00149 }
00150 
00151 
00152 static inline int fsign(const float f)          // Return the sign of f
00153 { return f < 0 ? -1 : f==0 ? 0 : 1; }
00154 
00155                                 // Construct the arranged tables for a
00156                                 // non-uniform grid
00157 ALInterp::ALInterp(const double _q, const Vector& x, const Vector& y)
00158         : arg(x.q_no_elems()), val(y.q_no_elems()), q(_q)
00159 {
00160   assure( y.q_no_elems() > 1,
00161           "Vector y (function values) must have at least 2 points");
00162   are_compatible(x,y);
00163 
00164                         // Selection is done by sorting x,y arrays
00165                         // in the way mentioned above. Fisrt an array
00166                         // of indices is created and sorted, then arg,
00167                         // val arrays are filled in using the sorted indices
00168   class index_permutation
00169   {
00170     struct El { int x_ind; float x_to_q; };     // x_to_q = |x[x_ind]-q|
00171     El * const permutation;
00172     const int n;
00173     static int comparison_func(const void * ip, const void * jp)
00174       { return fsign(((const El*)ip)->x_to_q - ((const El*)jp)->x_to_q); }
00175   public:
00176     index_permutation(const double q, const Vector& x) :
00177         permutation(new El[x.q_no_elems()]), n(x.q_no_elems())
00178     {
00179       register El * pp = permutation;
00180       for(register int i=x.q_lwb(); i<=x.q_upb(); i++,pp++)
00181         pp->x_ind = i, pp->x_to_q = fabs(q-x(i));
00182                                 // Sort indices so that
00183                                 // |q-x[x_ind[i]]| < |q-x[x_ind[j]]|
00184                                 // for all i<j
00185       qsort(permutation,n,sizeof(permutation[0]),comparison_func);
00186     }
00187     ~index_permutation(void) { delete permutation; }
00188                                 // Apply the permutation to x to get arg
00189                                 // and to y to get val
00190     void apply(Vector& arg, Vector& val, const Vector& x, const Vector& y)
00191     {
00192       register const El* pp = permutation;
00193       for(LAStreamOut args(arg), vals(val); !args.eof(); pp++)
00194         args.get() = x(pp->x_ind), vals.get() = y(pp->x_ind);
00195       assert(pp==permutation+n);
00196     }
00197   };
00198   
00199   index_permutation(q,x).apply(arg,val,x,y);
00200  }
00201 
00202 /*
00203  *------------------------------------------------------------------------
00204  *                      Aitken - Lagrange process
00205  *
00206  *  arg and val tables are assumed to be arranged in the proper way
00207  *
00208  */
00209 
00210 double ALInterp::interpolate()
00211 {
00212   LAStreamIn args(arg);
00213   LAStreamOut vals(val);
00214   register double valp = vals.peek();   // The best approximation found so far
00215   register double diffp = DBL_MAX;      // abs(valp - prev. to valp)
00216 
00217 #ifdef DEBUG
00218   arg.print("arg - interpolation nodes");
00219   val.print("Arranged table of function values");
00220 #endif
00221                         // Compute the j-th row of the Aitken scheme and
00222                         // place it in the 'val' array
00223   for(register int j=2; j<=val.q_upb(); j++)
00224   {
00225     register double argj = (args.get(), args.peek());
00226     register REAL&  valj = (vals.get(), vals.peek());
00227     args.rewind(); vals.rewind();
00228     for(register int i=1; i<=j-1; i++)
00229     {
00230       double argi = args.get();
00231       valj = ( vals.get()*(q-argj) - valj*(q-argi) ) / (argi - argj);
00232     }
00233 
00234 #ifdef DEBUG
00235     message("\nval(j) = %g, valp = %g, arg(j) = %g",valj,valp,argj);
00236 #endif
00237 
00238     register double diff = fabs( valj - valp );
00239 
00240     if( j>2 && diff == 0 )              // An exact result has been achieved
00241       break;
00242 
00243     if( j>4 && diff > diffp )           // Difference stoped diminishing
00244       break;                            // after the 4. step
00245 
00246     valp = valj;  diffp = diff;
00247   }
00248 
00249   return valp;
00250 }
00251 
00252 
00253 
00254 /* 
00255  *=======================================================================
00256  *                              Root modules
00257  */
00258 
00259                                 // Uniform mesh x[i] = x0 + s*(i-y.lwb)
00260 double ali(const double q, const double x0, const double s, const Vector& y)
00261 {
00262   ALInterp al(q,x0,s,y);
00263   return al.interpolate();
00264 }
00265 
00266                                 // Nonuniform grid with nodes in x[i]
00267 double ali(const double q, const Vector& x, const Vector& y)
00268 {
00269   ALInterp al(q,x,y);
00270   return al.interpolate();
00271 }
00272  
00273 }

Generated on Wed Dec 15 21:20:27 2004 for vuVolume by  doxygen 1.3.9.1