1*DECK CDRIV3 2 SUBROUTINE CDRIV3 (N, T, Y, F, NSTATE, TOUT, NTASK, NROOT, EPS, 3 8 EWT, IERROR, MINT, MITER, IMPL, ML, MU, MXORD, HMAX, WORK, 4 8 LENW, IWORK, LENIW, JACOBN, FA, NDE, MXSTEP, G, USERS, IERFLG) 5C***BEGIN PROLOGUE CDRIV3 6C***PURPOSE The function of CDRIV3 is to solve N ordinary differential 7C equations of the form dY(I)/dT = F(Y(I),T), given the 8C initial conditions Y(I) = YI. The program has options to 9C allow the solution of both stiff and non-stiff differential 10C equations. Other important options are available. CDRIV3 11C allows complex-valued differential equations. 12C***LIBRARY SLATEC (SDRIVE) 13C***CATEGORY I1A2, I1A1B 14C***TYPE COMPLEX (SDRIV3-S, DDRIV3-D, CDRIV3-C) 15C***KEYWORDS COMPLEX VALUED, GEAR'S METHOD, INITIAL VALUE PROBLEMS, 16C ODE, ORDINARY DIFFERENTIAL EQUATIONS, SDRIVE, STIFF 17C***AUTHOR Kahaner, D. K., (NIST) 18C National Institute of Standards and Technology 19C Gaithersburg, MD 20899 20C Sutherland, C. D., (LANL) 21C Mail Stop D466 22C Los Alamos National Laboratory 23C Los Alamos, NM 87545 24C***DESCRIPTION 25C 26C I. ABSTRACT ....................................................... 27C 28C The primary function of CDRIV3 is to solve N ordinary differential 29C equations of the form dY(I)/dT = F(Y(I),T), given the initial 30C conditions Y(I) = YI. The program has options to allow the 31C solution of both stiff and non-stiff differential equations. In 32C addition, CDRIV3 may be used to solve: 33C 1. The initial value problem, A*dY(I)/dT = F(Y(I),T), where A is 34C a non-singular matrix depending on Y and T. 35C 2. The hybrid differential/algebraic initial value problem, 36C A*dY(I)/dT = F(Y(I),T), where A is a vector (whose values may 37C depend upon Y and T) some of whose components will be zero 38C corresponding to those equations which are algebraic rather 39C than differential. 40C CDRIV3 is to be called once for each output point of T. 41C 42C II. PARAMETERS .................................................... 43C 44C The user should use parameter names in the call sequence of CDRIV3 45C for those quantities whose value may be altered by CDRIV3. The 46C parameters in the call sequence are: 47C 48C N = (Input) The number of dependent functions whose solution 49C is desired. N must not be altered during a problem. 50C 51C T = (Real) The independent variable. On input for the first 52C call, T is the initial point. On output, T is the point 53C at which the solution is given. 54C 55C Y = (Complex) The vector of dependent variables. Y is used as 56C input on the first call, to set the initial values. On 57C output, Y is the computed solution vector. This array Y 58C is passed in the call sequence of the user-provided 59C routines F, JACOBN, FA, USERS, and G. Thus parameters 60C required by those routines can be stored in this array in 61C components N+1 and above. (Note: Changes by the user to 62C the first N components of this array will take effect only 63C after a restart, i.e., after setting NSTATE to 1 .) 64C 65C F = A subroutine supplied by the user. The name must be 66C declared EXTERNAL in the user's calling program. This 67C subroutine is of the form: 68C SUBROUTINE F (N, T, Y, YDOT) 69C COMPLEX Y(*), YDOT(*) 70C . 71C . 72C YDOT(1) = ... 73C . 74C . 75C YDOT(N) = ... 76C END (Sample) 77C This computes YDOT = F(Y,T), the right hand side of the 78C differential equations. Here Y is a vector of length at 79C least N. The actual length of Y is determined by the 80C user's declaration in the program which calls CDRIV3. 81C Thus the dimensioning of Y in F, while required by FORTRAN 82C convention, does not actually allocate any storage. When 83C this subroutine is called, the first N components of Y are 84C intermediate approximations to the solution components. 85C The user should not alter these values. Here YDOT is a 86C vector of length N. The user should only compute YDOT(I) 87C for I from 1 to N. Normally a return from F passes 88C control back to CDRIV3. However, if the user would like 89C to abort the calculation, i.e., return control to the 90C program which calls CDRIV3, he should set N to zero. 91C CDRIV3 will signal this by returning a value of NSTATE 92C equal to 6 . Altering the value of N in F has no effect 93C on the value of N in the call sequence of CDRIV3. 94C 95C NSTATE = An integer describing the status of integration. The 96C meaning of NSTATE is as follows: 97C 1 (Input) Means the first call to the routine. This 98C value must be set by the user. On all subsequent 99C calls the value of NSTATE should be tested by the 100C user, but must not be altered. (As a convenience to 101C the user who may wish to put out the initial 102C conditions, CDRIV3 can be called with NSTATE=1, and 103C TOUT=T. In this case the program will return with 104C NSTATE unchanged, i.e., NSTATE=1.) 105C 2 (Output) Means a successful integration. If a normal 106C continuation is desired (i.e., a further integration 107C in the same direction), simply advance TOUT and call 108C again. All other parameters are automatically set. 109C 3 (Output)(Unsuccessful) Means the integrator has taken 110C MXSTEP steps without reaching TOUT. The user can 111C continue the integration by simply calling CDRIV3 112C again. 113C 4 (Output)(Unsuccessful) Means too much accuracy has 114C been requested. EPS has been increased to a value 115C the program estimates is appropriate. The user can 116C continue the integration by simply calling CDRIV3 117C again. 118C 5 (Output) A root was found at a point less than TOUT. 119C The user can continue the integration toward TOUT by 120C simply calling CDRIV3 again. 121C 6 (Output)(Unsuccessful) N has been set to zero in 122C SUBROUTINE F. 123C 7 (Output)(Unsuccessful) N has been set to zero in 124C FUNCTION G. See description of G below. 125C 8 (Output)(Unsuccessful) N has been set to zero in 126C SUBROUTINE JACOBN. See description of JACOBN below. 127C 9 (Output)(Unsuccessful) N has been set to zero in 128C SUBROUTINE FA. See description of FA below. 129C 10 (Output)(Unsuccessful) N has been set to zero in 130C SUBROUTINE USERS. See description of USERS below. 131C 11 (Output)(Successful) For NTASK = 2 or 3, T is beyond 132C TOUT. The solution was obtained by interpolation. 133C The user can continue the integration by simply 134C advancing TOUT and calling CDRIV3 again. 135C 12 (Output)(Unsuccessful) The solution could not be 136C obtained. The value of IERFLG (see description 137C below) for a "Recoverable" situation indicates the 138C type of difficulty encountered: either an illegal 139C value for a parameter or an inability to continue the 140C solution. For this condition the user should take 141C corrective action and reset NSTATE to 1 before 142C calling CDRIV3 again. Otherwise the program will 143C terminate the run. 144C 145C TOUT = (Input, Real) The point at which the solution is desired. 146C The position of TOUT relative to T on the first call 147C determines the direction of integration. 148C 149C NTASK = (Input) An index specifying the manner of returning the 150C solution, according to the following: 151C NTASK = 1 Means CDRIV3 will integrate past TOUT and 152C interpolate the solution. This is the most 153C efficient mode. 154C NTASK = 2 Means CDRIV3 will return the solution after 155C each internal integration step, or at TOUT, 156C whichever comes first. In the latter case, 157C the program integrates exactly to TOUT. 158C NTASK = 3 Means CDRIV3 will adjust its internal step to 159C reach TOUT exactly (useful if a singularity 160C exists beyond TOUT.) 161C 162C NROOT = (Input) The number of equations whose roots are desired. 163C If NROOT is zero, the root search is not active. This 164C option is useful for obtaining output at points which are 165C not known in advance, but depend upon the solution, e.g., 166C when some solution component takes on a specified value. 167C The root search is carried out using the user-written 168C function G (see description of G below.) CDRIV3 attempts 169C to find the value of T at which one of the equations 170C changes sign. CDRIV3 can find at most one root per 171C equation per internal integration step, and will then 172C return the solution either at TOUT or at a root, whichever 173C occurs first in the direction of integration. The initial 174C point is never reported as a root. The index of the 175C equation whose root is being reported is stored in the 176C sixth element of IWORK. 177C NOTE: NROOT is never altered by this program. 178C 179C EPS = (Real) On input, the requested relative accuracy in all 180C solution components. EPS = 0 is allowed. On output, the 181C adjusted relative accuracy if the input value was too 182C small. The value of EPS should be set as large as is 183C reasonable, because the amount of work done by CDRIV3 184C increases as EPS decreases. 185C 186C EWT = (Input, Real) Problem zero, i.e., the smallest, nonzero, 187C physically meaningful value for the solution. (Array, 188C possibly of length one. See following description of 189C IERROR.) Setting EWT smaller than necessary can adversely 190C affect the running time. 191C 192C IERROR = (Input) Error control indicator. A value of 3 is 193C suggested for most problems. Other choices and detailed 194C explanations of EWT and IERROR are given below for those 195C who may need extra flexibility. 196C 197C These last three input quantities EPS, EWT and IERROR 198C control the accuracy of the computed solution. EWT and 199C IERROR are used internally to compute an array YWT. One 200C step error estimates divided by YWT(I) are kept less than 201C EPS in root mean square norm. 202C IERROR (Set by the user) = 203C 1 Means YWT(I) = 1. (Absolute error control) 204C EWT is ignored. 205C 2 Means YWT(I) = ABS(Y(I)), (Relative error control) 206C EWT is ignored. 207C 3 Means YWT(I) = MAX(ABS(Y(I)), EWT(1)). 208C 4 Means YWT(I) = MAX(ABS(Y(I)), EWT(I)). 209C This choice is useful when the solution components 210C have differing scales. 211C 5 Means YWT(I) = EWT(I). 212C If IERROR is 3, EWT need only be dimensioned one. 213C If IERROR is 4 or 5, the user must dimension EWT at least 214C N, and set its values. 215C 216C MINT = (Input) The integration method indicator. 217C MINT = 1 Means the Adams methods, and is used for 218C non-stiff problems. 219C MINT = 2 Means the stiff methods of Gear (i.e., the 220C backward differentiation formulas), and is 221C used for stiff problems. 222C MINT = 3 Means the program dynamically selects the 223C Adams methods when the problem is non-stiff 224C and the Gear methods when the problem is 225C stiff. When using the Adams methods, the 226C program uses a value of MITER=0; when using 227C the Gear methods, the program uses the value 228C of MITER provided by the user. Only a value 229C of IMPL = 0 and a value of MITER = 1, 2, 4, or 230C 5 is allowed for this option. The user may 231C not alter the value of MINT or MITER without 232C restarting, i.e., setting NSTATE to 1. 233C 234C MITER = (Input) The iteration method indicator. 235C MITER = 0 Means functional iteration. This value is 236C suggested for non-stiff problems. 237C MITER = 1 Means chord method with analytic Jacobian. 238C In this case, the user supplies subroutine 239C JACOBN (see description below). 240C MITER = 2 Means chord method with Jacobian calculated 241C internally by finite differences. 242C MITER = 3 Means chord method with corrections computed 243C by the user-written routine USERS (see 244C description of USERS below.) This option 245C allows all matrix algebra and storage 246C decisions to be made by the user. When using 247C a value of MITER = 3, the subroutine FA is 248C not required, even if IMPL is not 0. For 249C further information on using this option, see 250C Section IV-E below. 251C MITER = 4 Means the same as MITER = 1 but the A and 252C Jacobian matrices are assumed to be banded. 253C MITER = 5 Means the same as MITER = 2 but the A and 254C Jacobian matrices are assumed to be banded. 255C 256C IMPL = (Input) The implicit method indicator. 257C IMPL = 0 Means solving dY(I)/dT = F(Y(I),T). 258C IMPL = 1 Means solving A*dY(I)/dT = F(Y(I),T), non- 259C singular A (see description of FA below.) 260C Only MINT = 1 or 2, and MITER = 1, 2, 3, 4, 261C or 5 are allowed for this option. 262C IMPL = 2,3 Means solving certain systems of hybrid 263C differential/algebraic equations (see 264C description of FA below.) Only MINT = 2 and 265C MITER = 1, 2, 3, 4, or 5, are allowed for 266C this option. 267C The value of IMPL must not be changed during a problem. 268C 269C ML = (Input) The lower half-bandwidth in the case of a banded 270C A or Jacobian matrix. (I.e., maximum(R-C) for nonzero 271C A(R,C).) 272C 273C MU = (Input) The upper half-bandwidth in the case of a banded 274C A or Jacobian matrix. (I.e., maximum(C-R).) 275C 276C MXORD = (Input) The maximum order desired. This is .LE. 12 for 277C the Adams methods and .LE. 5 for the Gear methods. Normal 278C value is 12 and 5, respectively. If MINT is 3, the 279C maximum order used will be MIN(MXORD, 12) when using the 280C Adams methods, and MIN(MXORD, 5) when using the Gear 281C methods. MXORD must not be altered during a problem. 282C 283C HMAX = (Input, Real) The maximum magnitude of the step size that 284C will be used for the problem. This is useful for ensuring 285C that important details are not missed. If this is not the 286C case, a large value, such as the interval length, is 287C suggested. 288C 289C WORK 290C LENW = (Input) 291C WORK is an array of LENW complex words used 292C internally for temporary storage. The user must allocate 293C space for this array in the calling program by a statement 294C such as 295C COMPLEX WORK(...) 296C The following table gives the required minimum value for 297C the length of WORK, depending on the value of IMPL and 298C MITER. LENW should be set to the value used. The 299C contents of WORK should not be disturbed between calls to 300C CDRIV3. 301C 302C IMPL = 0 1 2 3 303C --------------------------------------------------------- 304C MITER = 0 (MXORD+4)*N Not allowed Not allowed Not allowed 305C + 2*NROOT 306C + 250 307C 308C 1,2 N*N + 2*N*N + N*N + N*(N + NDE) 309C (MXORD+5)*N (MXORD+5)*N (MXORD+6)*N + (MXORD+5)*N 310C + 2*NROOT + 2*NROOT + 2*NROOT + 2*NROOT 311C + 250 + 250 + 250 + 250 312C 313C 3 (MXORD+4)*N (MXORD+4)*N (MXORD+4)*N (MXORD+4)*N 314C + 2*NROOT + 2*NROOT + 2*NROOT + 2*NROOT 315C + 250 + 250 + 250 + 250 316C 317C 4,5 (2*ML+MU+1) 2*(2*ML+MU+1) (2*ML+MU+1) (2*ML+MU+1)* 318C *N + *N + *N + (N+NDE) + 319C (MXORD+5)*N (MXORD+5)*N (MXORD+6)*N + (MXORD+5)*N 320C + 2*NROOT + 2*NROOT + 2*NROOT + 2*NROOT 321C + 250 + 250 + 250 + 250 322C --------------------------------------------------------- 323C 324C IWORK 325C LENIW = (Input) 326C IWORK is an integer array of length LENIW used internally 327C for temporary storage. The user must allocate space for 328C this array in the calling program by a statement such as 329C INTEGER IWORK(...) 330C The length of IWORK should be at least 331C 50 if MITER is 0 or 3, or 332C N+50 if MITER is 1, 2, 4, or 5, or MINT is 3, 333C and LENIW should be set to the value used. The contents 334C of IWORK should not be disturbed between calls to CDRIV3. 335C 336C JACOBN = A subroutine supplied by the user, if MITER is 1 or 4. 337C If this is the case, the name must be declared EXTERNAL in 338C the user's calling program. Given a system of N 339C differential equations, it is meaningful to speak about 340C the partial derivative of the I-th right hand side with 341C respect to the J-th dependent variable. In general there 342C are N*N such quantities. Often however the equations can 343C be ordered so that the I-th differential equation only 344C involves dependent variables with index near I, e.g., I+1, 345C I-2. Such a system is called banded. If, for all I, the 346C I-th equation depends on at most the variables 347C Y(I-ML), Y(I-ML+1), ... , Y(I), Y(I+1), ... , Y(I+MU) 348C then we call ML+MU+1 the bandwidth of the system. In a 349C banded system many of the partial derivatives above are 350C automatically zero. For the cases MITER = 1, 2, 4, and 5, 351C some of these partials are needed. For the cases 352C MITER = 2 and 5 the necessary derivatives are 353C approximated numerically by CDRIV3, and we only ask the 354C user to tell CDRIV3 the value of ML and MU if the system 355C is banded. For the cases MITER = 1 and 4 the user must 356C derive these partials algebraically and encode them in 357C subroutine JACOBN. By computing these derivatives the 358C user can often save 20-30 per cent of the computing time. 359C Usually, however, the accuracy is not much affected and 360C most users will probably forego this option. The optional 361C user-written subroutine JACOBN has the form: 362C SUBROUTINE JACOBN (N, T, Y, DFDY, MATDIM, ML, MU) 363C COMPLEX Y(*), DFDY(MATDIM,*) 364C . 365C . 366C Calculate values of DFDY 367C . 368C . 369C END (Sample) 370C Here Y is a vector of length at least N. The actual 371C length of Y is determined by the user's declaration in the 372C program which calls CDRIV3. Thus the dimensioning of Y in 373C JACOBN, while required by FORTRAN convention, does not 374C actually allocate any storage. When this subroutine is 375C called, the first N components of Y are intermediate 376C approximations to the solution components. The user 377C should not alter these values. If the system is not 378C banded (MITER=1), the partials of the I-th equation with 379C respect to the J-th dependent function are to be stored in 380C DFDY(I,J). Thus partials of the I-th equation are stored 381C in the I-th row of DFDY. If the system is banded 382C (MITER=4), then the partials of the I-th equation with 383C respect to Y(J) are to be stored in DFDY(K,J), where 384C K=I-J+MU+1 . Normally a return from JACOBN passes control 385C back to CDRIV3. However, if the user would like to abort 386C the calculation, i.e., return control to the program which 387C calls CDRIV3, he should set N to zero. CDRIV3 will signal 388C this by returning a value of NSTATE equal to +8(-8). 389C Altering the value of N in JACOBN has no effect on the 390C value of N in the call sequence of CDRIV3. 391C 392C FA = A subroutine supplied by the user if IMPL is not zero, and 393C MITER is not 3. If so, the name must be declared EXTERNAL 394C in the user's calling program. This subroutine computes 395C the array A, where A*dY(I)/dT = F(Y(I),T). 396C There are three cases: 397C 398C IMPL=1. 399C Subroutine FA is of the form: 400C SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE) 401C COMPLEX Y(*), A(MATDIM,*) 402C . 403C . 404C Calculate ALL values of A 405C . 406C . 407C END (Sample) 408C In this case A is assumed to be a nonsingular matrix, 409C with the same structure as DFDY (see JACOBN description 410C above). Programming considerations prevent complete 411C generality. If MITER is 1 or 2, A is assumed to be full 412C and the user must compute and store all values of 413C A(I,J), I,J=1, ... ,N. If MITER is 4 or 5, A is assumed 414C to be banded with lower and upper half bandwidth ML and 415C MU. The left hand side of the I-th equation is a linear 416C combination of dY(I-ML)/dT, dY(I-ML+1)/dT, ... , 417C dY(I)/dT, ... , dY(I+MU-1)/dT, dY(I+MU)/dT. Thus in the 418C I-th equation, the coefficient of dY(J)/dT is to be 419C stored in A(K,J), where K=I-J+MU+1. 420C NOTE: The array A will be altered between calls to FA. 421C 422C IMPL=2. 423C Subroutine FA is of the form: 424C SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE) 425C COMPLEX Y(*), A(*) 426C . 427C . 428C Calculate non-zero values of A(1),...,A(NDE) 429C . 430C . 431C END (Sample) 432C In this case it is assumed that the system is ordered by 433C the user so that the differential equations appear 434C first, and the algebraic equations appear last. The 435C algebraic equations must be written in the form: 436C 0 = F(Y(I),T). When using this option it is up to the 437C user to provide initial values for the Y(I) that satisfy 438C the algebraic equations as well as possible. It is 439C further assumed that A is a vector of length NDE. All 440C of the components of A, which may depend on T, Y(I), 441C etc., must be set by the user to non-zero values. 442C 443C IMPL=3. 444C Subroutine FA is of the form: 445C SUBROUTINE FA (N, T, Y, A, MATDIM, ML, MU, NDE) 446C COMPLEX Y(*), A(MATDIM,*) 447C . 448C . 449C Calculate ALL values of A 450C . 451C . 452C END (Sample) 453C In this case A is assumed to be a nonsingular NDE by NDE 454C matrix with the same structure as DFDY (see JACOBN 455C description above). Programming considerations prevent 456C complete generality. If MITER is 1 or 2, A is assumed 457C to be full and the user must compute and store all 458C values of A(I,J), I,J=1, ... ,NDE. If MITER is 4 or 5, 459C A is assumed to be banded with lower and upper half 460C bandwidths ML and MU. The left hand side of the I-th 461C equation is a linear combination of dY(I-ML)/dT, 462C dY(I-ML+1)/dT, ... , dY(I)/dT, ... , dY(I+MU-1)/dT, 463C dY(I+MU)/dT. Thus in the I-th equation, the coefficient 464C of dY(J)/dT is to be stored in A(K,J), where K=I-J+MU+1. 465C It is assumed that the system is ordered by the user so 466C that the differential equations appear first, and the 467C algebraic equations appear last. The algebraic 468C equations must be written in the form 0 = F(Y(I),T). 469C When using this option it is up to the user to provide 470C initial values for the Y(I) that satisfy the algebraic 471C equations as well as possible. 472C NOTE: For IMPL = 3, the array A will be altered between 473C calls to FA. 474C Here Y is a vector of length at least N. The actual 475C length of Y is determined by the user's declaration in the 476C program which calls CDRIV3. Thus the dimensioning of Y in 477C FA, while required by FORTRAN convention, does not 478C actually allocate any storage. When this subroutine is 479C called, the first N components of Y are intermediate 480C approximations to the solution components. The user 481C should not alter these values. FA is always called 482C immediately after calling F, with the same values of T 483C and Y. Normally a return from FA passes control back to 484C CDRIV3. However, if the user would like to abort the 485C calculation, i.e., return control to the program which 486C calls CDRIV3, he should set N to zero. CDRIV3 will signal 487C this by returning a value of NSTATE equal to +9(-9). 488C Altering the value of N in FA has no effect on the value 489C of N in the call sequence of CDRIV3. 490C 491C NDE = (Input) The number of differential equations. This is 492C required only for IMPL = 2 or 3, with NDE .LT. N. 493C 494C MXSTEP = (Input) The maximum number of internal steps allowed on 495C one call to CDRIV3. 496C 497C G = A real FORTRAN function supplied by the user 498C if NROOT is not 0. In this case, the name must be 499C declared EXTERNAL in the user's calling program. G is 500C repeatedly called with different values of IROOT to obtain 501C the value of each of the NROOT equations for which a root 502C is desired. G is of the form: 503C REAL FUNCTION G (N, T, Y, IROOT) 504C COMPLEX Y(*) 505C GO TO (10, ...), IROOT 506C 10 G = ... 507C . 508C . 509C END (Sample) 510C Here, Y is a vector of length at least N, whose first N 511C components are the solution components at the point T. 512C The user should not alter these values. The actual length 513C of Y is determined by the user's declaration in the 514C program which calls CDRIV3. Thus the dimensioning of Y in 515C G, while required by FORTRAN convention, does not actually 516C allocate any storage. Normally a return from G passes 517C control back to CDRIV3. However, if the user would like 518C to abort the calculation, i.e., return control to the 519C program which calls CDRIV3, he should set N to zero. 520C CDRIV3 will signal this by returning a value of NSTATE 521C equal to +7(-7). In this case, the index of the equation 522C being evaluated is stored in the sixth element of IWORK. 523C Altering the value of N in G has no effect on the value of 524C N in the call sequence of CDRIV3. 525C 526C USERS = A subroutine supplied by the user, if MITER is 3. 527C If this is the case, the name must be declared EXTERNAL in 528C the user's calling program. The routine USERS is called 529C by CDRIV3 when certain linear systems must be solved. The 530C user may choose any method to form, store and solve these 531C systems in order to obtain the solution result that is 532C returned to CDRIV3. In particular, this allows sparse 533C matrix methods to be used. The call sequence for this 534C routine is: 535C 536C SUBROUTINE USERS (Y, YH, YWT, SAVE1, SAVE2, T, H, EL, 537C 8 IMPL, N, NDE, IFLAG) 538C COMPLEX Y(*), YH(*), YWT(*), SAVE1(*), SAVE2(*) 539C REAL T, H, EL 540C 541C The input variable IFLAG indicates what action is to be 542C taken. Subroutine USERS should perform the following 543C operations, depending on the value of IFLAG and IMPL. 544C 545C IFLAG = 0 546C IMPL = 0. USERS is not called. 547C IMPL = 1, 2 or 3. Solve the system A*X = SAVE2, 548C returning the result in SAVE2. The array SAVE1 can 549C be used as a work array. For IMPL = 1, there are N 550C components to the system, and for IMPL = 2 or 3, 551C there are NDE components to the system. 552C 553C IFLAG = 1 554C IMPL = 0. Compute, decompose and store the matrix 555C (I - H*EL*J), where I is the identity matrix and J 556C is the Jacobian matrix of the right hand side. The 557C array SAVE1 can be used as a work array. 558C IMPL = 1, 2 or 3. Compute, decompose and store the 559C matrix (A - H*EL*J). The array SAVE1 can be used as 560C a work array. 561C 562C IFLAG = 2 563C IMPL = 0. Solve the system 564C (I - H*EL*J)*X = H*SAVE2 - YH - SAVE1, 565C returning the result in SAVE2. 566C IMPL = 1, 2 or 3. Solve the system 567C (A - H*EL*J)*X = H*SAVE2 - A*(YH + SAVE1) 568C returning the result in SAVE2. 569C The array SAVE1 should not be altered. 570C If IFLAG is 0 and IMPL is 1 or 2 and the matrix A is 571C singular, or if IFLAG is 1 and one of the matrices 572C (I - H*EL*J), (A - H*EL*J) is singular, the INTEGER 573C variable IFLAG is to be set to -1 before RETURNing. 574C Normally a return from USERS passes control back to 575C CDRIV3. However, if the user would like to abort the 576C calculation, i.e., return control to the program which 577C calls CDRIV3, he should set N to zero. CDRIV3 will signal 578C this by returning a value of NSTATE equal to +10(-10). 579C Altering the value of N in USERS has no effect on the 580C value of N in the call sequence of CDRIV3. 581C 582C IERFLG = An error flag. The error number associated with a 583C diagnostic message (see Section III-A below) is the same 584C as the corresponding value of IERFLG. The meaning of 585C IERFLG: 586C 0 The routine completed successfully. (No message is 587C issued.) 588C 3 (Warning) The number of steps required to reach TOUT 589C exceeds MXSTEP. 590C 4 (Warning) The value of EPS is too small. 591C 11 (Warning) For NTASK = 2 or 3, T is beyond TOUT. 592C The solution was obtained by interpolation. 593C 15 (Warning) The integration step size is below the 594C roundoff level of T. (The program issues this 595C message as a warning but does not return control to 596C the user.) 597C 22 (Recoverable) N is not positive. 598C 23 (Recoverable) MINT is less than 1 or greater than 3 . 599C 24 (Recoverable) MITER is less than 0 or greater than 600C 5 . 601C 25 (Recoverable) IMPL is less than 0 or greater than 3 . 602C 26 (Recoverable) The value of NSTATE is less than 1 or 603C greater than 12 . 604C 27 (Recoverable) EPS is less than zero. 605C 28 (Recoverable) MXORD is not positive. 606C 29 (Recoverable) For MINT = 3, either MITER = 0 or 3, or 607C IMPL = 0 . 608C 30 (Recoverable) For MITER = 0, IMPL is not 0 . 609C 31 (Recoverable) For MINT = 1, IMPL is 2 or 3 . 610C 32 (Recoverable) Insufficient storage has been allocated 611C for the WORK array. 612C 33 (Recoverable) Insufficient storage has been allocated 613C for the IWORK array. 614C 41 (Recoverable) The integration step size has gone 615C to zero. 616C 42 (Recoverable) The integration step size has been 617C reduced about 50 times without advancing the 618C solution. The problem setup may not be correct. 619C 43 (Recoverable) For IMPL greater than 0, the matrix A 620C is singular. 621C 999 (Fatal) The value of NSTATE is 12 . 622C 623C III. OTHER COMMUNICATION TO THE USER .............................. 624C 625C A. The solver communicates to the user through the parameters 626C above. In addition it writes diagnostic messages through the 627C standard error handling program XERMSG. A complete description 628C of XERMSG is given in "Guide to the SLATEC Common Mathematical 629C Library" by Kirby W. Fong et al.. At installations which do not 630C have this error handling package the short but serviceable 631C routine, XERMSG, available with this package, can be used. That 632C program uses the file named OUTPUT to transmit messages. 633C 634C B. The first three elements of WORK and the first five elements of 635C IWORK will contain the following statistical data: 636C AVGH The average step size used. 637C HUSED The step size last used (successfully). 638C AVGORD The average order used. 639C IMXERR The index of the element of the solution vector that 640C contributed most to the last error test. 641C NQUSED The order last used (successfully). 642C NSTEP The number of steps taken since last initialization. 643C NFE The number of evaluations of the right hand side. 644C NJE The number of evaluations of the Jacobian matrix. 645C 646C IV. REMARKS ....................................................... 647C 648C A. Other routines used: 649C CDNTP, CDZRO, CDSTP, CDNTL, CDPST, CDCOR, CDCST, 650C CDPSC, and CDSCL; 651C CGEFA, CGESL, CGBFA, CGBSL, and SCNRM2 (from LINPACK) 652C R1MACH (from the Bell Laboratories Machine Constants Package) 653C XERMSG (from the SLATEC Common Math Library) 654C The last seven routines above, not having been written by the 655C present authors, are not explicitly part of this package. 656C 657C B. On any return from CDRIV3 all information necessary to continue 658C the calculation is contained in the call sequence parameters, 659C including the work arrays. Thus it is possible to suspend one 660C problem, integrate another, and then return to the first. 661C 662C C. If this package is to be used in an overlay situation, the user 663C must declare in the primary overlay the variables in the call 664C sequence to CDRIV3. 665C 666C D. Changing parameters during an integration. 667C The value of NROOT, EPS, EWT, IERROR, MINT, MITER, or HMAX may 668C be altered by the user between calls to CDRIV3. For example, if 669C too much accuracy has been requested (the program returns with 670C NSTATE = 4 and an increased value of EPS) the user may wish to 671C increase EPS further. In general, prudence is necessary when 672C making changes in parameters since such changes are not 673C implemented until the next integration step, which is not 674C necessarily the next call to CDRIV3. This can happen if the 675C program has already integrated to a point which is beyond the 676C new point TOUT. 677C 678C E. As the price for complete control of matrix algebra, the CDRIV3 679C USERS option puts all responsibility for Jacobian matrix 680C evaluation on the user. It is often useful to approximate 681C numerically all or part of the Jacobian matrix. However this 682C must be done carefully. The FORTRAN sequence below illustrates 683C the method we recommend. It can be inserted directly into 684C subroutine USERS to approximate Jacobian elements in rows I1 685C to I2 and columns J1 to J2. 686C COMPLEX DFDY(N,N), R, SAVE1(N), SAVE2(N), Y(N), YJ, YWT(N) 687C REAL EPSJ, H, R1MACH, T, UROUND 688C UROUND = R1MACH(4) 689C EPSJ = SQRT(UROUND) 690C DO 30 J = J1,J2 691C IF (ABS(Y(J)) .GT. ABS(YWT(J))) THEN 692C R = EPSJ*Y(J) 693C ELSE 694C R = EPSJ*YWT(J) 695C END IF 696C IF (R .EQ. 0.E0) R = YWT(J) 697C YJ = Y(J) 698C Y(J) = Y(J) + R 699C CALL F (N, T, Y, SAVE1) 700C IF (N .EQ. 0) RETURN 701C Y(J) = YJ 702C DO 20 I = I1,I2 703C 20 DFDY(I,J) = (SAVE1(I) - SAVE2(I))/R 704C 30 CONTINUE 705C Many problems give rise to structured sparse Jacobians, e.g., 706C block banded. It is possible to approximate them with fewer 707C function evaluations than the above procedure uses; see Curtis, 708C Powell and Reid, J. Inst. Maths Applics, (1974), Vol. 13, 709C pp. 117-119. 710C 711C F. When any of the routines JACOBN, FA, G, or USERS, is not 712C required, difficulties associated with unsatisfied externals can 713C be avoided by using the name of the routine which calculates the 714C right hand side of the differential equations in place of the 715C corresponding name in the call sequence of CDRIV3. 716C 717C***REFERENCES C. W. Gear, Numerical Initial Value Problems in 718C Ordinary Differential Equations, Prentice-Hall, 1971. 719C***ROUTINES CALLED CDNTP, CDSTP, CDZRO, CGBFA, CGBSL, CGEFA, CGESL, 720C R1MACH, SCNRM2, XERMSG 721C***REVISION HISTORY (YYMMDD) 722C 790601 DATE WRITTEN 723C 900329 Initial submission to SLATEC. 724C***END PROLOGUE CDRIV3 725 EXTERNAL F, JACOBN, FA, G, USERS 726 COMPLEX WORK(*), Y(*) 727 REAL AE, AVGH, AVGORD, BIG, EL(13,12), EPS, EWT(*), 728 8 G, GLAST, GNOW, H, HMAX, HOLD, HSIGN, HUSED, NROUND, RC, RE, 729 8 RMAX, R1MACH, SCNRM2, SIZE, SUM, T, TLAST, TOUT, TQ(3,12), 730 8 TREND, TROOT, UROUND 731 INTEGER I, IA, IAVGH, IAVGRD, ICNVRG, IDFDY, IEL, IERFLG, IERROR, 732 8 IFAC, IFLAG, IGNOW, IH, IHMAX, IHOLD, IHSIGN, IHUSED, 733 8 IJROOT, IJSTPL, IJTASK, IMNT, IMNTLD, IMPL, IMTR, IMTRLD, 734 8 IMTRSV, IMXERR, IMXORD, IMXRDS, INDMXR, INDPRT, INDPVT, 735 8 INDTRT, INFE, INFO, INJE, INQ, INQUSE, INROOT, INRTLD, 736 8 INSTEP, INWAIT, IRC, IRMAX, IROOT, IMACH1, IMACH4, ISAVE1, 737 8 ISAVE2, IT, ITOUT, ITQ, ITREND, ITROOT, IWORK(*), IYH, 738 8 IYWT, J, JSTATE, JTROOT, LENCHK, LENIW, LENW, LIWCHK, 739 8 MATDIM, MAXORD, MINT, MITER, ML, MU, MXORD, MXSTEP, N, 740 8 NDE, NDECOM, NPAR, NROOT, NSTATE, NSTEPL, NTASK 741 LOGICAL CONVRG 742 CHARACTER INTGR1*8, INTGR2*8, RL1*16, RL2*16 743 PARAMETER(NROUND = 20.E0) 744 PARAMETER(IAVGH = 1, IHUSED = 2, IAVGRD = 3, 745 8 IEL = 4, IH = 160, IHMAX = 161, IHOLD = 162, 746 8 IHSIGN = 163, IRC = 164, IRMAX = 165, IT = 166, 747 8 ITOUT = 167, ITQ = 168, ITREND = 204, IMACH1 = 205, 748 8 IMACH4 = 206, IYH = 251, 749 8 INDMXR = 1, INQUSE = 2, INSTEP = 3, INFE = 4, INJE = 5, 750 8 INROOT = 6, ICNVRG = 7, IJROOT = 8, IJTASK = 9, 751 8 IMNTLD = 10, IMTRLD = 11, INQ = 12, INRTLD = 13, 752 8 INDTRT = 14, INWAIT = 15, IMNT = 16, IMTRSV = 17, 753 8 IMTR = 18, IMXRDS = 19, IMXORD = 20, INDPRT = 21, 754 8 IJSTPL = 22, INDPVT = 51) 755C***FIRST EXECUTABLE STATEMENT CDRIV3 756 IF (NSTATE .EQ. 12) THEN 757 IERFLG = 999 758 CALL XERMSG('SLATEC', 'CDRIV3', 759 8 'Illegal input. The value of NSTATE is 12 .', IERFLG, 2) 760 RETURN 761 ELSE IF (NSTATE .LT. 1 .OR. NSTATE .GT. 12) THEN 762 WRITE(INTGR1, '(I8)') NSTATE 763 IERFLG = 26 764 CALL XERMSG('SLATEC', 'CDRIV3', 765 8 'Illegal input. Improper value for NSTATE(= '//INTGR1//').', 766 8 IERFLG, 1) 767 NSTATE = 12 768 RETURN 769 END IF 770 NPAR = N 771 IF (EPS .LT. 0.E0) THEN 772 WRITE(RL1, '(E16.8)') EPS 773 IERFLG = 27 774 CALL XERMSG('SLATEC', 'CDRIV3', 775 8 'Illegal input. EPS, '//RL1//', is negative.', IERFLG, 1) 776 NSTATE = 12 777 RETURN 778 END IF 779 IF (N .LE. 0) THEN 780 WRITE(INTGR1, '(I8)') N 781 IERFLG = 22 782 CALL XERMSG('SLATEC', 'CDRIV3', 783 8 'Illegal input. Number of equations, '//INTGR1// 784 8 ', is not positive.', IERFLG, 1) 785 NSTATE = 12 786 RETURN 787 END IF 788 IF (MXORD .LE. 0) THEN 789 WRITE(INTGR1, '(I8)') MXORD 790 IERFLG = 28 791 CALL XERMSG('SLATEC', 'CDRIV3', 792 8 'Illegal input. Maximum order, '//INTGR1// 793 8 ', is not positive.', IERFLG, 1) 794 NSTATE = 12 795 RETURN 796 END IF 797 IF (MINT .LT. 1 .OR. MINT .GT. 3) THEN 798 WRITE(INTGR1, '(I8)') MINT 799 IERFLG = 23 800 CALL XERMSG('SLATEC', 'CDRIV3', 801 8 'Illegal input. Improper value for the integration method '// 802 8 'flag, '//INTGR1//' .', IERFLG, 1) 803 NSTATE = 12 804 RETURN 805 ELSE IF (MITER .LT. 0 .OR. MITER .GT. 5) THEN 806 WRITE(INTGR1, '(I8)') MITER 807 IERFLG = 24 808 CALL XERMSG('SLATEC', 'CDRIV3', 809 8 'Illegal input. Improper value for MITER(= '//INTGR1//').', 810 8 IERFLG, 1) 811 NSTATE = 12 812 RETURN 813 ELSE IF (IMPL .LT. 0 .OR. IMPL .GT. 3) THEN 814 WRITE(INTGR1, '(I8)') IMPL 815 IERFLG = 25 816 CALL XERMSG('SLATEC', 'CDRIV3', 817 8 'Illegal input. Improper value for IMPL(= '//INTGR1//').', 818 8 IERFLG, 1) 819 NSTATE = 12 820 RETURN 821 ELSE IF (MINT .EQ. 3 .AND. 822 8 (MITER .EQ. 0 .OR. MITER .EQ. 3 .OR. IMPL .NE. 0)) THEN 823 WRITE(INTGR1, '(I8)') MITER 824 WRITE(INTGR2, '(I8)') IMPL 825 IERFLG = 29 826 CALL XERMSG('SLATEC', 'CDRIV3', 827 8 'Illegal input. For MINT = 3, the value of MITER, '//INTGR1// 828 8 ', and/or IMPL, '//INTGR2//', is not allowed.', IERFLG, 1) 829 NSTATE = 12 830 RETURN 831 ELSE IF ((IMPL .GE. 1 .AND. IMPL .LE. 3) .AND. MITER .EQ. 0) THEN 832 WRITE(INTGR1, '(I8)') IMPL 833 IERFLG = 30 834 CALL XERMSG('SLATEC', 'CDRIV3', 835 8 'Illegal input. For MITER = 0, the value of IMPL, '//INTGR1// 836 8 ', is not allowed.', IERFLG, 1) 837 NSTATE = 12 838 RETURN 839 ELSE IF ((IMPL .EQ. 2 .OR. IMPL .EQ. 3) .AND. MINT .EQ. 1) THEN 840 WRITE(INTGR1, '(I8)') IMPL 841 IERFLG = 31 842 CALL XERMSG('SLATEC', 'CDRIV3', 843 8 'Illegal input. For MINT = 1, the value of IMPL, '//INTGR1// 844 8 ', is not allowed.', IERFLG, 1) 845 NSTATE = 12 846 RETURN 847 END IF 848 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) THEN 849 LIWCHK = INDPVT - 1 850 ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2 .OR. MITER .EQ. 4 .OR. 851 8 MITER .EQ. 5) THEN 852 LIWCHK = INDPVT + N - 1 853 END IF 854 IF (LENIW .LT. LIWCHK) THEN 855 WRITE(INTGR1, '(I8)') LIWCHK 856 IERFLG = 33 857 CALL XERMSG('SLATEC', 'CDRIV3', 858 8 'Illegal input. Insufficient storage allocated for the '// 859 8 'IWORK array. Based on the value of the input parameters '// 860 8 'involved, the required storage is '//INTGR1//' .', IERFLG, 1) 861 NSTATE = 12 862 RETURN 863 END IF 864C Allocate the WORK array 865C IYH is the index of YH in WORK 866 IF (MINT .EQ. 1 .OR. MINT .EQ. 3) THEN 867 MAXORD = MIN(MXORD, 12) 868 ELSE IF (MINT .EQ. 2) THEN 869 MAXORD = MIN(MXORD, 5) 870 END IF 871 IDFDY = IYH + (MAXORD + 1)*N 872C IDFDY is the index of DFDY 873C 874 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) THEN 875 IYWT = IDFDY 876 ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN 877 IYWT = IDFDY + N*N 878 ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN 879 IYWT = IDFDY + (2*ML + MU + 1)*N 880 END IF 881C IYWT is the index of YWT 882 ISAVE1 = IYWT + N 883C ISAVE1 is the index of SAVE1 884 ISAVE2 = ISAVE1 + N 885C ISAVE2 is the index of SAVE2 886 IGNOW = ISAVE2 + N 887C IGNOW is the index of GNOW 888 ITROOT = IGNOW + NROOT 889C ITROOT is the index of TROOT 890 IFAC = ITROOT + NROOT 891C IFAC is the index of FAC 892 IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. MINT .EQ. 3) THEN 893 IA = IFAC + N 894 ELSE 895 IA = IFAC 896 END IF 897C IA is the index of A 898 IF (IMPL .EQ. 0 .OR. MITER .EQ. 3) THEN 899 LENCHK = IA - 1 900 ELSE IF (IMPL .EQ. 1 .AND. (MITER .EQ. 1 .OR. MITER .EQ. 2)) THEN 901 LENCHK = IA - 1 + N*N 902 ELSE IF (IMPL .EQ. 1 .AND. (MITER .EQ. 4 .OR. MITER .EQ. 5)) THEN 903 LENCHK = IA - 1 + (2*ML + MU + 1)*N 904 ELSE IF (IMPL .EQ. 2 .AND. MITER .NE. 3) THEN 905 LENCHK = IA - 1 + N 906 ELSE IF (IMPL .EQ. 3 .AND. (MITER .EQ. 1 .OR. MITER .EQ. 2)) THEN 907 LENCHK = IA - 1 + N*NDE 908 ELSE IF (IMPL .EQ. 3 .AND. (MITER .EQ. 4 .OR. MITER .EQ. 5)) THEN 909 LENCHK = IA - 1 + (2*ML + MU + 1)*NDE 910 END IF 911 IF (LENW .LT. LENCHK) THEN 912 WRITE(INTGR1, '(I8)') LENCHK 913 IERFLG = 32 914 CALL XERMSG('SLATEC', 'CDRIV3', 915 8 'Illegal input. Insufficient storage allocated for the '// 916 8 'WORK array. Based on the value of the input parameters '// 917 8 'involved, the required storage is '//INTGR1//' .', IERFLG, 1) 918 NSTATE = 12 919 RETURN 920 END IF 921 IF (MITER .EQ. 0 .OR. MITER .EQ. 3) THEN 922 MATDIM = 1 923 ELSE IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN 924 MATDIM = N 925 ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN 926 MATDIM = 2*ML + MU + 1 927 END IF 928 IF (IMPL .EQ. 0 .OR. IMPL .EQ. 1) THEN 929 NDECOM = N 930 ELSE IF (IMPL .EQ. 2 .OR. IMPL .EQ. 3) THEN 931 NDECOM = NDE 932 END IF 933 IF (NSTATE .EQ. 1) THEN 934C Initialize parameters 935 IF (MINT .EQ. 1 .OR. MINT .EQ. 3) THEN 936 IWORK(IMXORD) = MIN(MXORD, 12) 937 ELSE IF (MINT .EQ. 2) THEN 938 IWORK(IMXORD) = MIN(MXORD, 5) 939 END IF 940 IWORK(IMXRDS) = MXORD 941 IF (MINT .EQ. 1 .OR. MINT .EQ. 2) THEN 942 IWORK(IMNT) = MINT 943 IWORK(IMTR) = MITER 944 IWORK(IMNTLD) = MINT 945 IWORK(IMTRLD) = MITER 946 ELSE IF (MINT .EQ. 3) THEN 947 IWORK(IMNT) = 1 948 IWORK(IMTR) = 0 949 IWORK(IMNTLD) = IWORK(IMNT) 950 IWORK(IMTRLD) = IWORK(IMTR) 951 IWORK(IMTRSV) = MITER 952 END IF 953 WORK(IHMAX) = HMAX 954 UROUND = R1MACH (4) 955 WORK(IMACH4) = UROUND 956 WORK(IMACH1) = R1MACH (1) 957 IF (NROOT .NE. 0) THEN 958 RE = UROUND 959 AE = WORK(IMACH1) 960 END IF 961 H = (TOUT - T)*(1.E0 - 4.E0*UROUND) 962 H = SIGN(MIN(ABS(H), HMAX), H) 963 WORK(IH) = H 964 HSIGN = SIGN(1.E0, H) 965 WORK(IHSIGN) = HSIGN 966 IWORK(IJTASK) = 0 967 AVGH = 0.E0 968 AVGORD = 0.E0 969 WORK(IAVGH) = 0.E0 970 WORK(IHUSED) = 0.E0 971 WORK(IAVGRD) = 0.E0 972 IWORK(INDMXR) = 0 973 IWORK(INQUSE) = 0 974 IWORK(INSTEP) = 0 975 IWORK(IJSTPL) = 0 976 IWORK(INFE) = 0 977 IWORK(INJE) = 0 978 IWORK(INROOT) = 0 979 WORK(IT) = T 980 IWORK(ICNVRG) = 0 981 IWORK(INDPRT) = 0 982C Set initial conditions 983 DO 30 I = 1,N 984 30 WORK(I+IYH-1) = Y(I) 985 IF (T .EQ. TOUT) RETURN 986 GO TO 180 987 ELSE 988 UROUND = WORK(IMACH4) 989 IF (NROOT .NE. 0) THEN 990 RE = UROUND 991 AE = WORK(IMACH1) 992 END IF 993 END IF 994C On a continuation, check 995C that output points have 996C been or will be overtaken. 997 IF (IWORK(ICNVRG) .EQ. 1) THEN 998 CONVRG = .TRUE. 999 ELSE 1000 CONVRG = .FALSE. 1001 END IF 1002 AVGH = WORK(IAVGH) 1003 AVGORD = WORK(IAVGRD) 1004 HOLD = WORK(IHOLD) 1005 RC = WORK(IRC) 1006 RMAX = WORK(IRMAX) 1007 TREND = WORK(ITREND) 1008 DO 35 J = 1,12 1009 DO 35 I = 1,13 1010 35 EL(I,J) = WORK(I+IEL+(J-1)*13-1) 1011 DO 40 J = 1,12 1012 DO 40 I = 1,3 1013 40 TQ(I,J) = WORK(I+ITQ+(J-1)*3-1) 1014 T = WORK(IT) 1015 H = WORK(IH) 1016 HSIGN = WORK(IHSIGN) 1017 IF (IWORK(IJTASK) .EQ. 0) GO TO 180 1018C 1019C IWORK(IJROOT) flags unreported 1020C roots, and is set to the value of 1021C NTASK when a root was last selected. 1022C It is set to zero when all roots 1023C have been reported. IWORK(INROOT) 1024C contains the index and WORK(ITOUT) 1025C contains the value of the root last 1026C selected to be reported. 1027C IWORK(INRTLD) contains the value of 1028C NROOT and IWORK(INDTRT) contains 1029C the value of ITROOT when the array 1030C of roots was last calculated. 1031 IF (NROOT .NE. 0) THEN 1032 IF (IWORK(IJROOT) .GT. 0) THEN 1033C TOUT has just been reported. 1034C If TROOT .LE. TOUT, report TROOT. 1035 IF (NSTATE .NE. 5) THEN 1036 IF (TOUT*HSIGN .GE. REAL(WORK(ITOUT))*HSIGN) THEN 1037 TROOT = WORK(ITOUT) 1038 CALL CDNTP (H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), Y) 1039 T = TROOT 1040 NSTATE = 5 1041 IERFLG = 0 1042 GO TO 580 1043 END IF 1044C A root has just been reported. 1045C Select the next root. 1046 ELSE 1047 TROOT = T 1048 IROOT = 0 1049 DO 50 I = 1,IWORK(INRTLD) 1050 JTROOT = I + IWORK(INDTRT) - 1 1051 IF (REAL(WORK(JTROOT))*HSIGN .LE. TROOT*HSIGN) THEN 1052C 1053C Check for multiple roots. 1054C 1055 IF (WORK(JTROOT) .EQ. WORK(ITOUT) .AND. 1056 8 I .GT. IWORK(INROOT)) THEN 1057 IROOT = I 1058 TROOT = WORK(JTROOT) 1059 GO TO 60 1060 END IF 1061 IF (REAL(WORK(JTROOT))*HSIGN .GT. 1062 8 REAL(WORK(ITOUT))*HSIGN) THEN 1063 IROOT = I 1064 TROOT = WORK(JTROOT) 1065 END IF 1066 END IF 1067 50 CONTINUE 1068 60 IWORK(INROOT) = IROOT 1069 WORK(ITOUT) = TROOT 1070 IWORK(IJROOT) = NTASK 1071 IF (NTASK .EQ. 1) THEN 1072 IF (IROOT .EQ. 0) THEN 1073 IWORK(IJROOT) = 0 1074 ELSE 1075 IF (TOUT*HSIGN .GE. TROOT*HSIGN) THEN 1076 CALL CDNTP (H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), 1077 8 Y) 1078 NSTATE = 5 1079 T = TROOT 1080 IERFLG = 0 1081 GO TO 580 1082 END IF 1083 END IF 1084 ELSE IF (NTASK .EQ. 2 .OR. NTASK .EQ. 3) THEN 1085C 1086C If there are no more roots, or the 1087C user has altered TOUT to be less 1088C than a root, set IJROOT to zero. 1089C 1090 IF (IROOT .EQ. 0 .OR. (TOUT*HSIGN .LT. TROOT*HSIGN)) THEN 1091 IWORK(IJROOT) = 0 1092 ELSE 1093 CALL CDNTP (H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), 1094 8 Y) 1095 NSTATE = 5 1096 T = TROOT 1097 IERFLG = 0 1098 GO TO 580 1099 END IF 1100 END IF 1101 END IF 1102 END IF 1103 END IF 1104C 1105 IF (NTASK .EQ. 1) THEN 1106 NSTATE = 2 1107 IF (T*HSIGN .GE. TOUT*HSIGN) THEN 1108 CALL CDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y) 1109 T = TOUT 1110 IERFLG = 0 1111 GO TO 580 1112 END IF 1113 ELSE IF (NTASK .EQ. 2) THEN 1114C Check if TOUT has 1115C been reset .LT. T 1116 IF (T*HSIGN .GT. TOUT*HSIGN) THEN 1117 WRITE(RL1, '(E16.8)') T 1118 WRITE(RL2, '(E16.8)') TOUT 1119 IERFLG = 11 1120 CALL XERMSG('SLATEC', 'CDRIV3', 1121 8 'While integrating exactly to TOUT, T, '//RL1// 1122 8 ', was beyond TOUT, '//RL2//' . Solution obtained by '// 1123 8 'interpolation.', IERFLG, 0) 1124 NSTATE = 11 1125 CALL CDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y) 1126 T = TOUT 1127 GO TO 580 1128 END IF 1129C Determine if TOUT has been overtaken 1130C 1131 IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN 1132 T = TOUT 1133 NSTATE = 2 1134 IERFLG = 0 1135 GO TO 560 1136 END IF 1137C If there are no more roots 1138C to report, report T. 1139 IF (NSTATE .EQ. 5) THEN 1140 NSTATE = 2 1141 IERFLG = 0 1142 GO TO 560 1143 END IF 1144 NSTATE = 2 1145C See if TOUT will 1146C be overtaken. 1147 IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN 1148 H = TOUT - T 1149 IF ((T + H)*HSIGN .GT. TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND) 1150 WORK(IH) = H 1151 IF (H .EQ. 0.E0) GO TO 670 1152 IWORK(IJTASK) = -1 1153 END IF 1154 ELSE IF (NTASK .EQ. 3) THEN 1155 NSTATE = 2 1156 IF (T*HSIGN .GT. TOUT*HSIGN) THEN 1157 WRITE(RL1, '(E16.8)') T 1158 WRITE(RL2, '(E16.8)') TOUT 1159 IERFLG = 11 1160 CALL XERMSG('SLATEC', 'CDRIV3', 1161 8 'While integrating exactly to TOUT, T, '//RL1// 1162 8 ', was beyond TOUT, '//RL2//' . Solution obtained by '// 1163 8 'interpolation.', IERFLG, 0) 1164 NSTATE = 11 1165 CALL CDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y) 1166 T = TOUT 1167 GO TO 580 1168 END IF 1169 IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN 1170 T = TOUT 1171 IERFLG = 0 1172 GO TO 560 1173 END IF 1174 IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN 1175 H = TOUT - T 1176 IF ((T + H)*HSIGN .GT. TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND) 1177 WORK(IH) = H 1178 IF (H .EQ. 0.E0) GO TO 670 1179 IWORK(IJTASK) = -1 1180 END IF 1181 END IF 1182C Implement changes in MINT, MITER, and/or HMAX. 1183C 1184 IF ((MINT .NE. IWORK(IMNTLD) .OR. MITER .NE. IWORK(IMTRLD)) .AND. 1185 8 MINT .NE. 3 .AND. IWORK(IMNTLD) .NE. 3) IWORK(IJTASK) = -1 1186 IF (HMAX .NE. WORK(IHMAX)) THEN 1187 H = SIGN(MIN(ABS(H), HMAX), H) 1188 IF (H .NE. WORK(IH)) THEN 1189 IWORK(IJTASK) = -1 1190 WORK(IH) = H 1191 END IF 1192 WORK(IHMAX) = HMAX 1193 END IF 1194C 1195 180 NSTEPL = IWORK(INSTEP) 1196 DO 190 I = 1,N 1197 190 Y(I) = WORK(I+IYH-1) 1198 IF (NROOT .NE. 0) THEN 1199 DO 200 I = 1,NROOT 1200 WORK(I+IGNOW-1) = G (NPAR, T, Y, I) 1201 IF (NPAR .EQ. 0) THEN 1202 IWORK(INROOT) = I 1203 NSTATE = 7 1204 RETURN 1205 END IF 1206 200 CONTINUE 1207 END IF 1208 IF (IERROR .EQ. 1) THEN 1209 DO 230 I = 1,N 1210 230 WORK(I+IYWT-1) = 1.E0 1211 GO TO 410 1212 ELSE IF (IERROR .EQ. 5) THEN 1213 DO 250 I = 1,N 1214 250 WORK(I+IYWT-1) = EWT(I) 1215 GO TO 410 1216 END IF 1217C Reset YWT array. Looping point. 1218 260 IF (IERROR .EQ. 2) THEN 1219 DO 280 I = 1,N 1220 IF (Y(I) .EQ. 0.E0) GO TO 290 1221 280 WORK(I+IYWT-1) = Y(I) 1222 GO TO 410 1223 290 IF (IWORK(IJTASK) .EQ. 0) THEN 1224 CALL F (NPAR, T, Y, WORK(ISAVE2)) 1225 IF (NPAR .EQ. 0) THEN 1226 NSTATE = 6 1227 RETURN 1228 END IF 1229 IWORK(INFE) = IWORK(INFE) + 1 1230 IF (MITER .EQ. 3 .AND. IMPL .NE. 0) THEN 1231 IFLAG = 0 1232 CALL USERS (Y, WORK(IYH), WORK(IYWT), WORK(ISAVE1), 1233 8 WORK(ISAVE2), T, H, REAL(WORK(IEL)), IMPL, NPAR, 1234 8 NDECOM, IFLAG) 1235 IF (IFLAG .EQ. -1) GO TO 690 1236 IF (NPAR .EQ. 0) THEN 1237 NSTATE = 10 1238 RETURN 1239 END IF 1240 ELSE IF (IMPL .EQ. 1) THEN 1241 IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN 1242 CALL FA (NPAR, T, Y, WORK(IA), MATDIM, ML, MU, NDECOM) 1243 IF (NPAR .EQ. 0) THEN 1244 NSTATE = 9 1245 RETURN 1246 END IF 1247 CALL CGEFA (WORK(IA), MATDIM, N, IWORK(INDPVT), INFO) 1248 IF (INFO .NE. 0) GO TO 690 1249 CALL CGESL (WORK(IA), MATDIM, N, IWORK(INDPVT), 1250 8 WORK(ISAVE2), 0) 1251 ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN 1252 CALL FA (NPAR, T, Y, WORK(IA+ML), MATDIM, ML, MU, NDECOM) 1253 IF (NPAR .EQ. 0) THEN 1254 NSTATE = 9 1255 RETURN 1256 END IF 1257 CALL CGBFA (WORK(IA), MATDIM, N, ML, MU, IWORK(INDPVT), 1258 8 INFO) 1259 IF (INFO .NE. 0) GO TO 690 1260 CALL CGBSL (WORK(IA), MATDIM, N, ML, MU, IWORK(INDPVT), 1261 8 WORK(ISAVE2), 0) 1262 END IF 1263 ELSE IF (IMPL .EQ. 2) THEN 1264 CALL FA (NPAR, T, Y, WORK(IA), MATDIM, ML, MU, NDECOM) 1265 IF (NPAR .EQ. 0) THEN 1266 NSTATE = 9 1267 RETURN 1268 END IF 1269 DO 340 I = 1,NDECOM 1270 IF (WORK(I+IA-1) .EQ. 0.E0) GO TO 690 1271 340 WORK(I+ISAVE2-1) = WORK(I+ISAVE2-1)/WORK(I+IA-1) 1272 ELSE IF (IMPL .EQ. 3) THEN 1273 IF (MITER .EQ. 1 .OR. MITER .EQ. 2) THEN 1274 CALL FA (NPAR, T, Y, WORK(IA), MATDIM, ML, MU, NDECOM) 1275 IF (NPAR .EQ. 0) THEN 1276 NSTATE = 9 1277 RETURN 1278 END IF 1279 CALL CGEFA (WORK(IA), MATDIM, NDE, IWORK(INDPVT), INFO) 1280 IF (INFO .NE. 0) GO TO 690 1281 CALL CGESL (WORK(IA), MATDIM, NDE, IWORK(INDPVT), 1282 8 WORK(ISAVE2), 0) 1283 ELSE IF (MITER .EQ. 4 .OR. MITER .EQ. 5) THEN 1284 CALL FA (NPAR, T, Y, WORK(IA+ML), MATDIM, ML, MU, NDECOM) 1285 IF (NPAR .EQ. 0) THEN 1286 NSTATE = 9 1287 RETURN 1288 END IF 1289 CALL CGBFA (WORK(IA), MATDIM, NDE, ML, MU, IWORK(INDPVT), 1290 8 INFO) 1291 IF (INFO .NE. 0) GO TO 690 1292 CALL CGBSL (WORK(IA), MATDIM, NDE, ML, MU, IWORK(INDPVT), 1293 8 WORK(ISAVE2), 0) 1294 END IF 1295 END IF 1296 END IF 1297 DO 360 J = I,N 1298 IF (Y(J) .NE. 0.E0) THEN 1299 WORK(J+IYWT-1) = Y(J) 1300 ELSE 1301 IF (IWORK(IJTASK) .EQ. 0) THEN 1302 WORK(J+IYWT-1) = H*WORK(J+ISAVE2-1) 1303 ELSE 1304 WORK(J+IYWT-1) = WORK(J+IYH+N-1) 1305 END IF 1306 END IF 1307 IF (WORK(J+IYWT-1) .EQ. 0.E0) WORK(J+IYWT-1) = UROUND 1308 360 CONTINUE 1309 ELSE IF (IERROR .EQ. 3) THEN 1310 DO 380 I = 1,N 1311 380 WORK(I+IYWT-1) = MAX(EWT(1), ABS(Y(I))) 1312 ELSE IF (IERROR .EQ. 4) THEN 1313 DO 400 I = 1,N 1314 400 WORK(I+IYWT-1) = MAX(EWT(I), ABS(Y(I))) 1315 END IF 1316C 1317 410 DO 420 I = 1,N 1318 420 WORK(I+ISAVE2-1) = Y(I)/WORK(I+IYWT-1) 1319 SUM = SCNRM2(N, WORK(ISAVE2), 1)/SQRT(REAL(N)) 1320 SUM = MAX(1.E0, SUM) 1321 IF (EPS .LT. SUM*UROUND) THEN 1322 EPS = SUM*UROUND*(1.E0 + 10.E0*UROUND) 1323 WRITE(RL1, '(E16.8)') T 1324 WRITE(RL2, '(E16.8)') EPS 1325 IERFLG = 4 1326 CALL XERMSG('SLATEC', 'CDRIV3', 1327 8 'At T, '//RL1//', the requested accuracy, EPS, was not '// 1328 8 'obtainable with the machine precision. EPS has been '// 1329 8 'increased to '//RL2//' .', IERFLG, 0) 1330 NSTATE = 4 1331 GO TO 560 1332 END IF 1333 IF (ABS(H) .GE. UROUND*ABS(T)) THEN 1334 IWORK(INDPRT) = 0 1335 ELSE IF (IWORK(INDPRT) .EQ. 0) THEN 1336 WRITE(RL1, '(E16.8)') T 1337 WRITE(RL2, '(E16.8)') H 1338 IERFLG = 15 1339 CALL XERMSG('SLATEC', 'CDRIV3', 1340 8 'At T, '//RL1//', the step size, '//RL2//', is smaller '// 1341 8 'than the roundoff level of T. This may occur if there is '// 1342 8 'an abrupt change in the right hand side of the '// 1343 8 'differential equations.', IERFLG, 0) 1344 IWORK(INDPRT) = 1 1345 END IF 1346 IF (NTASK.NE.2) THEN 1347 IF ((IWORK(INSTEP)-NSTEPL) .EQ. MXSTEP) THEN 1348 WRITE(RL1, '(E16.8)') T 1349 WRITE(INTGR1, '(I8)') MXSTEP 1350 WRITE(RL2, '(E16.8)') TOUT 1351 IERFLG = 3 1352 CALL XERMSG('SLATEC', 'CDRIV3', 1353 8 'At T, '//RL1//', '//INTGR1//' steps have been taken '// 1354 8 'without reaching TOUT, '//RL2//' .', IERFLG, 0) 1355 NSTATE = 3 1356 GO TO 560 1357 END IF 1358 END IF 1359C 1360C CALL CDSTP (EPS, F, FA, HMAX, IMPL, IERROR, JACOBN, MATDIM, 1361C 8 MAXORD, MINT, MITER, ML, MU, N, NDE, YWT, UROUND, 1362C 8 USERS, AVGH, AVGORD, H, HUSED, JTASK, MNTOLD, MTROLD, 1363C 8 NFE, NJE, NQUSED, NSTEP, T, Y, YH, A, CONVRG, 1364C 8 DFDY, EL, FAC, HOLD, IPVT, JSTATE, JSTEPL, NQ, NWAIT, 1365C 8 RC, RMAX, SAVE1, SAVE2, TQ, TREND, ISWFLG, MTRSV, 1366C 8 MXRDSV) 1367C 1368 CALL CDSTP (EPS, F, FA, HMAX, IMPL, IERROR, JACOBN, MATDIM, 1369 8 IWORK(IMXORD), IWORK(IMNT), IWORK(IMTR), ML, MU, NPAR, 1370 8 NDECOM, WORK(IYWT), UROUND, USERS, AVGH, AVGORD, H, 1371 8 HUSED, IWORK(IJTASK), IWORK(IMNTLD), IWORK(IMTRLD), 1372 8 IWORK(INFE), IWORK(INJE), IWORK(INQUSE), IWORK(INSTEP), 1373 8 T, Y, WORK(IYH), WORK(IA), CONVRG, WORK(IDFDY), EL, 1374 8 WORK(IFAC), HOLD, IWORK(INDPVT), JSTATE, IWORK(IJSTPL), 1375 8 IWORK(INQ), IWORK(INWAIT), RC, RMAX, WORK(ISAVE1), 1376 8 WORK(ISAVE2), TQ, TREND, MINT, IWORK(IMTRSV), 1377 8 IWORK(IMXRDS)) 1378C 1379 WORK(IH) = H 1380 WORK(IT) = T 1381 GO TO (470, 670, 680, 690, 690, 660, 660, 660, 660, 660), JSTATE 1382 470 IWORK(IJTASK) = 1 1383C Determine if a root has been overtaken 1384 IF (NROOT .NE. 0) THEN 1385 IROOT = 0 1386 DO 500 I = 1,NROOT 1387 GLAST = WORK(I+IGNOW-1) 1388 GNOW = G (NPAR, T, Y, I) 1389 IF (NPAR .EQ. 0) THEN 1390 IWORK(INROOT) = I 1391 NSTATE = 7 1392 RETURN 1393 END IF 1394 WORK(I+IGNOW-1) = GNOW 1395 IF (GLAST*GNOW .GT. 0.E0) THEN 1396 WORK(I+ITROOT-1) = T + H 1397 ELSE 1398 IF (GNOW .EQ. 0.E0) THEN 1399 WORK(I+ITROOT-1) = T 1400 IROOT = I 1401 ELSE 1402 IF (GLAST .EQ. 0.E0) THEN 1403 WORK(I+ITROOT-1) = T + H 1404 ELSE 1405 IF (ABS(HUSED) .GE. UROUND*ABS(T)) THEN 1406 TLAST = T - HUSED 1407 IROOT = I 1408 TROOT = T 1409 CALL CDZRO (AE, G, H, NPAR, IWORK(INQ), IROOT, RE, T, 1410 8 WORK(IYH), UROUND, TROOT, TLAST, 1411 8 GNOW, GLAST, Y) 1412 DO 480 J = 1,N 1413 480 Y(J) = WORK(IYH+J-1) 1414 IF (NPAR .EQ. 0) THEN 1415 IWORK(INROOT) = I 1416 NSTATE = 7 1417 RETURN 1418 END IF 1419 WORK(I+ITROOT-1) = TROOT 1420 ELSE 1421 WORK(I+ITROOT-1) = T 1422 IROOT = I 1423 END IF 1424 END IF 1425 END IF 1426 END IF 1427 500 CONTINUE 1428 IF (IROOT .EQ. 0) THEN 1429 IWORK(IJROOT) = 0 1430C Select the first root 1431 ELSE 1432 IWORK(IJROOT) = NTASK 1433 IWORK(INRTLD) = NROOT 1434 IWORK(INDTRT) = ITROOT 1435 TROOT = T + H 1436 DO 510 I = 1,NROOT 1437 IF (REAL(WORK(I+ITROOT-1))*HSIGN .LT. TROOT*HSIGN) THEN 1438 TROOT = WORK(I+ITROOT-1) 1439 IROOT = I 1440 END IF 1441 510 CONTINUE 1442 IWORK(INROOT) = IROOT 1443 WORK(ITOUT) = TROOT 1444 IF (TROOT*HSIGN .LE. TOUT*HSIGN) THEN 1445 CALL CDNTP (H, 0, N, IWORK(INQ), T, TROOT, WORK(IYH), Y) 1446 NSTATE = 5 1447 T = TROOT 1448 IERFLG = 0 1449 GO TO 580 1450 END IF 1451 END IF 1452 END IF 1453C Test for NTASK condition to be satisfied 1454 NSTATE = 2 1455 IF (NTASK .EQ. 1) THEN 1456 IF (T*HSIGN .LT. TOUT*HSIGN) GO TO 260 1457 CALL CDNTP (H, 0, N, IWORK(INQ), T, TOUT, WORK(IYH), Y) 1458 T = TOUT 1459 IERFLG = 0 1460 GO TO 580 1461C TOUT is assumed to have been attained 1462C exactly if T is within twenty roundoff 1463C units of TOUT, relative to MAX(TOUT, T). 1464C 1465 ELSE IF (NTASK .EQ. 2) THEN 1466 IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN 1467 T = TOUT 1468 ELSE 1469 IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN 1470 H = TOUT - T 1471 IF ((T + H)*HSIGN.GT.TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND) 1472 WORK(IH) = H 1473 IF (H .EQ. 0.E0) GO TO 670 1474 IWORK(IJTASK) = -1 1475 END IF 1476 END IF 1477 ELSE IF (NTASK .EQ. 3) THEN 1478 IF (ABS(TOUT - T).LE.NROUND*UROUND*MAX(ABS(T), ABS(TOUT))) THEN 1479 T = TOUT 1480 ELSE 1481 IF ((T + H)*HSIGN .GT. TOUT*HSIGN) THEN 1482 H = TOUT - T 1483 IF ((T + H)*HSIGN.GT.TOUT*HSIGN) H = H*(1.E0 - 4.E0*UROUND) 1484 WORK(IH) = H 1485 IF (H .EQ. 0.E0) GO TO 670 1486 IWORK(IJTASK) = -1 1487 END IF 1488 GO TO 260 1489 END IF 1490 END IF 1491 IERFLG = 0 1492C All returns are made through this 1493C section. IMXERR is determined. 1494 560 DO 570 I = 1,N 1495 570 Y(I) = WORK(I+IYH-1) 1496 580 IF (CONVRG) THEN 1497 IWORK(ICNVRG) = 1 1498 ELSE 1499 IWORK(ICNVRG) = 0 1500 END IF 1501 WORK(IAVGH) = AVGH 1502 WORK(IAVGRD) = AVGORD 1503 WORK(IHUSED) = HUSED 1504 WORK(IHOLD) = HOLD 1505 WORK(IRC) = RC 1506 WORK(IRMAX) = RMAX 1507 WORK(ITREND) = TREND 1508 DO 582 J = 1,12 1509 DO 582 I = 1,13 1510 582 WORK(I+IEL+(J-1)*13-1) = EL(I,J) 1511 DO 584 J = 1,12 1512 DO 584 I = 1,3 1513 584 WORK(I+ITQ+(J-1)*3-1) = TQ(I,J) 1514 IF (IWORK(IJTASK) .EQ. 0) RETURN 1515 BIG = 0.E0 1516 IMXERR = 1 1517 DO 590 I = 1,N 1518C SIZE = ABS(ERROR(I)/YWT(I)) 1519 SIZE = ABS(WORK(I+ISAVE1-1)/WORK(I+IYWT-1)) 1520 IF (BIG .LT. SIZE) THEN 1521 BIG = SIZE 1522 IMXERR = I 1523 END IF 1524 590 CONTINUE 1525 IWORK(INDMXR) = IMXERR 1526 RETURN 1527C 1528 660 NSTATE = JSTATE 1529 DO 662 I = 1,N 1530 662 Y(I) = WORK(I + IYH - 1) 1531 IF (CONVRG) THEN 1532 IWORK(ICNVRG) = 1 1533 ELSE 1534 IWORK(ICNVRG) = 0 1535 END IF 1536 WORK(IAVGH) = AVGH 1537 WORK(IAVGRD) = AVGORD 1538 WORK(IHUSED) = HUSED 1539 WORK(IHOLD) = HOLD 1540 WORK(IRC) = RC 1541 WORK(IRMAX) = RMAX 1542 WORK(ITREND) = TREND 1543 DO 664 J = 1,12 1544 DO 664 I = 1,13 1545 664 WORK(I+IEL+(J-1)*13-1) = EL(I,J) 1546 DO 666 J = 1,12 1547 DO 666 I = 1,3 1548 666 WORK(I+ITQ+(J-1)*3-1) = TQ(I,J) 1549 RETURN 1550C Fatal errors are processed here 1551C 1552 670 WRITE(RL1, '(E16.8)') T 1553 IERFLG = 41 1554 CALL XERMSG('SLATEC', 'CDRIV3', 1555 8 'At T, '//RL1//', the attempted step size has gone to '// 1556 8 'zero. Often this occurs if the problem setup is incorrect.', 1557 8 IERFLG, 1) 1558 NSTATE = 12 1559 RETURN 1560C 1561 680 WRITE(RL1, '(E16.8)') T 1562 IERFLG = 42 1563 CALL XERMSG('SLATEC', 'CDRIV3', 1564 8 'At T, '//RL1//', the step size has been reduced about 50 '// 1565 8 'times without advancing the solution. Often this occurs '// 1566 8 'if the problem setup is incorrect.', IERFLG, 1) 1567 NSTATE = 12 1568 RETURN 1569C 1570 690 WRITE(RL1, '(E16.8)') T 1571 IERFLG = 43 1572 CALL XERMSG('SLATEC', 'CDRIV3', 1573 8 'At T, '//RL1//', while solving A*YDOT = F, A is singular.', 1574 8 IERFLG, 1) 1575 NSTATE = 12 1576 RETURN 1577 END 1578