ADMB Documentation  11.1.2501
 All Classes Files Functions Variables Typedefs Friends Defines
vspline.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: vspline.cpp 1708 2014-02-28 21:43:41Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00013 #include <fvar.hpp>
00014 
00015 dvar_vector spline(const dvector &x,const dvar_vector&y,double yp1,double ypn);
00016 dvar_vector spline(const dvector &x,const dvar_vector&y,dvariable yp1,
00017   dvariable ypn);
00018 dvariable spline_cubic_val2(int n, const dvector& t, const prevariable tval,
00019   const dvar_vector& y, const dvar_vector& ypp);
00020 dvariable spline_cubic_val(int n, const dvector& t, double tval,
00021   const dvar_vector& y, const dvar_vector& ypp);
00022 
00023 dvar_vector spline_cubic_set (int n, const dvector& t, const dvar_vector& y,
00024   int ibcbeg, dvariable ybcbeg, int ibcend, dvariable ybcend );
00025 
00035 dvariable splint(const dvector& _xa,const dvar_vector& _ya,
00036     const dvar_vector& _y2a,double x)
00037 {
00038   RETURN_ARRAYS_INCREMENT();
00039   dvariable ret = spline_cubic_val(_xa.size(), _xa, x, _ya, _y2a);
00040   RETURN_ARRAYS_DECREMENT();
00041   return ret;
00042 }
00043 
00053 dvariable splint(const dvector& _xa,const dvar_vector& _ya,
00054     const dvar_vector& _y2a, const prevariable& _x)
00055 {
00056   RETURN_ARRAYS_INCREMENT();
00057   dvariable ret = spline_cubic_val2(_xa.size(), _xa, _x, _ya, _y2a);
00058   RETURN_ARRAYS_DECREMENT();
00059   return ret;
00060 }
00061 
00066 vcubic_spline_function::vcubic_spline_function(const dvector & _x,
00067   const dvar_vector& _y,dvariable yp1,dvariable ypn) : x(_x) , y(_y)
00068 {
00069   y2.allocate(x);
00070   y2=spline(x,y,yp1,ypn);
00071 }
00072 
00077 vcubic_spline_function::vcubic_spline_function(const dvector & _x,
00078   const dvar_vector& _y,double yp1,double ypn) : x(_x) , y(_y)
00079 {
00080   y2.allocate(x);
00081   y2=spline(x,y,yp1,ypn);
00082 }
00083 
00088 dvariable vcubic_spline_function::operator () (double u)
00089 {
00090   // need to deal with u<x(1) or u>x(2)
00091   return splint(x,y,y2,u);
00092 }
00093 
00098 dvar_vector vcubic_spline_function::operator () (const dvector& u)
00099 {
00100   int mmin=u.indexmin();
00101   int mmax=u.indexmax();
00102   dvar_vector z(mmin,mmax);
00103   for (int i=mmin;i<=mmax;i++)
00104   {
00105     z(i)=splint(x,y,y2,u(i));
00106   }
00107   return z;
00108 }
00109 
00114 dvar_vector vcubic_spline_function::operator () (const dvar_vector& u)
00115 {
00116   int mmin=u.indexmin();
00117   int mmax=u.indexmax();
00118   dvar_vector z(mmin,mmax);
00119   for (int i=mmin;i<=mmax;i++)
00120   {
00121     z(i)=splint(x,y,y2,u(i));
00122   }
00123   return z;
00124 }
00125 
00137 dvar_vector spline(const dvector &_x,const dvar_vector&_y,dvariable yp1,
00138     dvariable ypn)
00139 {
00140   int ibcbeg, ibcend;
00141   dvariable ybcbeg, ybcend;
00142   dvector x = _x;
00143   x.shift(0);
00144   dvar_vector y = _y;
00145   y.shift(0);
00146   if(value(yp1) > 0.99e30 )
00147   {
00148     ibcbeg = 2;
00149     ybcbeg = 0.0;
00150   }
00151   else
00152   {
00153     ibcbeg = 1;
00154     ybcbeg = yp1;
00155   }
00156   if(value(ypn) > 0.99e30 )
00157   {
00158     ibcend = 2;
00159     ybcend = 0.0;
00160   }
00161   else
00162   {
00163     ibcend = 1;
00164     ybcend = ypn;
00165   }
00166 
00167   dvar_vector ret = spline_cubic_set(x.size(), x, y, ibcbeg, ybcbeg, ibcend,
00168     ybcend);
00169   ret.shift(_x.indexmin());
00170   return ret;
00171 }
00172 
00184 dvar_vector spline(const dvector &_x,const dvar_vector&_y,double yp1,
00185     double ypn)
00186 {
00187   int ibcbeg, ibcend;
00188   dvariable ybcbeg, ybcend;
00189   dvector x = _x;
00190   x.shift(0);
00191   dvar_vector y = _y;
00192   y.shift(0);
00193   if(yp1 > 0.99e30 )
00194   {
00195     ibcbeg = 2;
00196     ybcbeg = 0.0;
00197   }
00198   else
00199   {
00200     ibcbeg = 1;
00201     ybcbeg = yp1;
00202   }
00203   if(ypn > 0.99e30 )
00204   {
00205     ibcend = 2;
00206     ybcend = 0.0;
00207   }
00208   else
00209   {
00210     ibcend = 1;
00211     ybcend = ypn;
00212   }
00213 
00214   dvar_vector ret = spline_cubic_set(x.size(), x, y, ibcbeg, ybcbeg, ibcend,
00215     ybcend);
00216   ret.shift(_x.indexmin());
00217   return ret;
00218 }
00219 
00231 dvar_vector spline(const dvector &_x,const dvar_vector&_y,dvariable yp1,
00232     double ypn)
00233 {
00234   int ibcbeg, ibcend;
00235   dvariable ybcbeg, ybcend;
00236   dvector x = _x;
00237   x.shift(0);
00238   dvar_vector y = _y;
00239   y.shift(0);
00240   if(value(yp1) > 0.99e30 )
00241   {
00242     ibcbeg = 2;
00243     ybcbeg = 0.0;
00244   }
00245   else
00246   {
00247     ibcbeg = 1;
00248     ybcbeg = yp1;
00249   }
00250   if(ypn > 0.99e30 )
00251   {
00252     ibcend = 2;
00253     ybcend = 0.0;
00254   }
00255   else
00256   {
00257     ibcend = 1;
00258     ybcend = ypn;
00259   }
00260 
00261   dvar_vector ret = spline_cubic_set(x.size(), x, y, ibcbeg, ybcbeg, ibcend,
00262     ybcend);
00263   ret.shift(_x.indexmin());
00264   return ret;
00265 }
00266 
00279 dvariable spline_cubic_val(int n,  const dvector& _t, double tval,
00280   const dvar_vector& _y, const dvar_vector& _ypp)
00281 //
00282 //  Purpose:
00283 //
00284 //    SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
00285 //
00286 //  Discussion:
00287 //
00288 //    SPLINE_CUBIC_SET must have already been called to define the values of
00289 //    YPP.
00290 //
00291 //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00292 //    the spline is
00293 //
00294 //      SPL(T) = A
00295 //             + B * ( T - T(IVAL) )
00296 //             + C * ( T - T(IVAL) )**2
00297 //             + D * ( T - T(IVAL) )**3
00298 //
00299 //    Here:
00300 //      A = Y(IVAL)
00301 //      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00302 //        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00303 //      C = YPP(IVAL) / 2
00304 //      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00305 //
00306 //  Licensing:
00307 //
00308 //    This code is distributed under the GNU LGPL license.
00309 //
00310 //  Modified:
00311 //
00312 //    04 February 1999
00313 //
00314 //  Author:
00315 //
00316 //    John Burkardt
00317 //
00318 //  Modified By:
00319 //
00320 //    Derek Seiple (20 August 2010)
00321 //
00322 //  Parameters:
00323 //
00324 //    Input, int N, the number of knots.
00325 //
00326 //    Input, double Y[N], the data values at the knots.
00327 //
00328 //    Input, double T[N], the knot values.
00329 //
00330 //    Input, double TVAL, a point, typically between T[0] and T[n-1], at
00331 //    which the spline is to be evalulated.  If TVAL lies outside
00332 //    this range, extrapolation is used.
00333 //
00334 //    Input, double YPP[N], the second derivatives of the spline at
00335 //    the knots.
00336 //
00337 //    Output, double SPLINE_VAL, the value of the spline at TVAL.
00338 //
00339 {
00340   dvariable dt = 0.0;
00341   dvariable h = 0.0;
00342   int i = 0;
00343   int ival = 0;
00344   dvariable yval = 0.0;
00345   dvector& t = (dvector&)_t;
00346   int lbt = t.indexmin();
00347   t.shift(0);
00348   dvar_vector& y = (dvar_vector&)_y;
00349   int lby = y.indexmin();
00350   y.shift(0);
00351   dvar_vector& ypp = (dvar_vector&)_ypp;
00352   int lbypp = ypp.indexmin();
00353   ypp.shift(0);
00354 
00355   //  Determine the interval [ T(I), T(I+1) ] that contains TVAL.
00356   //  Values below T[0] or above T[N-1] use extrapolation.
00357   ival = n - 2;
00358 
00359   for ( i = 0; i < n-1; i++ )
00360   {
00361     if ( tval < t[i+1] )
00362     {
00363       ival = i;
00364       break;
00365     }
00366   }
00367 
00368   //  In the interval I, the polynomial is in terms of a normalized
00369   //  coordinate between 0 and 1.
00370   dt = tval - t[ival];
00371   h = t[ival+1] - t[ival];
00372 
00373   yval = y[ival]
00374     + dt * ( ( y[ival+1] - y[ival] ) / h
00375            - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
00376     + dt * ( 0.5 * ypp[ival]
00377     + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );
00378 
00379   /* *ypval = ( y[ival+1] - y[ival] ) / h
00380     - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
00381     + dt * ( ypp[ival]
00382     + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );
00383 
00384   *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;*/
00385   t.shift(lbt);
00386   y.shift(lby);
00387   ypp.shift(lbypp);
00388 
00389   return yval;
00390 }
00391 
00404 dvariable spline_cubic_val2(int n, const dvector& _t, const prevariable tval,
00405   const dvar_vector& _y, const dvar_vector& _ypp)
00406 //
00407 //  Purpose:
00408 //
00409 //    SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
00410 //
00411 //  Discussion:
00412 //
00413 //    SPLINE_CUBIC_SET must have already been called to define the values of
00414 //    YPP.
00415 //
00416 //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00417 //    the spline is
00418 //
00419 //      SPL(T) = A
00420 //             + B * ( T - T(IVAL) )
00421 //             + C * ( T - T(IVAL) )**2
00422 //             + D * ( T - T(IVAL) )**3
00423 //
00424 //    Here:
00425 //      A = Y(IVAL)
00426 //      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00427 //        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00428 //      C = YPP(IVAL) / 2
00429 //      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00430 //
00431 //  Licensing:
00432 //
00433 //    This code is distributed under the GNU LGPL license.
00434 //
00435 //  Modified:
00436 //
00437 //    04 February 1999
00438 //
00439 //  Author:
00440 //
00441 //    John Burkardt
00442 //
00443 //  Modified By:
00444 //
00445 //    Derek Seiple (20 August 2010)
00446 //
00447 //  Parameters:
00448 //
00449 //    Input, int N, the number of knots.
00450 //
00451 //    Input, double Y[N], the data values at the knots.
00452 //
00453 //    Input, double T[N], the knot values.
00454 //
00455 //    Input, double TVAL, a point, typically between T[0] and T[n-1], at
00456 //    which the spline is to be evalulated.  If TVAL lies outside
00457 //    this range, extrapolation is used.
00458 //
00459 //    Input, double YPP[N], the second derivatives of the spline at
00460 //    the knots.
00461 //
00462 //    Output, double SPLINE_VAL, the value of the spline at TVAL.
00463 //
00464 {
00465   dvariable dt = 0.0;
00466   dvariable h = 0.0;
00467   int i = 0;
00468   int ival = 0;
00469   dvariable yval = 0.0;
00470   dvector& t = (dvector&)_t;
00471   int lbt = t.indexmin();
00472   t.shift(0);
00473   dvar_vector& y = (dvar_vector&)_y;
00474   int lby = y.indexmin();
00475   y.shift(0);
00476   dvar_vector& ypp = (dvar_vector&)_ypp;
00477   int lbypp = ypp.indexmin();
00478   ypp.shift(0);
00479 
00480   //  Determine the interval [ T(I), T(I+1) ] that contains TVAL.
00481   //  Values below T[0] or above T[N-1] use extrapolation.
00482   ival = n - 2;
00483 
00484   for ( i = 0; i < n-1; i++ )
00485   {
00486     if ( tval < t[i+1] )
00487     {
00488       ival = i;
00489       break;
00490     }
00491   }
00492 
00493   //  In the interval I, the polynomial is in terms of a normalized
00494   //  coordinate between 0 and 1.
00495   dt = tval - t[ival];
00496   h = t[ival+1] - t[ival];
00497 
00498   yval = y[ival]
00499     + dt * ( ( y[ival+1] - y[ival] ) / h
00500            - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
00501     + dt * ( 0.5 * ypp[ival]
00502     + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );
00503 
00504   /* *ypval = ( y[ival+1] - y[ival] ) / h
00505     - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
00506     + dt * ( ypp[ival]
00507     + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );
00508 
00509   *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;*/
00510   t.shift(lbt);
00511   y.shift(lby);
00512   ypp.shift(lbypp);
00513 
00514   return yval;
00515 }
00516 
00525 dvar_vector d3_np_fs ( int n, const dvar_vector& _a, const dvar_vector& _b)
00526 //
00527 //  Purpose:
00528 //
00529 //    D3_NP_FS factors and solves a D3 system.
00530 //
00531 //  Discussion:
00532 //
00533 //    The D3 storage format is used for a tridiagonal matrix.
00534 //    The superdiagonal is stored in entries (1,2:N), the diagonal in
00535 //    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
00536 //    original matrix is "collapsed" vertically into the array.
00537 //
00538 //    This algorithm requires that each diagonal entry be nonzero.
00539 //    It does not use pivoting, and so can fail on systems that
00540 //    are actually nonsingular.
00541 //
00542 //  Example:
00543 //
00544 //    Here is how a D3 matrix of order 5 would be stored:
00545 //
00546 //       *  A12 A23 A34 A45
00547 //      A11 A22 A33 A44 A55
00548 //      A21 A32 A43 A54  *
00549 //
00550 //  Licensing:
00551 //
00552 //    This code is distributed under the GNU LGPL license.
00553 //
00554 //  Modified:
00555 //
00556 //    15 November 2003
00557 //
00558 //  Author:
00559 //
00560 //    John Burkardt
00561 //
00562 //  Modified By:
00563 //
00564 //    Derek Seiple (20 August 2010)
00565 //
00566 //  Parameters:
00567 //
00568 //    Input, int N, the order of the linear system.
00569 //
00570 //    Input/output, double A[3*N].
00571 //    On input, the nonzero diagonals of the linear system.
00572 //    On output, the data in these vectors has been overwritten
00573 //    by factorization information.
00574 //
00575 //    Input, double B[N], the right hand side.
00576 //
00577 //    Output, double D3_NP_FS[N], the solution of the linear system.
00578 //    This is NULL if there was an error because one of the diagonal
00579 //    entries was zero.
00580 //
00581 {
00582   ADUNCONST(dvar_vector,a)
00583   ADUNCONST(dvar_vector,b)
00584   int i;
00585   dvariable xmult;
00586 
00587   //  Check.
00588   for ( i = 0; i < n; i++ )
00589   {
00590     if ( a[1+i*3] == 0.0 )
00591     {
00592       return NULL;
00593     }
00594   }
00595   dvar_vector x(0,n-1);
00596 
00597   for ( i = 0; i < n; i++ )
00598   {
00599     x[i] = b[i];
00600   }
00601 
00602   for ( i = 1; i < n; i++ )
00603   {
00604     xmult = a[2+(i-1)*3] / a[1+(i-1)*3];
00605     a[1+i*3] = a[1+i*3] - xmult * a[0+i*3];
00606     x[i] = x[i] - xmult * x[i-1];
00607   }
00608 
00609   x[n-1] = x[n-1] / a[1+(n-1)*3];
00610   for ( i = n-2; 0 <= i; i-- )
00611   {
00612     x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3];
00613   }
00614 
00615   return x;
00616 }
00617 
00638 dvar_vector spline_cubic_set (int n, const dvector& t, const dvar_vector& y,
00639   int ibcbeg, dvariable ybcbeg, int ibcend, dvariable ybcend )
00640 //
00641 //  Purpose:
00642 //
00643 //    SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic
00644 //    spline.
00645 //
00646 //  Discussion:
00647 //
00648 //    For data interpolation, the user must call SPLINE_SET to determine
00649 //    the second derivative data, passing in the data to be interpolated,
00650 //    and the desired boundary conditions.
00651 //
00652 //    The data to be interpolated, plus the SPLINE_SET output, defines
00653 //    the spline.  The user may then call SPLINE_VAL to evaluate the
00654 //    spline at any point.
00655 //
00656 //    The cubic spline is a piecewise cubic polynomial.  The intervals
00657 //    are determined by the "knots" or abscissas of the data to be
00658 //    interpolated.  The cubic spline has continous first and second
00659 //    derivatives over the entire interval of interpolation.
00660 //
00661 //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00662 //    the spline is
00663 //
00664 //      SPL(T) = A(IVAL)
00665 //             + B(IVAL) * ( T - T(IVAL) )
00666 //             + C(IVAL) * ( T - T(IVAL) )**2
00667 //             + D(IVAL) * ( T - T(IVAL) )**3
00668 //
00669 //    If we assume that we know the values Y(*) and YPP(*), which represent
00670 //    the values and second derivatives of the spline at each knot, then
00671 //    the coefficients can be computed as:
00672 //
00673 //      A(IVAL) = Y(IVAL)
00674 //      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00675 //        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00676 //      C(IVAL) = YPP(IVAL) / 2
00677 //      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00678 //
00679 //    Since the first derivative of the spline is
00680 //
00681 //      SPL'(T) =     B(IVAL)
00682 //              + 2 * C(IVAL) * ( T - T(IVAL) )
00683 //              + 3 * D(IVAL) * ( T - T(IVAL) )**2,
00684 //
00685 //    the requirement that the first derivative be continuous at interior
00686 //    knot I results in a total of N-2 equations, of the form:
00687 //
00688 //      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
00689 //      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
00690 //
00691 //    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
00692 //
00693 //      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
00694 //      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
00695 //      + YPP(IVAL-1) * H(IVAL-1)
00696 //      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
00697 //      =
00698 //      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
00699 //      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
00700 //
00701 //    or
00702 //
00703 //      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
00704 //      + YPP(IVAL) * H(IVAL)
00705 //      =
00706 //      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
00707 //      - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
00708 //
00709 //    Boundary conditions must be applied at the first and last knots.
00710 //    The resulting tridiagonal system can be solved for the YPP values.
00711 //
00712 //  Licensing:
00713 //
00714 //    This code is distributed under the GNU LGPL license.
00715 //
00716 //  Modified:
00717 //
00718 //    06 February 2004
00719 //
00720 //  Author:
00721 //
00722 //    John Burkardt
00723 //
00724 //  Modified By:
00725 //
00726 //    Derek Seiple (20 August 2010)
00727 //
00728 //  Parameters:
00729 //
00730 //    Input, int N, the number of data points.  N must be at least 2.
00731 //    In the special case where N = 2 and IBCBEG = IBCEND = 0, the
00732 //    spline will actually be linear.
00733 //
00734 //    Input, double T[N], the knot values, that is, the points were data is
00735 //    specified.  The knot values should be distinct, and increasing.
00736 //
00737 //    Input, double Y[N], the data values to be interpolated.
00738 //
00739 //    Input, int IBCBEG, left boundary condition flag:
00740 //      0: the cubic spline should be a quadratic over the first interval;
00741 //      1: the first derivative at the left endpoint should be YBCBEG;
00742 //      2: the second derivative at the left endpoint should be YBCBEG.
00743 //
00744 //    Input, double YBCBEG, the values to be used in the boundary
00745 //    conditions if IBCBEG is equal to 1 or 2.
00746 //
00747 //    Input, int IBCEND, right boundary condition flag:
00748 //      0: the cubic spline should be a quadratic over the last interval;
00749 //      1: the first derivative at the right endpoint should be YBCEND;
00750 //      2: the second derivative at the right endpoint should be YBCEND.
00751 //
00752 //    Input, double YBCEND, the values to be used in the boundary
00753 //    conditions if IBCEND is equal to 1 or 2.
00754 //
00755 //    Output, double SPLINE_CUBIC_SET[N], the second derivatives of the cubic
00756 //    spline.
00757 //
00758 {
00759   dvar_vector a(0,3*n-1);
00760   dvar_vector b(0,n-1);
00761   int i;
00762 
00763   //  Check.
00764   if ( n <= 1 )
00765   {
00766     cout << "\n";
00767     cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00768     cout << "  The number of data points N must be at least 2.\n";
00769     cout << "  The input value is " << n << ".\n";
00770     //return NULL;
00771   }
00772 
00773   for ( i = 0; i < n - 1; i++ )
00774   {
00775     if ( t[i+1] <= t[i] )
00776     {
00777       cout << "\n";
00778       cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00779       cout << "  The knots must be strictly increasing, but\n";
00780       cout << "  T(" << i   << ") = " << t[i]   << "\n";
00781       cout << "  T(" << i+1 << ") = " << t[i+1] << "\n";
00782     }
00783   }
00784 
00785   //  Set up the first equation.
00786   if ( ibcbeg == 0 )
00787   {
00788     b[0] = 0.0;
00789     a[1+0*3] = 1.0;
00790     a[0+1*3] = -1.0;
00791   }
00792   else if ( ibcbeg == 1 )
00793   {
00794     b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg;
00795     a[1+0*3] = ( t[1] - t[0] ) / 3.0;
00796     a[0+1*3] = ( t[1] - t[0] ) / 6.0;
00797   }
00798   else if ( ibcbeg == 2 )
00799   {
00800     b[0] = ybcbeg;
00801     a[1+0*3] = 1.0;
00802     a[0+1*3] = 0.0;
00803   }
00804   else
00805   {
00806     cout << "\n";
00807     cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00808     cout << "  IBCBEG must be 0, 1 or 2.\n";
00809     cout << "  The input value is " << ibcbeg << ".\n";
00810   }
00811 
00812   //  Set up the intermediate equations.
00813   for ( i = 1; i < n-1; i++ )
00814   {
00815     b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] )
00816       - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] );
00817     a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0;
00818     a[1+ i   *3] = ( t[i+1] - t[i-1] ) / 3.0;
00819     a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0;
00820   }
00821 
00822   //  Set up the last equation.
00823   if ( ibcend == 0 )
00824   {
00825     b[n-1] = 0.0;
00826     a[2+(n-2)*3] = -1.0;
00827     a[1+(n-1)*3] = 1.0;
00828   }
00829   else if ( ibcend == 1 )
00830   {
00831     b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] );
00832     a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0;
00833     a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0;
00834   }
00835   else if ( ibcend == 2 )
00836   {
00837     b[n-1] = ybcend;
00838     a[2+(n-2)*3] = 0.0;
00839     a[1+(n-1)*3] = 1.0;
00840   }
00841   else
00842   {
00843     cout << "\n";
00844     cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00845     cout << "  IBCEND must be 0, 1 or 2.\n";
00846     cout << "  The input value is " << ibcend << ".\n";
00847   }
00848 
00849   //  Solve the linear system.
00850   if ( n == 2 && ibcbeg == 0 && ibcend == 0 )
00851   {
00852     dvar_vector ret(0,1);
00853     ret(0) = 0.0;
00854     ret(1) = 0.0;
00855     return ret;
00856   }
00857   else
00858   {
00859     dvar_vector ypp = d3_np_fs ( n, a, b );
00860     dvar_vector ret(0,n-1);
00861     if ( !ypp )
00862     {
00863       cout << "\n";
00864       cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00865       cout << "  The linear system could not be solved.\n";
00866     }
00867     ret = ypp;
00868     return ret;
00869   }
00870 }