Mercurial > hg > octave-max

changeset 3992:53b4eab68976

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-07-16 19:36:52 by jwe]
author jwe
date Tue, 16 Jul 2002 19:36:53 +0000
parents 48d2bc4a3729
children f23bc69132cc
files liboctave/ChangeLog liboctave/DASRT.cc liboctave/DASRT.h liboctave/DASSL.h src/ChangeLog src/DLD-FUNCTIONS/dasrt.cc
diffstat 6 files changed, 50 insertions(+), 116 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,5 +1,8 @@
 2002-07-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* DASRT.h (DASRT::set_ng, DASRT::get_ng): Delete
+	* DASRT.cc (DASRT::DASRT): Set ng here.
+
 	* DAEFunc.h (DAEFunc): Jacobian function now follows format of DASSL.
 	* DASSL.cc (ddassl_j): Make it work.
 	* DASPK.cc (ddaspk_j): Likewise.
--- a/liboctave/DASRT.cc
+++ b/liboctave/DASRT.cc
@@ -163,13 +163,14 @@
   info.resize (30, 0);
 
   npar = 0;
+  ng = 0;
 
   liw = 0;
   lrw = 0;
 }
 
-DASRT::DASRT (const int& ng, const ColumnVector& state, 
-	      const ColumnVector& deriv, double time, DAERTFunc& f)
+DASRT::DASRT (const ColumnVector& state, const ColumnVector& deriv,
+	      double time, DAERTFunc& f)
   : DAERT (state, deriv, time, f)
 {
   n = size ();
@@ -180,8 +181,6 @@
   stop_time_set = false;
   stop_time = 0.0;
 
-  DAERTFunc::operator = (f);
-
   sanity_checked = false;
 
   info.resize (30, 0);
@@ -189,6 +188,16 @@
 
   npar = 0;
 
+  DAERTFunc::DAERTConstrFunc tmp_csub = DAERTFunc::constraint_function ();
+  
+  if (tmp_csub)
+    {
+      ColumnVector tmp = tmp_csub (state, time);
+      ng = tmp.length ();
+    }
+  else
+    ng = 0;
+
   rpar.resize (npar+1);
   ipar.resize (npar+1);
 
@@ -226,7 +235,6 @@
   double *prwork = rwork.fortran_vec ();
   int *piwork = iwork.fortran_vec ();
 
-
   F77_FUNC (ddasrt, DASRT) (ddasrt_f, n, t, py, pydot, t, pinfo,
 			    &rel_tol, &abs_tol, idid, prwork, lrw,
 			    piwork, liw, prpar, pipar, ddasrt_j,
@@ -234,7 +242,6 @@
 
   int iwadd = iwork(18);
 
-
   if (iwadd > 0)
     liw += iwadd;
 
@@ -251,14 +258,12 @@
 
   int rwadd = iwork(19);
 
-
   if (rwadd > 0)
     lrw += rwadd;
 
   rwork.resize (lrw, 0.0);
 
   info(0) = info_zero;
-
 }
 
 void
@@ -276,17 +281,6 @@
 }
 
 void
-DASRT::set_ng (int the_ng)
-{
-  ng = the_ng;
-}
-
-int DASRT::get_ng (void)
-{
-  return ng;
-}
-
-void
 DASRT::clear_stop_time (void)
 {
   stop_time_set = false;
@@ -302,10 +296,10 @@
       info(0) = 0;
 
       for (int i = 0; i < n; i++)
-       {
-	y(i,0) = x(i);
-        ydot(i,0) = xdot(i);
-       }
+	{
+	  y(i,0) = x(i);
+	  ydot(i,0) = xdot(i);
+	}
 
       integration_error = false;
 
@@ -335,13 +329,9 @@
 
 	  sanity_checked = true;
 	}
-
   
       init_work_size (info(0));
 
-
-
-
       if (iwork.length () != liw)
 	iwork.resize (liw);
 
@@ -351,7 +341,6 @@
       abs_tol = absolute_tolerance ();
       rel_tol = relative_tolerance ();
 
