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