home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fs.zip / octave / octave-2.1.23 / liboctave / DASSL.cc < prev    next >
C/C++ Source or Header  |  2000-01-15  |  10KB  |  479 lines
  1. /*
  2.  
  3. Copyright (C) 1996, 1997 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #if defined (__GNUG__)
  24. #pragma implementation
  25. #endif
  26.  
  27. #ifdef HAVE_CONFIG_H
  28. #include <config.h>
  29. #endif
  30.  
  31. #include <cfloat>
  32. #include <cmath>
  33.  
  34. #include "DASSL.h"
  35. #include "f77-fcn.h"
  36. #include "lo-error.h"
  37.  
  38. extern "C"
  39. {
  40.   int F77_FCN (ddassl, DDASSL) (int (*)(const double&, double*,
  41.                     double*, double*, int&,
  42.                     double*, int*),
  43.                 const int&, double&, double*, double*,
  44.                 double&, const int*, const double&,
  45.                 const double&, int&, double*,
  46.                 const int&, int*, const int&,
  47.                 const double*, const int*, 
  48.                 int (*)(const double&, double*,
  49.                     double*, double*, const
  50.                     double&, double*, int*));
  51. }
  52.  
  53. static DAEFunc::DAERHSFunc user_fun;
  54. static DAEFunc::DAEJacFunc user_jac;
  55. static int nn;
  56.  
  57. DASSL::DASSL (void) : DAE ()
  58. {
  59.   stop_time_set = 0;
  60.   stop_time = 0.0;
  61.  
  62.   liw = 0;
  63.   lrw = 0;
  64.  
  65.   sanity_checked = 0;
  66.  
  67.   info.resize (15);
  68.  
  69.   for (int i = 0; i < 15; i++)
  70.     info.elem (i) = 0;
  71. }
  72.  
  73. DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f)
  74.   : DAE (state, time, f)
  75. {
  76.   n = size ();
  77.  
  78.   stop_time_set = 0;
  79.   stop_time = 0.0;
  80.  
  81.   liw = 20 + n;
  82.   lrw = 40 + 9*n + n*n;
  83.  
  84.   sanity_checked = 0;
  85.  
  86.   info.resize (15);
  87.  
  88.   for (int i = 0; i < 15; i++)
  89.     info.elem (i) = 0;
  90. }
  91.  
  92. DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv,
  93.       double time, DAEFunc& f)
  94.   : DAE (state, deriv, time, f)
  95. {
  96.   n = size ();
  97.  
  98.   stop_time_set = 0;
  99.   stop_time = 0.0;
  100.  
  101.   DAEFunc::set_function (f.function ());
  102.   DAEFunc::set_jacobian_function (f.jacobian_function ());
  103.  
  104.   liw = 20 + n;
  105.   lrw = 40 + 9*n + n*n;
  106.  
  107.   sanity_checked = 0;
  108.  
  109.   info.resize (15);
  110.  
  111.   for (int i = 0; i < 15; i++)
  112.     info.elem (i) = 0;
  113. }
  114.  
  115. void
  116. DASSL::force_restart (void)
  117. {
  118.   restart = 1;
  119.   integration_error = 0;
  120. }
  121.  
  122. void
  123. DASSL::set_stop_time (double t)
  124. {
  125.   stop_time_set = 1;
  126.   stop_time = t;
  127. }
  128.  
  129. void
  130. DASSL::clear_stop_time (void)
  131. {
  132.   stop_time_set = 0;
  133. }
  134.  
  135. int
  136. ddassl_f (const double& time, double *state, double *deriv,
  137.       double *delta, int& ires, double *, int *)
  138. {
  139.   ColumnVector tmp_deriv (nn);
  140.   ColumnVector tmp_state (nn);
  141.   ColumnVector tmp_delta (nn);
  142.  
  143.   for (int i = 0; i < nn; i++)
  144.     {
  145.       tmp_deriv.elem (i) = deriv [i];
  146.       tmp_state.elem (i) = state [i];
  147.     }
  148.  
  149.   tmp_delta = user_fun (tmp_state, tmp_deriv, time);
  150.  
  151.   if (tmp_delta.length () == 0)
  152.     ires = -2;
  153.   else
  154.     {
  155.       for (int i = 0; i < nn; i++)
  156.     delta [i] = tmp_delta.elem (i);
  157.     }
  158.  
  159.   return 0;
  160. }
  161.  
  162. int
  163. ddassl_j (const double& time, double *, double *, double *pd, const
  164.       double& cj, double *, int *)
  165. {
  166.   ColumnVector tmp_state (nn);
  167.   ColumnVector tmp_deriv (nn);
  168.  
  169.   // XXX FIXME XXX
  170.  
  171.   Matrix tmp_dfdxdot (nn, nn);
  172.   Matrix tmp_dfdx (nn, nn);
  173.  
  174.   DAEFunc::DAEJac tmp_jac;
  175.   tmp_jac.dfdxdot = &tmp_dfdxdot;
  176.   tmp_jac.dfdx    = &tmp_dfdx;
  177.  
  178.   tmp_jac = user_jac (tmp_state, tmp_deriv, time);
  179.  
  180.   // Fix up the matrix of partial derivatives for dassl.
  181.  
  182.   tmp_dfdx = tmp_dfdx +  cj * tmp_dfdxdot;
  183.  
  184.   for (int j = 0; j < nn; j++)
  185.     for (int i = 0; i < nn; i++)
  186.       pd [nn * j + i] = tmp_dfdx.elem (i, j);
  187.  
  188.   return 0;
  189. }
  190.  
  191. ColumnVector
  192. DASSL::do_integrate (double tout)
  193. {
  194.   ColumnVector retval;
  195.  
  196.   if (restart)
  197.     {
  198.       restart = 0;
  199.       info.elem (0) = 0;
  200.     }
  201.  
  202.   if (iwork.length () != liw)
  203.     iwork.resize (liw);
  204.  
  205.   if (rwork.length () != lrw)
  206.     rwork.resize (lrw);
  207.  
  208.   integration_error = 0;
  209.  
  210.   if (DAEFunc::jacobian_function ())
  211.     info.elem (4) = 1;
  212.   else
  213.     info.elem (4) = 0;
  214.  
  215.   double *px    = x.fortran_vec ();
  216.   double *pxdot = xdot.fortran_vec ();
  217.  
  218.   nn = n;
  219.   user_fun = DAEFunc::fun;
  220.   user_jac = DAEFunc::jac;
  221.  
  222.   if (! sanity_checked)
  223.     {
  224.       ColumnVector res = (*user_fun) (x, xdot, t);
  225.  
  226.       if (res.length () != x.length ())
  227.     {
  228.       (*current_liboctave_error_handler)
  229.         ("dassl: inconsistent sizes for state and residual vectors");
  230.  
  231.       integration_error = 1;
  232.       return retval;
  233.     }
  234.  
  235.       sanity_checked = 1;
  236.     }
  237.   
  238.   if (stop_time_set)
  239.     {
  240.       rwork.elem (0) = stop_time;
  241.       info.elem (3) = 1;
  242.     }
  243.   else
  244.     info.elem (3) = 0;
  245.  
  246.   double abs_tol = absolute_tolerance ();
  247.   double rel_tol = relative_tolerance ();
  248.  
  249.   if (initial_step_size () >= 0.0)
  250.     {
  251.       rwork.elem (2) = initial_step_size ();
  252.       info.elem (7) = 1;
  253.     }
  254.   else
  255.     info.elem (7) = 0;
  256.  
  257.   if (maximum_step_size () >= 0.0)
  258.     {
  259.       rwork.elem (1) = maximum_step_size ();
  260.       info.elem (6) = 1;
  261.     }
  262.   else
  263.     info.elem (6) = 0;
  264.  
  265.   double *dummy = 0;
  266.   int *idummy = 0;
  267.  
  268.   int *pinfo = info.fortran_vec ();
  269.   int *piwork = iwork.fortran_vec ();
  270.   double *prwork = rwork.fortran_vec ();
  271.  
  272. // again:
  273.  
  274.   F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo,
  275.                  rel_tol, abs_tol, idid, prwork, lrw,
  276.                  piwork, liw, dummy, idummy, ddassl_j));
  277.  
  278.   if (f77_exception_encountered)
  279.     {
  280.       integration_error = 1;
  281.       (*current_liboctave_error_handler) ("unrecoverable error in dassl");
  282.     }
  283.   else
  284.     {
  285.       switch (idid)
  286.     {
  287.     case 1: // A step was successfully taken in intermediate-output
  288.             // mode. The code has not yet reached TOUT.
  289.     case 2: // The integration to TSTOP was successfully completed
  290.             // (T=TSTOP) by stepping exactly to TSTOP.
  291.     case 3: // The integration to TOUT was successfully completed
  292.             // (T=TOUT) by stepping past TOUT.  Y(*) is obtained by
  293.             // interpolation.  YPRIME(*) is obtained by interpolation.
  294.  
  295.       retval = x;
  296.       t = tout;
  297.       break;
  298.  
  299.     case -1: // A large amount of work has been expended.  (~500 steps).
  300.     case -2: // The error tolerances are too stringent.
  301.     case -3: // The local error test cannot be satisfied because you
  302.              // specified a zero component in ATOL and the
  303.          // corresponding computed solution component is zero.
  304.          // Thus, a pure relative error test is impossible for
  305.          // this component.
  306.     case -6: // DDASSL had repeated error test failures on the last
  307.          // attempted step.
  308.     case -7: // The corrector could not converge.
  309.     case -8: // The matrix of partial derivatives is singular.
  310.     case -9: // The corrector could not converge.  There were repeated
  311.          // error test failures in this step.
  312.     case -10: // The corrector could not converge because IRES was
  313.           // equal to minus one.
  314.     case -11: // IRES equal to -2 was encountered and control is being
  315.           // returned to the calling program.
  316.     case -12: // DDASSL failed to compute the initial YPRIME.
  317.     case -33: // The code has encountered trouble from which it cannot
  318.           // recover. A message is printed explaining the trouble
  319.           // and control is returned to the calling program. For
  320.           // example, this occurs when invalid input is detected.
  321.     default:
  322.       integration_error = 1;
  323.       break;
  324.     }
  325.     }
  326.  
  327.   return retval;
  328. }
  329.  
  330. Matrix
  331. DASSL::do_integrate (const ColumnVector& tout)
  332. {
  333.   Matrix dummy;
  334.   return integrate (tout, dummy);
  335. }
  336.  
  337. Matrix
  338. DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out)
  339. {
  340.   Matrix retval;
  341.   int n_out = tout.capacity ();
  342.  
  343.   if (n_out > 0 && n > 0)
  344.     {
  345.       retval.resize (n_out, n);
  346.       xdot_out.resize (n_out, n);
  347.  
  348.       for (int i = 0; i < n; i++)
  349.     {
  350.       retval.elem (0, i) = x.elem (i);
  351.       xdot_out.elem (0, i) = xdot.elem (i);
  352.     }
  353.  
  354.       for (int j = 1; j < n_out; j++)
  355.     {
  356.       ColumnVector x_next = do_integrate (tout.elem (j));
  357.  
  358.       if (integration_error)
  359.         return retval;
  360.  
  361.       for (int i = 0; i < n; i++)
  362.         {
  363.           retval.elem (j, i) = x_next.elem (i);
  364.           xdot_out.elem (j, i) = xdot.elem (i);
  365.         }
  366.     }
  367.     }
  368.  
  369.   return retval;
  370. }
  371.  
  372. Matrix
  373. DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out,
  374.           const ColumnVector& tcrit) 
  375. {
  376.   Matrix retval;
  377.   int n_out = tout.capacity ();
  378.  
  379.   if (n_out > 0 && n > 0)
  380.     {
  381.       retval.resize (n_out, n);
  382.       xdot_out.resize (n_out, n);
  383.  
  384.       for (int i = 0; i < n; i++)
  385.     {
  386.       retval.elem (0, i) = x.elem (i);
  387.       xdot_out.elem (0, i) = xdot.elem (i);
  388.     }
  389.  
  390.       int n_crit = tcrit.capacity ();
  391.  
  392.       if (n_crit > 0)
  393.     {
  394.       int i_crit = 0;
  395.       int i_out = 1;
  396.       double next_crit = tcrit.elem (0);
  397.       double next_out;
  398.       while (i_out < n_out)
  399.         {
  400.           int do_restart = 0;
  401.  
  402.           next_out = tout.elem (i_out);
  403.           if (i_crit < n_crit)
  404.         next_crit = tcrit.elem (i_crit);
  405.  
  406.           int save_output;
  407.           double t_out;
  408.  
  409.           if (next_crit == next_out)
  410.         {
  411.           set_stop_time (next_crit);
  412.           t_out = next_out;
  413.           save_output = 1;
  414.           i_out++;
  415.           i_crit++;
  416.           do_restart = 1;
  417.         }
  418.           else if (next_crit < next_out)
  419.         {
  420.           if (i_crit < n_crit)
  421.             {
  422.               set_stop_time (next_crit);
  423.               t_out = next_crit;
  424.               save_output = 0;
  425.               i_crit++;
  426.               do_restart = 1;
  427.             }
  428.           else
  429.             {
  430.               clear_stop_time ();
  431.               t_out = next_out;
  432.               save_output = 1;
  433.               i_out++;
  434.             }
  435.         }
  436.           else
  437.         {
  438.           set_stop_time (next_crit);
  439.           t_out = next_out;
  440.           save_output = 1;
  441.           i_out++;
  442.         }
  443.  
  444.           ColumnVector x_next = do_integrate (t_out);
  445.  
  446.           if (integration_error)
  447.         return retval;
  448.  
  449.           if (save_output)
  450.         {
  451.           for (int i = 0; i < n; i++)
  452.             {
  453.               retval.elem (i_out-1, i) = x_next.elem (i);
  454.               xdot_out.elem (i_out-1, i) = xdot.elem (i);
  455.             }
  456.         }
  457.  
  458.           if (do_restart)
  459.         force_restart ();
  460.         }
  461.     }
  462.       else
  463.     {
  464.       retval = integrate (tout, xdot_out);
  465.  
  466.       if (integration_error)
  467.         return retval;
  468.     }
  469.     }
  470.  
  471.   return retval;
  472. }
  473.  
  474. /*
  475. ;;; Local Variables: ***
  476. ;;; mode: C++ ***
  477. ;;; End: ***
  478. */
  479.