lsoda.c

GESPI 2.0动态系统模拟工具

/*
From [e-mail] Wed Apr 24 01:35:52 1991
Return-Path: <tam>
Date: Wed, 24 Apr 91 03:35:24 CDT
From: [e-mail]
To: [e-mail]
Subject: lsoda.c
Cc: [e-mail]


I'm told by Steve Nichols at Georgia Tech that you are interested in
a stiff integrator.  Here's a translation of the fortran code LSODA.

Please note
that there is no comment.  The interface is the same as the FORTRAN
code and I believe the documentation in LSODA will suffice.
As usual, a free software comes with no guarantee.

Hon Wah Tam
Wolfram Research, Inc.
[e-mail]
*/

/*
I have done some additions to lsoda.c . These were mainly to fill the
gap of some features that were available in the fortran code and were
missing in the C version.

Changes are:
	     * all messages printed by lsoda routines will start with
	       a '\n' (instead of ending with one).
	     * xsetf() was added so that user can control printing
	       of messages by lsoda. xsetf should be called before
	       calling lsoda: xsetf(0) switches printing off
			      xsetf(1) swithces printing on (default)
	       this implies one new global variable prfl (print flag).
	       xsetf(0) will stop *any* printing by lsoda functions.
	       Turning printing off means losing valuable information
	       but will not scramble your stderr output ;-)
	       This function is part of the original FORTRAN version.
	     * xsetf() and intdy() are now callable from outside the
	       library as assumed in the FORTRAN version.
	     * created lsoda.h that should be included in blocks
	       calling functions in lsoda's library.
	     * iwork5 can now have an extra value:
		      0 - no extra printing (default)
		      1 - print data on each method switch
		  ->  2 - print useful information at *each* stoda
			  step (one lsoda call has performs many stoda
			  calls) and also data on each method switch
	       note that a xsetf(0) call will prevent any printing even
	       if iwork5 > 0.
	     * hu, tn were made available as extern variables.
	     * the variable _ml has been renamed to _ml

					Pedro Mendes

		      eMail:      INTERNET: [e-mail]
		      sMail:                Pedro Mendes,
					    Dept. of Biological Sciences,
					    University of Wales, Aberystwyth
					    Dyfed, SY23 3DA,
					    United Kingdom.

		      phone: international: + 44 970 622353
				      U.K.:     0970 622353
*/

#include <stdio.h>
#include <math.h>
/*#include <conio.h>
#include "globals.h"
#include "globvar.h"
#include "datab.h"*/
#ifndef _ZTC
#include <malloc.h>
#define mem_malloc malloc
#define mem_free free
#define mem_realloc realloc
#else
#include <stdlib.h>
#define MEM_DEBUG 1
#include "mem.h"
#endif

#define max( a , b )  ( (a) > (b) ? (a) : (b) )
#define min( a , b )  ( (a) < (b) ? (a) : (b) )

#define ETA 2.2204460492503131e-16

/* functions defined in other blocks */
extern void daxpy(), dgefa(), dgesl(), dscal();
extern double ddot();
extern int idamax();

/* functions accessible from other blocks */
void xsetf( int  mflag );
void intdy();

/* internal functions */
static void terminate( int *istate );
static void terminate2( double *y, double *t );
void freevectors( void );
static void successreturn( double *y, double *t, int itask,
			   int ihit, double tcrit, int *istate );
static void stoda( int neq, double *y, void (*f)() );
static void endstoda( void );
static void prja( int neq, double *y, void (*f)() );
static void correction( int neq, double *y, void (*f)(), int *corflag,
			double pnorm, double *del, double *delp,
			double *told, int *ncf, double *rh, int *m );
static void resetcoeff( void );
static void methodswitch( double dsm, double pnorm, double *pdh, double *rh );

static void
   solsy(),
   cfode(),
   ewset(),
   scaleh(),
   orderswitch(),
   corfailure();

static double
   vmnorm(), bnorm(), fnorm();

/* new static variable to enable optional printing (see comments) */
static int prfl = 1;

/* newly added static variables */
/* _ml was originally ml - changed to avoid conflict with one external ml */
static int _ml, mu;
static mord[3] = { 0, 12, 5 };
static double sqrteta, *yp1, *yp2;
static double sm1[13] = { 0., 0.5, 0.575, 0.55, 0.45, 0.35, 0.25,
			      0.2, 0.15, 0.1, 0.075, 0.05, 0.025 };

/* variables accessible from other routines */
double hu, tn, h, tsw, tolsf;
int nst, nfe, nje, nqu, nq, imxer, mused, meth;