-
       if (initial_step_size () >= 0.0)
 	{
 	  rwork(2) = initial_step_size ();
@@ -376,7 +365,6 @@
       else
 	info(6) = 0;
 
-
       py = y.fortran_vec ();
       pydot = ydot.fortran_vec ();
       pinfo = info.fortran_vec ();
@@ -404,9 +392,6 @@
 	info(3) = 0;
     }
 
-
-
-
   F77_XFCN (ddasrt, DASRT, (ddasrt_f, n, t, py, pydot, tout, pinfo,
 			    &rel_tol, &abs_tol, idid, prwork, lrw,
 			    piwork, liw, prpar, pipar, ddasrt_j,
@@ -486,7 +471,6 @@
 
   int n_out = tout.capacity ();
 
-
   if (n_out > 0 && n > 0)
     {
       x_out.resize (n_out, n);
@@ -501,25 +485,26 @@
 	      retval = DASRT_result (x_out, xdot_out, t_out);
 	      return retval;
 	    }
+
           if (idid == 4)
             t_out(j) = t;
           else
             t_out(j) = tout(j);
 
-
 	  for (int i = 0; i < n; i++)
 	    {
 	      x_out(j,i) = y(i,0);
 	      xdot_out(j,i) = ydot(i,0);
 	    }
-          if (idid ==4)
-           {
-            oldj = j;
-            j = n_out;
-            x_out.resize (oldj+1, n);
-            xdot_out.resize (oldj+1, n);
-            t_out.resize (oldj+1);
-           }
+
+          if (idid == 4)
+	    {
+	      oldj = j;
+	      j = n_out;
+	      x_out.resize (oldj+1, n);
+	      xdot_out.resize (oldj+1, n);
+	      t_out.resize (oldj+1);
+	    }
 	}
     }
 
@@ -606,6 +591,7 @@
 		  retval = DASRT_result (x_out, xdot_out, t_outs);
 		  return retval;
 		}
+
               if (idid == 4)
                 t_out = t;
 
@@ -624,7 +610,6 @@
                       t_outs.resize (i_out);
                       i_out = n_out;
                     }
-
 		}
 
 	      if (do_restart)
--- a/liboctave/DASRT.h
+++ b/liboctave/DASRT.h
@@ -158,7 +158,7 @@
 
   DASRT (void);
 
