Mercurial > forge
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) */