/* static variables for lsoda() */
static double ccmax, el0, hmin, hmxi, rc;
static int illin = 0, init = 0, mxstep, mxhnil, nhnil, ntrep = 0,
	   nslast, nyh, ierpj, iersl, jcur, jstart, kflag, l,
	   miter, maxord, maxcor, msbp, mxncf, n;
static double pdnorm;
static int ixpr = 0, jtyp, mxordn, mxords;

/* no static variable for prja(), solsy() */
/* static variables for stoda() */

static double conit, crate, el[14], elco[13][14], hold, rmax,
       tesco[13][4];
static int ialth, ipup, lmax, meo, nslp;
static double pdest, pdlast, ratio, cm1[13], cm2[6];
static int icount, irflag;

/* static variable for block data */
static int mesflg = 1;

/* static variables for various vectors and the Jacobian. */
static double **yh, **wm, *ewt, *savf, *acor;
static int *ipvt;

/*
   The following are useful statistics.

   hu,
   h,
   tn,
   tolsf,
   tsw,
   nst,
   nfe,
   nje,
   nqu,
   nq,
   imxer,
   mused,
   meth
*/

/*
   free allocated vectors
*/

void freevectors( void )
{
 int i, lenyh;

   lenyh = 1 + max( mxordn, mxords );
   for( i=1; i<=lenyh; i++ ) mem_free( (void *) yh[i] );
   for( i=1; i<=nyh; i++ ) mem_free( (void *) wm[i] );
   mem_free( (void *) yh );
   mem_free( (void *) wm );
   mem_free( (void *) ewt );
   mem_free( (void *) savf );
   mem_free( (void *) acor );
   mem_free( (void *) ipvt );
}     /*   end freevectors   */

/*
   Terminate lsoda due to illegal input.
*/

static void terminate( int *istate )
{
 if ( ( illin == 5 ) && prfl )
 {
  printf( "\nlsoda -- repeated occurrence of illegal input\n" );
  printf(   "         run aborted.. apparent infinite loop" );
 }
 else
 {
  illin++;
  *istate = -3;
 }
}         /*   end terminate   */


/*
   Terminate lsoda due to various error conditions.
*/


static void terminate2( double *y, double *t )
{
 int i;

 yp1 = yh[1];
 for ( i = 1 ; i <= n ; i++ ) y[i] = yp1[i];
 *t = tn;
 illin = 0;
 freevectors();
 return;
}         /*   end terminate2   */


/*
   The following block handles all successful returns from lsoda.
   If itask != 1, y is loaded from yh and t is set accordingly.
   *Istate is set to 2, the illegal input counter is zeroed, and the
   optional outputs are loaded into the work arrays before returning.
*/

static void successreturn( double *y, double *t, int itask,
			   int ihit, double tcrit, int *istate )
{
 int i;

 yp1 = yh[1];
 for ( i = 1 ; i <= n ; i++ ) y[i] = yp1[i];
 *t = tn;
 if ( itask == 4 || itask == 5 )
  if ( ihit )
   *t = tcrit;
 *istate = 2;
 illin = 0;
 freevectors();
}   /*   end successreturn   */


/*
   In this version all memory allocated using malloc is freed upon exit.
   Therefore *istate = 2 and *istate = 3 will not work.
*/

void lsoda( f, neq, y, t, tout, itol, rtol, atol, itask, istate,
	    iopt, jt, iwork1, iwork2, iwork5, iwork6, iwork7, iwork8,
	    iwork9, rwork1, rwork5, rwork6, rwork7 )

void (*f)();
int neq, itol, itask, *istate, iopt, jt;
int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
double *y, *t, tout, *rtol, *atol;
double rwork1, rwork5, rwork6, rwork7;

/*
   If the user does not supply any of these values, the calling program
   should initialize those untouched working variables to zero.

   _ml = iwork1
   mu = iwork2
   ixpr = iwork5
   mxstep = iwork6
   mxhnil = iwork7
   mxordn = iwork8
   mxords = iwork9

   tcrit = rwork1
   h0 = rwork5
   hmax = rwork6
   hmin = rwork7
*/


