changeset 3100:0fce7ea511ec octave-forge

Finished implementation of the interface function for stiff-problem (ODE/DAE) solver radau.
author treichl
date Sat, 10 Feb 2007 18:18:29 +0000
parents 453456a6c246
children 3e4ed9d2b97a
files main/odepkg/src/odepkg_mexsolver_radau.c
diffstat 1 files changed, 11 insertions(+), 9 deletions(-) [+]
line wrap: on
line diff
--- a/main/odepkg/src/odepkg_mexsolver_radau.c	Sat Feb 10 17:25:41 2007 +0000
+++ b/main/odepkg/src/odepkg_mexsolver_radau.c	Sat Feb 10 18:18:29 2007 +0000
@@ -152,8 +152,9 @@
   }
 }
 
-void F77_FUNC (fmas, FMAS) (int *N, double *X, double *Y, double *AM, int *LMAS,
-  double *RPAR, int *IPAR) {
+void F77_FUNC (fmas, FMAS) (int *N, double *X, double *Y, double *AM, 
+  GCC_ATTR_UNUSED int *LMAS, GCC_ATTR_UNUSED double *RPAR,
+  GCC_ATTR_UNUSED int *IPAR) {
 
   int  vcnt = 0;
   int  vnum = 0;
@@ -553,8 +554,9 @@
   fodepkgvar (2, "OdeOptions", &vtmp);
   fodepkgvar (2, "DefaultOptions", &vtem);
   if (!mxIsEqual (mxGetField (vtmp, 0, "Jacobian"), mxGetField (vtem, 0, "Jacobian"))) {
-    IJAC = 1;  /* Tell the solver that we have a Jacobian matrix */
-    MLJAC = N; /* Tell the solver that the matrix is full */
+    IJAC = 1;     /* Tell the solver that we have a Jacobian matrix */
+    MLJAC = N;    /* Tell the solver that the matrix is full */
+    WORK[2] = -1; /* Tell the solver to recompute Jacobian after every succesful step */
     vtem = mxGetField (vtmp, 0, "Jacobian");
     if (mxGetClassID (vtem) == mxFUNCTION_CLASS) { /* function handle */
       if (mexCallMATLAB (1, &vtmp, 1, &vtem, "func2str"))
@@ -711,7 +713,7 @@
   mexPrintf ("ODEPKGDEBUG: ----- STARTING SOLVER CALCULATION PROCEDURE\n");
 #endif
 
-  IWORK[00] = 0;          /* Switch for transformation of Jacobian into Hessenberg form */
+  IWORK[00] = 1;          /* Switch for transformation of Jacobian into Hessenberg form */
   /* IWORK[01] = 0; */    /* Maximal number of allowed steps */
   /* IWORK[02] = 7; */    /* Maximal number of newton iterations (default value = 7) */
   /* IWORK[03] = 0; */    /* Take extrapolated collocation solution at startup */
@@ -731,7 +733,7 @@
   /* FORTRAN SUBROUTINE RADAU(N,FCN,X,Y,XEND,H,RTOL,ATOL,ITOL,JAC,IJAC,MLJAC,MUJAC, */
   /*   MAS,IMAS,MLMAS,MUMAS,SOLOUT,IOUT,WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) */
 
-  vnum = 2; /* Needed to call output at every successful step */
+  vnum = 1; /* Needed to call output at every successful step */
   F77_FUNC (radau, RADAU) (&N, &F77_FUNC (ffcn, FFCN), &SLOT[0],
     INIT, &SLOT[1], &H, RTOL, ATOL, &ITOL,
     &F77_FUNC (fjac, FJAC), &IJAC, &MLJAC, &MUJAC,
@@ -744,7 +746,7 @@
 
   switch (VRET) {
     case -4:
-      mexPrintf ("The computation has been stopped because the problem is probably stiff\n");
+      mexPrintf ("The matrix is repeatedly singular\n");
       break;
     case -3:
       mexPrintf ("The step size grew too small, reduce InitialStep and/or MaxStep\n");
@@ -757,7 +759,7 @@
       break;
     case 0:  break; /* Not treated */
     case 1:  break; /* Computation has been successful */
-    case 2:  break; /* Computation has been successful, stopped by fsolout */
+    case 2:  break; /* Computation has been successful, stopped by fsol */
     default: break;
   }
 #ifdef __ODEPKGDEBUG__
@@ -801,7 +803,7 @@
     vnum = IWORK[19]; /* A dopri solver result */
     vtem = mxCreateDoubleScalar ((double) vnum);
     fodepkgvar (1, "vlinsol", &vtem);
-    mexPrintf ("Number of fb-substitutions:  %d\n", vnum);
+    mexPrintf ("Number of forw/back. subst.: %d\n", vnum);
   }
 
   if (nlhs == 1) { /* Handle the PLHS array (1 output argument) */