lsoda.c

GESPI 2.0动态系统模拟工具

	       mxords = 100;
	    mxords = min( mxords, mord[2] );
	    if ( ( tout - *t ) * h0 < 0. ) {
	       if ( prfl ) {
		  printf( "\nlsoda -- tout = %g behind t = %g", tout, *t );
		  printf( "\n         integration direction is given by %g",
			   h0 );
	       }
	       terminate( istate );
	       return;
	    }
	 }         /*  end if ( *istate == 1 )  */
	 hmax = rwork6;
	 if ( hmax < 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- hmax < 0." );
	    terminate( istate );
	    return;
	 }
	 hmxi = 0.;
	 if ( hmax > 0 )
	    hmxi = 1. / hmax;
	 hmin = rwork7;
	 if ( hmin < 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- hmin < 0." );
	    terminate( istate );
	    return;
	 }
      }      /*   end else   */      /*   end iopt = 1   */
   }    /*   end if ( *istate == 1 || *istate == 3 )   */
/*
   If *istate = 1, meth is initialized to 1.

   Also allocate memory for yh, wm, ewt, savf, acor, ipvt.
*/
   if ( *istate == 1 ) {
/*
   If memory were not freed, *istate = 3 need not reallocate memory.
   Hence this section is not executed by *istate = 3.
*/
      sqrteta = sqrt( ETA );
      meth = 1;
      nyh = n;
      lenyh = 1 + max( mxordn, mxords );

      yh = ( double ** ) mem_malloc( ( 1 + lenyh ) * sizeof( *yh ) );
      if ( yh == NULL ) {
	 if ( prfl )
	    printf( "\nlsoda -- insufficient memory for your problem" );
	 terminate( istate );
	 return;
      }
      for ( i = 1 ; i <= lenyh ; i++ )
	 yh[i] = ( double * ) mem_malloc( ( 1 + nyh ) * sizeof( double ) );

      wm = ( double ** ) mem_malloc( ( 1 + nyh ) * sizeof( *wm ) );
      if ( wm == NULL ) {
	 mem_free( yh );
	 if ( prfl )
	    printf( "\nlsoda -- insufficient memory for your problem" );
	 terminate( istate );
	 return;
      }
      for ( i = 1 ; i <= nyh ; i++ )
	 wm[i] = ( double * ) mem_malloc( ( 1 + nyh ) * sizeof( double ) );

      ewt = ( double * ) mem_malloc( ( 1 + nyh ) * sizeof( double ) );
      if ( ewt == NULL ) {
	 mem_free( yh );
	 mem_free( wm );
	 if ( prfl )
	    printf( "\nlsoda -- insufficient memory for your problem" );
	 terminate( istate );
	 return;
      }

      savf = ( double * ) mem_malloc( ( 1 + nyh ) * sizeof( double ) );
      if ( savf == NULL ) {
	 mem_free( yh );
	 mem_free( wm );
	 mem_free( ewt );
	 if ( prfl )
	    printf( "\nlsoda -- insufficient memory for your problem" );
	 terminate( istate );
	 return;
      }

/*
printf( "\nLSODA Block B - before acor alloc");
printf( "\n -  low = %lf, ampl = %lf, dens = %d", sparam[24].low, sparam[24].ampl, sparam[24].dens);
printf("\nlow = %Fp", &sparam[24].low);
printf("\nacor = %Fp", acor);
*/
      acor = ( double * ) mem_malloc( ( 1 + nyh ) * sizeof( double ) );
      if ( acor == NULL ) {
	 mem_free( yh );
	 mem_free( wm );
	 mem_free( ewt );
	 mem_free( savf );
	 if ( prfl )
	    printf( "\nlsoda -- insufficient memory for your problem" );
	 terminate( istate );
	 return;
      }

      ipvt = ( int * ) mem_malloc( ( 1 + nyh ) * sizeof( int ) );
      if ( ipvt == NULL ) {
	 mem_free( yh );
	 mem_free( wm );
	 mem_free( ewt );
	 mem_free( savf );
	 mem_free( acor );
	 if ( prfl )
	    printf( "\nlsoda -- insufficient memory for your problem" );
	 terminate( istate );
	 return;
      }
   }
