1 2 SUBROUTINE DDRIV3 (N,T,Y,F,NSTATE,TOUT,NTASK,NROOT,EPS,EWT,IERROR, 3 8 MINT,MITER,IMPL,ML,MU,MXORD,HMAX,WORK,LENW,IWORK,LENIW,JACOBN, 4 8 FA,NDE,MXSTEP,G) 5C***BEGIN PROLOGUE DDRIV3 6C***DATE WRITTEN 790601 (YYMMDD) 7C***REVISION DATE 850924 (YYMMDD) 8C***CATEGORY NO. I1A2,I1A1B 9C***KEYWORDS ODE,STIFF,ORDINARY DIFFERENTIAL EQUATIONS, 10C INITIAL VALUE PROBLEMS,GEAR'S METHOD, 11C DOUBLE PRECISION 12C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS, 13C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY 14C***PURPOSE The function of DDRIV3 is to solve N ordinary differential 15C equations of the form dY(I)/dT = F(Y(I),T), given the 16C initial conditions Y(I) = YI. The program has options to 17C allow the solution of both stiff and non-stiff differential 18C equations. Other important options are available. DDRIV3 19C uses double precision arithmetic. 20C***DESCRIPTION 21C 22C I. ABSTRACT ....................................................... 23C 24C The primary function of DDRIV3 is to solve N ordinary differential 25C equations of the form dY(I)/dT = F(Y(I),T), given the initial 26C conditions Y(I) = YI. The program has options to allow the 27C solution of both stiff and non-stiff differential equations. In 28C addition, DDRIV3 may be used to solve: 29C 1. The initial value problem, A*dY(I)/dT = F(Y(I),T), where A is 30C a non-singular matrix depending on Y and T. 31C 2. The hybrid differential/algebraic initial value problem, 32C A*dY(I)/dT = F(Y(I),T), where A is a vector (whose values may 33C depend upon Y and T) some of whose components will be zero 34C corresponding to those equations which are algebraic rather 35C than differential. 36C DDRIV3 is to be called once for each output point of T. 37C 38C II. PARAMETERS .................................................... 39C 40C (REMEMBER--To run DDRIV3 correctly in double precision, ALL 41C non-integer arguments in the call sequence, including 42C arrays, MUST be declared double precision.) 43C The user should use parameter names in the call sequence of DDRIV3 44C for those quantities whose value may be altered by DDRIV3. The 45C parameters in the call sequence are: 46C 47C N = (Input) The number of dependent functions whose solution 48C is desired. N must not be altered during a problem. 49C 50C T = The independent variable. On input for the first call, T 51C is the initial point. On output, T is the point at which 52C the solution is given. 53C 54C Y = The vector of dependent variables. Y is used as input on 55C the first call, to set the initial values. On output, Y 56C is the computed solution vector. This array Y is passed 57C in the call sequence of the user-provided routines F, 58C JACOBN, FA, USERS, and G. 59C 60C F = A subroutine supplied by the user. The name must be 61C declared EXTERNAL in the user's calling program. This 62C subroutine is of the form: 63C SUBROUTINE F (N, T, Y, YDOT) 64C DOUBLE PRECISION Y(*), YDOT(*) 65C . 66C . 67C YDOT(1) = ... 68C . 69C . 70C YDOT(N) = ... 71C END (Sample) 72C This computes YDOT = F(Y,T), the right hand side of the 73C differential equations. Here Y is a vector of length at 74C least N. The actual length of Y is determined by the 75C user's declaration in the program which calls DDRIV3. 76C Thus the dimensioning of Y in F, while required by FORTRAN 77C convention, does not actually allocate any storage. When 78C this subroutine is called, the first N components of Y are 79C intermediate approximations to the solution components. 80C The user should not alter these values. Here YDOT is a 81C vector of length N. The user should only compute YDOT(I) 82C for I from 1 to N. 83C 84C NSTATE = An integer describing the status of integration. The 85C meaning of NSTATE is as follows: 86C 1 (Input) Means the first call to the routine. This 87C value must be set by the user. On all subsequent 88C calls the value of NSTATE should be tested by the 89C user, but must not be altered. (As a convenience to 90C the user who may wish to put out the initial 91C conditions, DDRIV3 can be called with NSTATE=1, and 92C TOUT=T. In this case the program will return with 93C NSTATE unchanged, i.e., NSTATE=1.) 94C 2 (Output) Means a successful integration. If a normal 95C continuation is desired (i.e., a further integration 96C in the same direction), simply advance TOUT and call 97C again. All other parameters are automatically set. 98C 3 (Output)(Unsuccessful) Means the integrator has taken 99C MXSTEP steps without reaching TOUT. The user can 100C continue the integration by simply calling DDRIV3 101C again. 102C 4 (Output)(Unsuccessful) Means too much accuracy has 103C been requested. EPS has been increased to a value 104C the program estimates is appropriate. The user can 105C continue the integration by simply calling DDRIV3 106C again. 107C 5 (Output) A root was found at a point less than TOUT. 108C The user can continue the integration toward TOUT by 109C simply calling DDRIV3 again. 110C 111C TOUT = (Input) The point at which the solution is desired. The 112C position of TOUT relative to T on the first call 113C determines the direction of integration. 114C 115C NTASK = (Input) An index specifying the manner of returning the 116C solution, according to the following: 117C NTASK = 1 Means DDRIV3 will integrate past TOUT and 118C interpolate the solution. This is the most 119C efficient mode. 120C NTASK = 2 Means DDRIV3 will return the solution after 121C each internal integration step, or at TOUT, 122C whichever comes first. In the latter case, 123C the program integrates exactly to TOUT. 124C NTASK = 3 Means DDRIV3 will adjust its internal step to 125C reach TOUT exactly (useful if a singularity 126C exists beyond TOUT.) 127C 128C NROOT = (Input) The number of equations whose roots are desired. 129C If NROOT is zero, the root search is not active. This 130C option is useful for obtaining output at points which are 131C not known in advance, but depend upon the solution, e.g., 132C when some solution component takes on a specified value. 133C The root search is carried out using the user-written 134C function G (see description of G below.) DDRIV3 attempts 135C to find the value of T at which one of the equations 136C changes sign. DDRIV3 can find at most one root per 137C equation per internal integration step, and will then 138C return the solution either at TOUT or at a root, whichever 139C occurs first in the direction of integration. The index 140C of the equation whose root is being reported is stored in 141C the sixth element of IWORK. 142C NOTE: NROOT is never altered by this program. 143C 144C EPS = On input, the requested relative accuracy in all solution 145C components. EPS = 0 is allowed. On output, the adjusted 146C relative accuracy if the input value was too small. The 147C value of EPS should be set as large as is reasonable, 148C because the amount of work done by DDRIV3 increases as EPS 149C decreases. 150C 151C EWT = (Input) Problem zero, i.e., the smallest, nonzero, 152C physically meaningful value for the solution. (Array, 153C possibly of length one. See following description of 154C IERROR.) Setting EWT smaller than necessary can adversely 155C affect the running time. 156C 157C IERROR = (Input) Error control indicator. A value of 3 is 158C suggested for most problems. Other choices and detailed 159C explanations of EWT and IERROR are given below for those 160C who may need extra flexibility. 161C 162C These last three input quantities EPS, EWT and IERROR 163C control the accuracy of the computed solution. EWT and 164C IERROR are used internally to compute an array YWT. One 165C step error estimates divided by YWT(I) are kept less than 166C EPS in root mean square norm. 167C IERROR (Set by the user) = 168C 1 Means YWT(I) = 1. (Absolute error control) 169C EWT is ignored. 170C 2 Means YWT(I) = ABS(Y(I)), (Relative error control) 171C EWT is ignored. 172C 3 Means YWT(I) = MAX(ABS(Y(I)), EWT(1)). 173C 4 Means YWT(I) = MAX(ABS(Y(I)), EWT(I)). 174C This choice is useful when the solution components 175C have differing scales. 176C 5 Means YWT(I) = EWT(I). 177C If IERROR is 3, EWT need only be dimensioned one. 178C If IERROR is 4 or 5, the user must dimension EWT at least 179C N, and set its values. 180C 181C MINT = (Input) The integration method indicator. 182C MINT = 1 Means the Adams methods, and is used for 183C non-stiff problems. 184C MINT = 2 Means the stiff methods of Gear (i.e., the 185C backward differentiation formulas), and is 186C used for stiff problems. 187C MINT = 3 Means the program dynamically selects the 188C Adams methods when the problem is non-stiff 189C and the Gear methods when the problem is 190C stiff. When using the Adams methods, the 191C program uses a value of MITER=0; when using 192C the Gear methods, the program uses the value 193C of MITER provided by the user. Only a value 194C of IMPL = 0 and a value of MITER = 1, 2, 4, or 195C 5 is allowed for this option. The user may 196C not alter the value of MINT or MITER without 197C restarting, i.e., setting NSTATE to 1. 198C 199C MITER = (Input) The iteration method indicator. 200C MITER = 0 Means functional iteration. This value is 201C suggested for non-stiff problems. 202C MITER = 1 Means chord method with analytic Jacobian. 203C In this case, the user supplies subroutine 204C JACOBN (see description below). 205C MITER = 2 Means chord method with Jacobian calculated 206C internally by finite differences. 207C MITER = 3 Means chord method with corrections computed 208C by the user-written routine named USERS. 209C This option allows all matrix algebra and 210C storage decisions to be made by the user. 211C The routine USERS is called by DDRIV3 when 212C certain linear systems must be solved. The 213C user may choose any method to form, store and 214C solve these systems in order to obtain the 215C solution result that is returned to DDRIV3. 216C In particular, this allows sparse matrix 217C methods to be used. 218C The call sequence for this routine is 219C 220C SUBROUTINE USERS (Y, YH, YWT, SAVE1, SAVE2, 221C 8 T, H, EL, IMPL, N, NDE, IFLAG) 222C DOUBLE PRECISION Y(*), YH(*), YWT(*), 223C 8 SAVE1(*), SAVE2(*), T, H, EL 224C 225C The input variable IFLAG indicates what 226C action is to be taken. Subroutine USERS 227C should perform the following operations, 228C depending on the value of IFLAG and IMPL. 229C 230C IFLAG = 0 231C IMPL = 0. USERS is not called. 232C IMPL = 1 or 2. Solve the system 233C A*X = SAVE2, 234C returning the result in SAVE2. The array 235C SAVE1 can be used as a work array. 236C 237C IFLAG = 1 238C IMPL = 0. Compute, decompose and store the 239C matrix (I - H*EL*J), where I is the 240C identity matrix and J is the Jacobian 241C matrix of the right hand side. The array 242C SAVE1 can be used as a work array. 243C IMPL = 1 or 2. Compute, decompose and store 244C the matrix (A - H*EL*J). The array SAVE1 245C can be used as a work array. 246C 247C IFLAG = 2 248C IMPL = 0. Solve the system 249C (I - H*EL*J)*X = H*SAVE2 - YH - SAVE1, 250C returning the result in SAVE2. 251C IMPL = 1, or 2. Solve the system 252C (A - H*EL*J)*X = H*SAVE2 - A*(YH + SAVE1) 253C returning the result in SAVE2. 254C The array SAVE1 should not be altered. 255C 256C When using a value of MITER = 3, the 257C subroutine FA is not required, even if IMPL 258C is not 0. For further information on using 259C this option, see section IV-F below. 260C 261C MITER = 4 Means the same as MITER = 1 but the A and 262C Jacobian matrices are assumed to be banded. 263C MITER = 5 Means the same as MITER = 2 but the A and 264C Jacobian matrices are assumed to be banded. 265C 266C IMPL = (Input) The implicit method indicator. 267C IMPL = 0 Means solving dY(I)/dT = F(Y(I),T). 268C IMPL = 1 Means solving A*dY(I)/dT = F(Y(I),T), 269C non-singular A (see description of FA below.) 270C Only MINT = 1 or 2, and MITER = 1, 2, 3, 4, or 271C 5 are allowed for this option. 272C IMPL = 2 Means solving certain systems of hybrid 273C differential/algebraic equations (see 274C description of FA below.) Only MINT = 2 and 275C MITER = 1, 2, 3, 4, or 5, are allowed for this 276C option. 277C The value of IMPL must not be changed during a problem. 278C 279C ML = (Input) The lower half-bandwidth in the case of a banded 280C A or Jacobian matrix. (I.e., maximum(R-C) for nonzero 281C A(R,C).) 282C 283C MU = (Input) The upper half-bandwidth in the case of a banded 284C A or Jacobian matrix. (I.e., maximum(C-R).) 285C 286C MXORD = (Input) The maximum order desired. This is .LE. 12 for 287C the Adams methods and .LE. 5 for the Gear methods. Normal 288C value is 12 and 5, respectively. If MINT is 3, the 289C maximum order used will be MIN(MXORD, 12) when using the 290C Adams methods, and MIN(MXORD, 5) when using the Gear 291C methods. MXORD must not be altered during a problem. 292C 293C HMAX = (Input) The maximum magnitude of the step size that will 294C be used for the problem. This is useful for ensuring that 295C important details are not missed. If this is not the 296C case, a large value, such as the interval length, is 297C suggested. 298C 299C WORK 300C LENW = (Input) 301C WORK is an array of LENW double precision words used 302C internally for temporary storage. The user must allocate 303C space for this array in the calling program by a statement 304C such as 305C DOUBLE PRECISION WORK(...) 306C The following table gives the required minimum value for 307C the length of WORK, depending on the value of IMPL and 308C MITER. LENW should be set to the value used. The 309C contents of WORK should not be disturbed between calls to 310C DDRIV3. 311C 312C IMPL = 0 1 2 313C --------------------------------------------------------- 314C MITER = 0 (MXORD+4)*N + Not allowed Not allowed 315C 2*NROOT + 204 316C 317C 1,2 N*N+(MXORD+4)*N 2*N*N+(MXORD+4)*N N*N+(MXORD+5)*N 318C + 2*NROOT + 204 + 2*NROOT + 204 + 2*NROOT + 204 319C 320C 3 (MXORD+4)*N + (MXORD+4)*N + (MXORD+4)*N + 321C 2*NROOT + 204 2*NROOT + 204 2*NROOT + 204 322C 323C 4,5 (2*ML+MU)*N + (4*ML+2*MU)*N + (2*ML+MU)*N + 324C (MXORD+5)*N + (MXORD+6)*N + (MXORD+6)*N + 325C 2*NROOT + 204 2*NROOT + 204 2*NROOT + 204 326C --------------------------------------------------------- 327C 328C IWORK 329C LENIW = (Input) 330C IWORK is an integer array of length LENIW used internally 331C for temporary storage. The user must allocate space for 332C this array in the calling program by a statement such as 333C INTEGER IWORK(...) 334C The length of IWORK should be at least 335C 21 if MITER is 0 or 3, or 336C N+21 if MITER is 1, 2, 4, or 5, or MINT is 3, 337C and LENIW should be set to the value used. The contents 338C of IWORK should not be disturbed between calls to DDRIV3. 339C 340C JACOBN = A subroutine supplied by the user, if MITER is 1 or 4. 341C If this is the case, the name must be declared EXTERNAL in 342C the user's calling program. Given a system of N 343C differential equations, it is meaningful to speak about 344C the partial derivative of the I-th right hand side with 345C respect to the J-th dependent variable. In general there 346C are N*N such quantities. Often however the equations can 347C be ordered so that the I-th differential equation only 348C involves dependent variables with index near I, e.g., I+1, 349C I-2. Such a system is called banded. If, for all I, the 350C I-th equation depends on at most the variables 351C Y(I-ML), Y(I-ML+1), ... , Y(I), Y(I+1), ... , Y(I+MU) 352C then we call ML+MU+1 the bandwith of the system. In a 353C banded system many of the partial derivatives above are 354C automatically zero. For the cases MITER = 1, 2, 4, and 5, 355C some of these partials are needed. For the cases 356C MITER = 2 and 5 the necessary derivatives are 357C approximated numerically by DDRIV3, and we only ask the 358C user to tell DDRIV3 the value of ML and MU if the system 359C is banded. For the cases MITER = 1 and 4 the user must 360C derive these partials algebraically and encode them in 361C subroutine JACOBN. By computing these derivatives the 362C user can often save 20-30 per cent of the computing time. 363C Usually, however, the accuracy is not much affected and 364C most users will probably forego this option. The optional 365C user-written subroutine JACOBN has the form: 366C SUBROUTINE JACOBN (N, T, Y, DFDY, MATDIM, ML, MU) 367C DOUBLE PRECISION Y(*), DFDY(MATDIM,*) 368C . 369C . 370C Calculate values of DFDY 371C . 372C . 373C END (Sample) 374C Here Y is a vector of length at least N. The actual 375C length of Y is determined by the user's declaration in the 376C program which calls DDRIV3. Thus the dimensioning of Y in 377C JACOBN, while required by FORTRAN convention, does not 378C actually allocate any storage. When this subroutine is 379C called, the first N components of Y are intermediate 380C approximations to the solution components. The user 381C should not alter these values. If the system is not 382C banded (MITER=1), the partials of the I-th equation with 383C respect to the J-th dependent function are to be stored in 384C DFDY(I,J). Thus partials of the I-th equation are stored 385C in the I-th row of DFDY. If the system is banded 386C (MITER=4), then the partials of the I-th equation with 387C respect to Y(J) are to be stored in DFDY(K,J), where 388C K=I-J+MU+1. 389C 390C FA = A subroutine supplied by the user if IMPL is 1 or 2, and 391C MITER is not 3. If so, the name must be declared EXTERNAL 392C in the user's calling program. This subroutine computes 393C the array A, where A*dY(I)/dT = F(Y(I),T). 394C There are two cases: 395C 396C IMPL=1. 397C Subroutine FA is of the form: 398C SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE) 399C DOUBLE PRECISION Y(*), A(MATDIM,*) 400C . 401C . 402C Calculate ALL values of A 403C . 404C . 405C END (Sample) 406C In this case A is assumed to be a nonsingular matrix, 407C with the same structure as DFDY (see JACOBN description 408C above). Programming considerations prevent complete 409C generality. If MITER is 1 or 2, A is assumed to be full 410C and the user must compute and store all values of 411C A(I,J), I,J=1, ... ,N. If MITER is 4 or 5, A is assumed 412C to be banded with lower and upper half bandwidth ML and 413C MU. The left hand side of the I-th equation is a linear 414C combination of dY(I-ML)/dT, dY(I-ML+1)/dT, ... , 415C dY(I)/dT, ... , dY(I+MU-1)/dT, dY(I+MU)/dT. Thus in the 416C I-th equation, the coefficient of dY(J)/dT is to be 417C stored in A(K,J), where K=I-J+MU+1. 418C NOTE: The array A will be altered between calls to FA. 419C 420C IMPL=2. 421C Subroutine FA is of the form: 422C SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE) 423C DOUBLE PRECISION Y(*), A(*) 424C . 425C . 426C Calculate non-zero values of A(1),...,A(NDE) 427C . 428C . 429C END (Sample) 430C In this case it is assumed that the system is ordered by 431C the user so that the differential equations appear 432C first, and the algebraic equations appear last. The 433C algebraic equations must be written in the form: 434C 0 = F(Y(I),T). When using this option it is up to the 435C user to provide initial values for the Y(I) that satisfy 436C the algebraic equations as well as possible. It is 437C further assumed that A is a vector of length NDE. All 438C of the components of A, which may depend on T, Y(I), 439C etc., must be set by the user to non-zero values. 440C Here Y is a vector of length at least N. The actual 441C length of Y is determined by the user's declaration in the 442C program which calls DDRIV3. Thus the dimensioning of Y in 443C FA, while required by FORTRAN convention, does not 444C actually allocate any storage. When this subroutine is 445C called, the first N components of Y are intermediate 446C approximations to the solution components. The user 447C should not alter these values. FA is always called 448C immediately after calling F, with the same values of T 449C and Y. 450C 451C NDE = (Input) The number of differential equations. This is 452C required only for IMPL = 2, with NDE .LT. N. 453C 454C MXSTEP = (Input) The maximum number of internal steps allowed on 455C one call to DDRIV3. 456C 457C G = A double precision FORTRAN function supplied by the user 458C if NROOT is not 0. In this case, the name must be 459C declared EXTERNAL in the user's calling program. G is 460C repeatedly called with different values of IROOT to obtain 461C the value of each of the NROOT equations for which a root 462C is desired. G is of the form: 463C DOUBLE PRECISION FUNCTION G (N, T, Y, IROOT) 464C DOUBLE PRECISION Y(*) 465C GO TO (10, ...), IROOT 466C 10 G = ... 467C . 468C . 469C END (Sample) 470C Here, Y is a vector of length at least N, whose first N 471C components are the solution components at the point T. 472C The user should not alter these values. The actual length 473C of Y is determined by the user's declaration in the 474C program which calls DDRIV3. Thus the dimensioning of Y in 475C G, while required by FORTRAN convention, does not actually 476C allocate any storage. 477C 478C***LONG DESCRIPTION 479C 480C III. OTHER COMMUNICATION TO THE USER .............................. 481C 482C A. The solver communicates to the user through the parameters 483C above. In addition it writes diagnostic messages through the 484C standard error handling program XERROR. That program will 485C terminate the user's run if it detects a probable problem setup 486C error, e.g., insufficient storage allocated by the user for the 487C WORK array. Messages are written on the standard error message 488C file. At installations which have this error handling package 489C the user should determine the standard error handling file from 490C the local documentation. Otherwise the short but serviceable 491C routine, XERROR, available with this package, can be used. That 492C program writes on logical unit 6 to transmit messages. A 493C complete description of XERROR is given in the Sandia 494C Laboratories report SAND78-1189 by R. E. Jones. Following is a 495C list of possible errors. Unless otherwise noted, all messages 496C come from DDRIV3: 497C 498C No. Type Message 499C --- ---- ------- 500C 1 Fatal From DDRIV2: The integration method flag has 501C an illegal value. 502C 2 Warning The output point is inconsistent with the 503C value of NTASK and T. 504C 3 Warning Number of steps to reach TOUT exceeds MXSTEP. 505C 4 Recoverable Requested accuracy is too stringent. 506C 5 Warning Step size is below the roundoff level. 507C 6 Fatal EPS is less than zero. 508C 7 Fatal N is not positive. 509C 8 Fatal Insufficient work space provided. 510C 9 Fatal Improper value for MINT, MITER and/or IMPL. 511C 10 Fatal The IWORK array is too small. 512C 11 Fatal The step size has gone to zero. 513C 12 Fatal Excessive amount of work. 514C 13 Fatal For IMPL=1 or 2, the matrix A is singular. 515C 14 Fatal MXORD is not positive. 516C 15 Fatal From DDRIV1: N is greater than 200. 517C 16 Fatal From DDRIV1: The WORK array is too small. 518C 519C B. The first three elements of WORK and the first five elements of 520C IWORK will contain the following statistical data: 521C AVGH The average step size used. 522C HUSED The step size last used (successfully). 523C AVGORD The average order used. 524C IMXERR The index of the element of the solution vector that 525C contributed most to the last error test. 526C NQUSED The order last used (successfully). 527C NSTEP The number of steps taken. 528C NFE The number of evaluations of the right hand side. 529C NJE The number of evaluations of the Jacobian matrix. 530C 531C IV. REMARKS ....................................................... 532C 533C A. Other routines used: 534C DDNTP, DDZRO, DDSTP, DDNTL, DDPST, DDCOR, DDCST, 535C DDPSC, and DDSCL; 536C DGEFA, DGESL, DGBFA, DGBSL, and DNRM2 (from LINPACK) 537C D1MACH (from the Bell Laboratories Machine Constants Package) 538C XERROR (from the SLATEC Common Math Library) 539C The last seven routines above, not having been written by the 540C present authors, are not explicitly part of this package. 541C 542C B. On any return from DDRIV3 all information necessary to continue 543C the calculation is contained in the call sequence parameters, 544C including the work arrays. Thus it is possible to suspend one 545C problem, integrate another, and then return to the first. 546C 547C C. There are user-written routines which are only required by 548C DDRIV3 when certain parameters are set. Thus a message warning 549C of unsatisfied externals may be issued during the load or link 550C phase. This message should never refer to F. This message can 551C be ignored if: it refers to JACOBN and MITER is not 1 or 4, or 552C it refers to FA and IMPL is 0 or MITER is 3, or it refers to 553C USERS and MITER is not 3, or it refers to G and NROOT is 0. 554C 555C D. If this package is to be used in an overlay situation, the user 556C must declare in the primary overlay the variables in the call 557C sequence to DDRIV3. 558C 559C E. Changing parameters during an integration. 560C The value of NROOT, EPS, EWT, IERROR, MINT, MITER, or HMAX may 561C be altered by the user between calls to DDRIV3. For example, if 562C too much accuracy has been requested (the program returns with 563C NSTATE = 4 and an increased value of EPS) the user may wish to 564C increase EPS further. In general, prudence is necessary when 565C making changes in parameters since such changes are not 566C implemented until the next integration step, which is not 567C necessarily the next call to DDRIV3. This can happen if the 568C program has already integrated to a point which is beyond the 569C new point TOUT. 570C 571C F. As the price for complete control of matrix algebra, the DDRIV3 572C USERS option puts all responsibility for Jacobian matrix 573C evaluation on the user. It is often useful to approximate 574C numerically all or part of the Jacobian matrix. However this 575C must be done carefully. The FORTRAN sequence below illustrates 576C the method we recommend. It can be inserted directly into 577C subroutine USERS to approximate Jacobian elements in rows I1 578C to I2 and columns J1 to J2. 579C DOUBLE PRECISION DFDY(N,N), EPSJ, H, R, D1MACH, 580C 8 SAVE1(N), SAVE2(N), T, UROUND, Y(N), YJ, YWT(N) 581C UROUND = D1MACH(4) 582C EPSJ = UROUND**(1.D0/3.D0) 583C DO 30 J = J1,J2 584C R = EPSJ*MAX(ABS(YWT(J)), ABS(Y(J))) 585C IF (R .EQ. 0.D0) R = EPSJ 586C YJ = Y(J) 587C Y(J) = Y(J) + R 588C CALL F (N, T, Y, SAVE1) 589C Y(J) = YJ 590C DO 20 I = I1,I2 591C 20 DFDY(I,J) = (SAVE1(I) - SAVE2(I))/R 592C 30 CONTINUE 593C Many problems give rise to structured sparse Jacobians, e.g., 594C block banded. It is possible to approximate them with fewer 595C function evaluations than the above procedure uses; see Curtis, 596C Powell and Reid, J. Inst. Maths Applics, (1974), Vol. 13, 597C pp. 117-119. 598C 599C***REFERENCES GEAR, C. W., "NUMERICAL INITIAL VALUE PROBLEMS IN 600C ORDINARY DIFFERENTIAL EQUATIONS", PRENTICE-HALL, 1971. 601C***ROUTINES CALLED DDSTP,SDNTP,SDZRO,DGEFA,DGESL,DGBFA,DGBSL,DNRM2, 602C D1MACH,XERROR 603C***END PROLOGUE DDRIV3 604 605 606