-  DASRT (const int& ng, const ColumnVector& x, const ColumnVector& xdot, 
+  DASRT (const ColumnVector& state, const ColumnVector& deriv,
 	 double time, DAERTFunc& f);
 
   ~DASRT (void) { }
@@ -167,8 +167,6 @@
 
   void set_stop_time (double t);
   void clear_stop_time (void);
-  void set_ng (int the_ng);
-  int get_ng (void);
 
   DASRT_result integrate (const ColumnVector& tout);
 
--- a/liboctave/DASSL.h
+++ b/liboctave/DASSL.h
@@ -113,9 +113,9 @@
 
   DASSL (void);
 
-  DASSL (const ColumnVector& x, double time, DAEFunc& f);
+  DASSL (const ColumnVector& state, double time, DAEFunc& f);
 
-  DASSL (const ColumnVector& x, const ColumnVector& xdot,
+  DASSL (const ColumnVector& state, const ColumnVector& xdot,
 	 double time, DAEFunc& f);
 
   ~DASSL (void) { }
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,5 +1,7 @@
 2002-07-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* DLD-FUNCTIONS/dasrt.cc (Fdasrt): No need to find ng here.
+
 	* DLD-FUNCTIONS/dassl.cc (dassl_user_jacobian): New function.
 	(Fdassl): Handle Jacobian function.
 
--- a/src/DLD-FUNCTIONS/dasrt.cc
+++ b/src/DLD-FUNCTIONS/dasrt.cc
@@ -327,7 +327,6 @@
     DASRT_ABORT1 ("invalid recursive call");
 
   int argp = 0;
-  int ng = 0;
 
   int nargin = args.length ();
 
@@ -393,76 +392,25 @@
   DAERTFunc func (dasrt_user_f);
   
   argp++;
-  double t0;
-  ColumnVector x0;
   
   if (args(1).is_string ())
     {
       dasrt_cf = is_valid_function (args(1), "dasrt", true);
 
       if (! dasrt_cf)
-	  DASRT_ABORT1 ("expecting function name as argument 2");
-      else
-        {
-	  octave_value_list blah;
-	  octave_value_list inputs;
-	  ColumnVector the_ts;
-	  x0 = ColumnVector (args(2).vector_value ());
-
-	  if (error_state)
-	    DASRT_ABORT1 ("bad x0");
-
-	  the_ts = ColumnVector (args(4).vector_value ());
-
-	  if (error_state)
-	    DASRT_ABORT1 ("bad tout");
-
-	  t0 = the_ts(0);
+	DASRT_ABORT1 ("expecting function name as argument 2");
 
-	  inputs(0) = x0;
-	  inputs(1) = t0;
-
-	  string g_name = args(1).string_value();
-
-	  if (error_state)
-	    DASRT_ABORT1 ("bad function name of stopping condition");
-
-	  blah = feval (g_name, inputs, 2);
-	  ColumnVector gout = ColumnVector (blah(0).vector_value ());
-
-	  ng = gout.length();
-	  int testg;
-	  for (testg = 0; testg < ng; testg++)
-	    {
-	      if (gout(testg) == 0)
-		DASRT_ABORT1 ("stopping condition satisfied at initial point");
-	    }
-        }
       argp++;
 
       func.set_constraint_function (dasrt_user_cf);
     }
-  else
-    {
-      // Now this second argument is not a string.  It has to be x0.
-      // we call the dummy g function and set ng = 1;
-
-      x0 = ColumnVector (args(1).vector_value ());
-
-      if (error_state)
-	DASRT_ABORT1 ("bad x0");
-
-      func.set_constraint_function (dasrt_dumb_cf);
-
-      ng = 1;
-    }
 
   ColumnVector state (args(argp++).vector_value ());
 
   if (error_state)
     DASRT_ABORT2 ("expecting state vector as argument %d", argp);
 
-  ColumnVector stateprime (args(argp++).vector_value());
+  ColumnVector stateprime (args(argp++).vector_value ());
 
   if (error_state)
     DASRT_ABORT2 
@@ -477,12 +425,11 @@
   double tzero = old_out_times (0);
 
   int ol = old_out_times.length ();
-  int ijk;
 
   ColumnVector out_times (ol-1, 0.0);
 
-  for (ijk = 1; ijk < ol; ijk++)
-    out_times(ijk-1) = old_out_times(ijk);
+  for (int i = 1; i < ol; i++)
+    out_times(i-1) = old_out_times(i);
 
   ColumnVector crit_times;
 
@@ -506,10 +453,9 @@
 
   DASRT_result output;
 
-  DASRT dae = DASRT (ng, state, stateprime, tzero, func);
+  DASRT dae = DASRT (state, stateprime, tzero, func);
 
   dae.copy (dasrt_opts);
-  dae.set_ng (ng);
 
   if (error_state)
     DASRT_ABORT1 ("something is wrong");
@@ -527,7 +473,7 @@
       Matrix old_output_state = output.state ();
 
       int lstuff = old_output_times.length ();
-      int lstate = x0.length ();
+      int lstate = state.length ();
 
       ColumnVector output_times (lstuff+1, 0.0);
 
@@ -536,20 +482,20 @@
 
       output_times(0) = tzero;
 
-      for (ijk = 0; ijk < lstate; ijk++)
+      for (int i = 0; i < lstate; i++)
 	{
-	  output_deriv(0,ijk) = stateprime(ijk);
-	  output_state(0,ijk) = state(ijk);
+	  output_deriv(0,i) = stateprime(i);
+	  output_state(0,i) = state(i);
 	}
 
-      for (ijk = 0; ijk < lstuff; ijk++)
+      for (int i = 0; i < lstuff; i++)
 	{
-	  output_times(ijk+1) = old_output_times(ijk);
+	  output_times(i+1) = old_output_times(i);
 
-	  for (int lmnop=0; lmnop < lstate; lmnop++)
+	  for (int j = 0; j < lstate; j++)
 	    {
-	      output_deriv(ijk+1,lmnop) = old_output_deriv(ijk,lmnop);
-	      output_state(ijk+1,lmnop) = old_output_state(ijk,lmnop);
+	      output_deriv(i+1,j) = old_output_deriv(i,j);
+	      output_state(i+1,j) = old_output_state(i,j);
 	    }
 	}