/*
   Check rtol and atol for legality.
*/
   if ( *istate == 1 || *istate == 3 ) {
      rtoli = rtol[1];
      atoli = atol[1];
      for ( i = 1 ; i <= n ; i++ ) {
	 if ( itol >= 3 )
	    rtoli = rtol[i];
	 if ( itol == 2 || itol == 4 )
	    atoli = atol[i];
	 if ( rtoli < 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- rtol = %g is less than 0.", rtoli );
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 if ( atoli < 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- atol = %g is less than 0.", atoli );
	    terminate( istate );
	    freevectors();
	    return;
	 }
      }     /*   end for   */
   }   /*   end if ( *istate == 1 || *istate == 3 )   */
/*
   If *istate = 3, set flag to signal parameter changes to stoda.
*/
   if ( *istate == 3 ) {
      jstart = -1;
   }

/*
   Block c.
   The next block is for the initial call only ( *istate = 1 ).
   It contains all remaining initializations, the initial call to f,
   and the calculation of the initial step size.
   The error weights in ewt are inverted after being loaded.
*/
   if ( *istate == 1 ) {
      tn = *t;
      tsw = *t;
      maxord = mxordn;
      if ( itask == 4 || itask == 5 ) {
	 tcrit = rwork1;
	 if ( ( tcrit - tout ) * ( tout - *t )  < 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- itask = 4 or 5 and tcrit behind tout" );
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 if ( h0 != 0. && ( *t + h0 - tcrit ) * h0 > 0. )
	    h0 = tcrit - *t;
      }
      jstart = 0;
      nhnil = 0;
      nst = 0;
      nje = 0;
      nslast = 0;
      hu = 0.;
      nqu = 0;
      mused = 0;
      miter = 0;
      ccmax = 0.3;
      maxcor = 3;
      msbp = 20;
      mxncf = 10;
/*
   Initial call to f.
*/
      (*f)( neq, *t, y, yh[2] );
      nfe = 1;
/*
   Load the initial value vector in yh.
*/
      yp1 = yh[1];
      for ( i = 1 ; i <= n ; i++)
	 yp1[i] = y[i];
/*
   Load and invert the ewt array.  ( h is temporarily set to 1. )
*/
      nq = 1;
      h = 1.;
      ewset( itol, rtol, atol, y );
      for ( i=1 ; i<=n ; i++ )
      {
       if ( ewt[i]<=0. )
       {
	if ( prfl )
	 printf( "\nlsoda -- ewt[%d] = %g <= 0.", i, ewt[i] );
	terminate( istate );
	return;
       }
       ewt[i] = 1. / ewt[i];
      }

/*
   The coding below computes the step size, h0, to be attempted on the
   first step, unless the user has supplied a value for this.
   First check that tout - *t differs significantly from zero.
   A scalar tolerance quantity tol is computed, as max(rtol[i])
   if this is positive, or max(atol[i]/fabs(y[i])) otherwise, adjusted
   so as to be between 100*ETA and 0.001.
   Then the computed value h0 is given by

      h0^(-2) = 1. / ( tol * w0^2 ) + tol * ( norm(f) )^2

   where   w0     = max( fabs(*t), fabs(tout) ),
	   f      = the initial value of the vector f(t,y), and
	   norm() = the weighted vector norm used throughout, given by
		    the vmnorm function routine, and weighted by the
		    tolerances initially loaded into the ewt array.

   The sign of h0 is inferred from the initial values of tout and *t.
   fabs(h0) is made < fabs(tout-*t) in any case.
*/
      if ( h0 == 0. ) {
	 tdist = fabs( tout - *t );
	 w0 = max( fabs( *t ), fabs( tout ) );
	 if ( tdist < 2. * ETA * w0 ) {
	    if ( prfl )
	       printf( "\nlsoda -- tout too close to t to start integration\n ");
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 tol = rtol[1];
	 if ( itol > 2 ) {
	    for ( i = 2 ; i <= n ; i++ )
	       tol = max( tol, rtol[i] );
	 }
	 if ( tol <= 0. ) {
	    atoli = atol[1];
	    for ( i = 1 ; i <= n ; i++ ) {
	       if ( itol == 2 || itol == 4 )
		  atoli = atol[i];
	       ayi = fabs( y[i] );
	       if ( ayi != 0. )
		  tol = max( tol, atoli / ayi );
	    }
	 }
	 tol = max( tol, 100. * ETA );
	 tol = min( tol, 0.001 );
	 sum = vmnorm( n, yh[2], ewt );
	 sum = 1. / ( tol * w0 * w0 ) + tol * sum * sum;
	 h0 = 1. / sqrt( sum );
	 h0 = min( h0, tdist );
	 h0 = h0 * ( ( tout - *t >= 0. ) ? 1. : -1. );
      }                 /*   end if ( h0 == 0. )   */
/*
   Adjust h0 if necessary to meet hmax bound.
*/
      rh = fabs( h0 ) * hmxi;
      if ( rh > 1. )
	 h0 /= rh;
/*
   Load h with h0 and scale yh[2] by h0.
*/
      h = h0;
      yp1 = yh[2];
      for ( i = 1 ; i <= n ; i++ )
	 yp1[i] *= h0;
   }         /* if ( *istate == 1 )   */

/*
   Block d.
   The next code block is for continuation calls only ( *istate = 2 or 3 )
   and is to check stop conditions before taking a step.
*/
   if ( *istate == 2 || *istate == 3 ) {
      nslast = nst;
      switch ( itask ) {
      case 1 :
	 if ( ( tn - tout ) * h >= 0. ) {
	    intdy( tout, 0, y, &iflag );
	    if ( iflag != 0 ) {
	       if ( prfl )
		  printf( "\nlsoda -- trouble from intdy, itask = %d, tout = %g",
			   itask, tout );
	       terminate( istate );
	       freevectors();
	       return;
	    }
	    *t = tout;
	    *istate = 2;
	    illin = 0;
	    freevectors();
	    return;
	 }
	 break;
      case 2 :
	 break;
      case 3 :
	 tp = tn - hu * ( 1. + 100. * ETA );
	 if ( ( tp - tout ) * h > 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- itask = %d and tout behind tcur - hu", itask );
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 if ( ( tn - tout ) * h < 0. )
	    break;
	 successreturn( y, t, itask, ihit, tcrit, istate );
	 return;
      case 4 :
	 tcrit = rwork1;
	 if ( ( tn - tcrit ) * h > 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- itask = 4 or 5 and tcrit behind tcur" );
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 if ( ( tcrit - tout ) * h < 0. ) {
	    if ( prfl )
	       printf( "\nlsoda -- itask = 4 or 5 and tcrit behind tout" );
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 if ( ( tn - tout ) * h >= 0. ) {
	    intdy( tout, 0, y, &iflag );
	    if ( iflag != 0 ) {
	       if ( prfl )
		  printf( "\nlsoda -- trouble from intdy, itask = %d, tout = %g",
			   itask, tout );
	       terminate( istate );
	       freevectors();
	       return;
	    }
	    *t = tout;
	    *istate = 2;
	    illin = 0;
	    freevectors();
	    return;
	 }
      case 5 :
	 if ( itask == 5 ) {
	    tcrit = rwork1;
	    if ( ( tn - tcrit ) * h > 0. ) {
	       if ( prfl )
		  printf( "\nlsoda -- itask = 4 or 5 and tcrit behind tcur" );
	       terminate( istate );
	       freevectors();
	       return;
	    }
	 }
	 hmx = fabs( tn ) + fabs( h );
	 ihit = fabs( tn - tcrit ) <= ( 100. * ETA * hmx );
	 if ( ihit ) {
	    *t = tcrit;
	    successreturn( y, t, itask, ihit, tcrit, istate );
	    return;
	 }
	 tnext = tn + h * ( 1. + 4. * ETA );
	 if ( ( tnext - tcrit ) * h <= 0. )
	    break;
	 h = ( tcrit - tn ) * ( 1. - 4. * ETA );
	 if ( *istate == 2 )
	    jstart = -2;
	 break;
      }      /*   end switch   */
   }      /*   end if ( *istate == 2 || *istate == 3 )   */

/*
   Block e.
   The next block is normally executed for all calls and contains
   the call to the one-step core integrator stoda.

   This is a looping point for the integration steps.

   First check for too many steps being taken, update ewt ( if not at
   start of problem).  Check for too much accuracy being requested, and
   check for h below the roundoff level in *t.
*/
   while ( 1 ) {
      if ( *istate != 1 || nst != 0 ) {
	 if ( ( nst - nslast ) >= mxstep ) {
	    if ( prfl )
	       printf( "\nlsoda -- %d steps taken before reaching tout", mxstep );
	    *istate = -1;
	    terminate2( y, t );
	    return;
	 }
	 ewset( itol, rtol, atol, yh[1] );
	 for ( i = 1 ; i <= n ; i++ ) {
	    if ( ewt[i] <= 0. ) {
	       if ( prfl )
		  printf( "\nlsoda -- ewt[%d] = %g <= 0.", i, ewt[i] );
	       *istate = -6;
	       terminate2( y, t );
	       return;
	    }
	    ewt[i] = 1. / ewt[i];
	 }
      }
      tolsf = ETA * vmnorm( n, yh[1], ewt );
      if ( tolsf > 0.01 ) {
	 tolsf = tolsf * 200.;
	 if ( nst == 0 ) {
	    if ( prfl ) {
	       printf( "\nlsoda -- at start of problem, too much accuracy" );
	       printf( "\n         requested for precision of machine," );
	       printf( "\n         suggested scaling factor = %g", tolsf );
	    }
	    terminate( istate );
	    freevectors();
	    return;
	 }
	 if ( prfl ) {
	    printf( "\nlsoda -- at t = %g, too much accuracy requested", *t );
	    printf( "\n         for precision of machine, suggested" );
	    printf( "\n         scaling factor = %g", tolsf );
	 }
	 *istate = -2;
	 terminate2( y, t );
	 return;
      }
      if ( ( tn + h ) == tn ) {
	 nhnil++;
	 if ( nhnil <= mxhnil ) {
	    if ( prfl ) {
	       printf( "\nlsoda -- warning..internal t = %g and h = %g are", tn, h );
	       printf( "\n         such that in the machine, t + h = t on the next step" );
	       printf( "\n         solver will continue anyway." );
	    }
	    if ( ( nhnil == mxhnil ) && prfl ) {
	       printf( "\nlsoda -- above warning has been issued %d times,",
		  nhnil );
	       printf( "\n         it will not be issued again for this problem" );
	    }
	 }
      }

/*
   Call stoda
*/
      stoda( neq, y, f );

/*
   Print extra information
*/
      if ( ( ixpr == 2 ) && prfl )
      {
       printf( "\nmeth= %d,   order= %d,   nfe= %d,   nje= %d",