{
   int mxstp0 = 500, mxhnl0 = 10;

   int i, i1, i2, iflag, kgo, lf0, lenyh, ihit;
   double atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
	  tcrit, tdist, tnext, tol, tp, size, sum, w0;

/*
   Block a.
   This code block is executed on every call.
   It tests *istate and itask for legality and branches appropriately.
   If *istate > 1 but the flag init shows that initialization has not
   yet been done, an error return occurs.
   If *istate = 1 and tout = t, return immediately.
*/

   if ( *istate < 1 || *istate > 3 ) {
      if ( prfl )
	 printf( "\nlsoda -- illegal istate = %d", *istate );
      terminate( istate );
      return;
   }
   if ( itask < 1 || itask > 5 ) {
      if ( prfl )
	 printf( "\nlsoda -- illegal itask = %d", itask );
      terminate( istate );
      return;
   }
   if ( init == 0 && ( *istate == 2 || *istate == 3 ) ) {
      if ( prfl )
	 printf( "\nlsoda -- istate > 1 but lsoda not initialized" );
      terminate( istate );
      return;
   }
   if ( *istate == 1 ) {
      init = 0;
      if ( tout == *t ) {
	 ntrep++;
	 if ( ntrep < 5 )
	    return;
	 if ( prfl ) {
	    printf( "\nlsoda -- repeated calls with istate = 1 and tout = t" );
	    printf( "\n         run aborted.. apparent infinite loop" );
	 }
	 return;
      }
   }

/*
   Block b.
   The next code block is executed for the initial call ( *istate = 1 ),
   or for a continuation call with parameter changes ( *istate = 3 ).
   It contains checking of all inputs and various initializations.

   First check legality of the non-optional inputs neq, itol, iopt,
   jt, _ml, and mu.
*/

   if ( *istate == 1 || *istate == 3 ) {
      ntrep = 0;
      if ( neq <= 0 ) {
	 if ( prfl )
	    printf( "\nlsoda -- neq = %d is less than 1", neq );
	 terminate( istate );
	 return;
      }
      if ( *istate == 3 && neq > n ) {
	 if ( prfl )
	    printf( "\nlsoda -- istate = 3 and neq increased" );
	 terminate( istate );
	 return;
      }
      n = neq;
      if ( itol < 1 || itol > 4 ) {
	 if ( prfl )
	    printf( "\nlsoda -- itol = %d illegal", itol );
	 terminate( istate );
	 return;
      }
      if ( iopt < 0 || iopt > 1 ) {
	 if ( prfl )
	    printf( "\nlsoda -- iopt = %d illegal", iopt );
	 terminate( istate );
	 return;
      }
      if ( jt == 3 || jt < 1 || jt > 5 ) {
	 if ( prfl )
	    printf( "\nlsoda -- jt = %d illegal", jt );
	 terminate( istate );
	 return;
      }
      jtyp = jt;
      if ( jt > 2 ) {
	 _ml = iwork1;
	 mu = iwork2;
	 if ( _ml < 0 || _ml >= n ) {
	    if ( prfl )
	       printf( "\nlsoda -- _ml = %d not between 1 and neq", _ml );
	    terminate( istate );
	    return;
	 }
	 if ( mu < 0 || mu >= n ) {
	    if ( prfl )
	       printf( "\nlsoda -- mu = %d not between 1 and neq", mu );
	    terminate( istate );
	    return;
	 }
      }

/* Next process and check the optional inpus.   */

/* Default options.   */

      if ( iopt == 0 ) {
	 ixpr = 0;
	 mxstep = mxstp0;
	 mxhnil = mxhnl0;
	 hmxi = 0.;
	 hmin = 0.;
	 if ( *istate == 1 ) {
	    h0 = 0.;
	    mxordn = mord[1];
	    mxords = mord[2];
	 }
      }        /*   end if ( iopt == 0 )   */

/* Optional inputs.   */

      else {             /*   if ( iopt = 1 )  */
	 ixpr = iwork5;
	 if ( ixpr < 0 || ixpr > 2 ) {
	    if ( prfl )
	       printf( "\nlsoda -- ixpr = %d is illegal", ixpr );
	    terminate( istate );
	    return;
	 }
	 mxstep = iwork6;
	 if ( mxstep < 0 ) {
	    if ( prfl )
	       printf( "\nlsoda -- mxstep < 0" );
	    terminate( istate );
	    return;
	 }
	 if ( mxstep == 0 )
	    mxstep = mxstp0;
	 mxhnil = iwork7;
	 if ( mxhnil < 0 ) {
	    if ( prfl )
	       printf( "\nlsoda -- mxhnil < 0" );
	    terminate( istate );
	    return;
	 }
	 if ( *istate == 1 ) {
	    h0 = rwork5;
	    mxordn = iwork8;
	    if ( mxordn < 0 ) {
	       if ( prfl )
		  printf( "\nlsoda -- mxordn = %d is less than 0", mxordn );
	       terminate( istate );
	       return;
	    }
	    if ( mxordn == 0 )
	       mxordn = 100;
	    mxordn = min( mxordn, mord[1] );
	    mxords = iwork9;
	    if ( mxords < 0 ) {
	       if ( prfl )
		  printf( "\nlsoda -- mxords = %d is less than 0", mxords );
	       terminate( istate );
	       return;
	    }
	    if ( mxords == 0 )