1! 2! SLATEC: public domain 3! 4*DECK DDEBDF 5 SUBROUTINE DDEBDF (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, 6 + RWORK, LRW, IWORK, LIW, RPAR, IPAR, DJAC) 7C***BEGIN PROLOGUE DDEBDF 8C***PURPOSE Solve an initial value problem in ordinary differential 9C equations using backward differentiation formulas. It is 10C intended primarily for stiff problems. 11C***LIBRARY SLATEC (DEPAC) 12C***CATEGORY I1A2 13C***TYPE DOUBLE PRECISION (DEBDF-S, DDEBDF-D) 14C***KEYWORDS BACKWARD DIFFERENTIATION FORMULAS, DEPAC, 15C INITIAL VALUE PROBLEMS, ODE, 16C ORDINARY DIFFERENTIAL EQUATIONS, STIFF 17C***AUTHOR Shampine, L. F., (SNLA) 18C Watts, H. A., (SNLA) 19C***DESCRIPTION 20C 21C This is the backward differentiation code in the package of 22C differential equation solvers DEPAC, consisting of the codes 23C DDERKF, DDEABM, and DDEBDF. Design of the package was by 24C L. F. Shampine and H. A. Watts. It is documented in 25C SAND-79-2374 , DEPAC - Design of a User Oriented Package of ODE 26C Solvers. 27C DDEBDF is a driver for a modification of the code LSODE written by 28C A. C. Hindmarsh 29C Lawrence Livermore Laboratory 30C Livermore, California 94550 31C 32C ********************************************************************** 33C ** DEPAC PACKAGE OVERVIEW ** 34C ********************************************************************** 35C 36C You have a choice of three differential equation solvers from 37C DEPAC. The following brief descriptions are meant to aid you 38C in choosing the most appropriate code for your problem. 39C 40C DDERKF is a fifth order Runge-Kutta code. It is the simplest of 41C the three choices, both algorithmically and in the use of the 42C code. DDERKF is primarily designed to solve non-stiff and mild- 43C ly stiff differential equations when derivative evaluations are 44C not expensive. It should generally not be used to get high 45C accuracy results nor answers at a great many specific points. 46C Because DDERKF has very low overhead costs, it will usually 47C result in the least expensive integration when solving 48C problems requiring a modest amount of accuracy and having 49C equations that are not costly to evaluate. DDERKF attempts to 50C discover when it is not suitable for the task posed. 51C 52C DDEABM is a variable order (one through twelve) Adams code. Its 53C complexity lies somewhere between that of DDERKF and DDEBDF. 54C DDEABM is primarily designed to solve non-stiff and mildly 55C stiff differential equations when derivative evaluations are 56C expensive, high accuracy results are needed or answers at 57C many specific points are required. DDEABM attempts to discover 58C when it is not suitable for the task posed. 59C 60C DDEBDF is a variable order (one through five) backward 61C differentiation formula code. It is the most complicated of 62C the three choices. DDEBDF is primarily designed to solve stiff 63C differential equations at crude to moderate tolerances. 64C If the problem is very stiff at all, DDERKF and DDEABM will be 65C quite inefficient compared to DDEBDF. However, DDEBDF will be 66C inefficient compared to DDERKF and DDEABM on non-stiff problems 67C because it uses much more storage, has a much larger overhead, 68C and the low order formulas will not give high accuracies 69C efficiently. 70C 71C The concept of stiffness cannot be described in a few words. 72C If you do not know the problem to be stiff, try either DDERKF 73C or DDEABM. Both of these codes will inform you of stiffness 74C when the cost of solving such problems becomes important. 75C 76C ********************************************************************** 77C ** ABSTRACT ** 78C ********************************************************************** 79C 80C Subroutine DDEBDF uses the backward differentiation formulas of 81C orders one through five to integrate a system of NEQ first order 82C ordinary differential equations of the form 83C DU/DX = DF(X,U) 84C when the vector Y(*) of initial values for U(*) at X=T is given. 85C The subroutine integrates from T to TOUT. It is easy to continue the 86C integration to get results at additional TOUT. This is the interval 87C mode of operation. It is also easy for the routine to return with 88C the solution at each intermediate step on the way to TOUT. This is 89C the intermediate-output mode of operation. 90C 91C ********************************************************************** 92C * Description of The Arguments To DDEBDF (An Overview) * 93C ********************************************************************** 94C 95C The Parameters are: 96C 97C DF -- This is the name of a subroutine which you provide to 98C define the differential equations. 99C 100C NEQ -- This is the number of (first order) differential 101C equations to be integrated. 102C 103C T -- This is a DOUBLE PRECISION value of the independent 104C variable. 105C 106C Y(*) -- This DOUBLE PRECISION array contains the solution 107C components at T. 108C 109C TOUT -- This is a DOUBLE PRECISION point at which a solution is 110C desired. 111C 112C INFO(*) -- The basic task of the code is to integrate the 113C differential equations from T to TOUT and return an 114C answer at TOUT. INFO(*) is an INTEGER array which is used 115C to communicate exactly how you want this task to be 116C carried out. 117C 118C RTOL, ATOL -- These DOUBLE PRECISION quantities 119C represent relative and absolute error tolerances which you 120C provide to indicate how accurately you wish the solution 121C to be computed. You may choose them to be both scalars 122C or else both vectors. 123C 124C IDID -- This scalar quantity is an indicator reporting what 125C the code did. You must monitor this INTEGER variable to 126C decide what action to take next. 127C 128C RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of 129C length LRW which provides the code with needed storage 130C space. 131C 132C IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW 133C which provides the code with needed storage space and an 134C across call flag. 135C 136C RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter 137C arrays which you can use for communication between your 138C calling program and the DF subroutine (and the DJAC 139C subroutine). 140C 141C DJAC -- This is the name of a subroutine which you may choose to 142C provide for defining the Jacobian matrix of partial 143C derivatives DF/DU. 144C 145C Quantities which are used as input items are 146C NEQ, T, Y(*), TOUT, INFO(*), 147C RTOL, ATOL, RWORK(1), LRW, 148C IWORK(1), IWORK(2), and LIW. 149C 150C Quantities which may be altered by the code are 151C T, Y(*), INFO(1), RTOL, ATOL, 152C IDID, RWORK(*) and IWORK(*). 153C 154C ********************************************************************** 155C * INPUT -- What To Do On The First Call To DDEBDF * 156C ********************************************************************** 157C 158C The first call of the code is defined to be the start of each new 159C problem. Read through the descriptions of all the following items, 160C provide sufficient storage space for designated arrays, set 161C appropriate variables for the initialization of the problem, and 162C give information about how you want the problem to be solved. 163C 164C 165C DF -- Provide a subroutine of the form 166C DF(X,U,UPRIME,RPAR,IPAR) 167C to define the system of first order differential equations 168C which is to be solved. For the given values of X and the 169C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must 170C evaluate the NEQ components of the system of differential 171C equations DU/DX=DF(X,U) and store the derivatives in the 172C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for 173C equations I=1,...,NEQ. 174C 175C Subroutine DF must not alter X or U(*). You must declare 176C the name DF in an external statement in your program that 177C calls DDEBDF. You must dimension U and UPRIME in DF. 178C 179C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter 180C arrays which you can use for communication between your 181C calling program and subroutine DF. They are not used or 182C altered by DDEBDF. If you do not need RPAR or IPAR, 183C ignore these parameters by treating them as dummy 184C arguments. If you do choose to use them, dimension them in 185C your calling program and in DF as arrays of appropriate 186C length. 187C 188C NEQ -- Set it to the number of differential equations. 189C (NEQ .GE. 1) 190C 191C T -- Set it to the initial point of the integration. 192C You must use a program variable for T because the code 193C changes its value. 194C 195C Y(*) -- Set this vector to the initial values of the NEQ solution 196C components at the initial point. You must dimension Y at 197C least NEQ in your calling program. 198C 199C TOUT -- Set it to the first point at which a solution is desired. 200C You can take TOUT = T, in which case the code 201C will evaluate the derivative of the solution at T and 202C return. Integration either forward in T (TOUT .GT. T) 203C or backward in T (TOUT .LT. T) is permitted. 204C 205C The code advances the solution from T to TOUT using 206C step sizes which are automatically selected so as to 207C achieve the desired accuracy. If you wish, the code will 208C return with the solution and its derivative following 209C each intermediate step (intermediate-output mode) so that 210C you can monitor them, but you still must provide TOUT in 211C accord with the basic aim of the code. 212C 213C The first step taken by the code is a critical one 214C because it must reflect how fast the solution changes near 215C the initial point. The code automatically selects an 216C initial step size which is practically always suitable for 217C the problem. By using the fact that the code will not 218C step past TOUT in the first step, you could, if necessary, 219C restrict the length of the initial step size. 220C 221C For some problems it may not be permissible to integrate 222C past a point TSTOP because a discontinuity occurs there 223C or the solution or its derivative is not defined beyond 224C TSTOP. When you have declared a TSTOP point (see INFO(4) 225C and RWORK(1)), you have told the code not to integrate 226C past TSTOP. In this case any TOUT beyond TSTOP is invalid 227C input. 228C 229C INFO(*) -- Use the INFO array to give the code more details about 230C how you want your problem solved. This array should be 231C dimensioned of length 15 to accommodate other members of 232C DEPAC or possible future extensions, though DDEBDF uses 233C only the first six entries. You must respond to all of 234C the following items which are arranged as questions. The 235C simplest use of the code corresponds to answering all 236C questions as YES ,i.e. setting all entries of INFO to 0. 237C 238C INFO(1) -- This parameter enables the code to initialize 239C itself. You must set it to indicate the start of every 240C new problem. 241C 242C **** Is this the first call for this problem ... 243C YES -- Set INFO(1) = 0 244C NO -- Not applicable here. 245C See below for continuation calls. **** 246C 247C INFO(2) -- How much accuracy you want of your solution 248C is specified by the error tolerances RTOL and ATOL. 249C The simplest use is to take them both to be scalars. 250C To obtain more flexibility, they can both be vectors. 251C The code must be told your choice. 252C 253C **** Are both error tolerances RTOL, ATOL scalars ... 254C YES -- Set INFO(2) = 0 255C and input scalars for both RTOL and ATOL 256C NO -- Set INFO(2) = 1 257C and input arrays for both RTOL and ATOL **** 258C 259C INFO(3) -- The code integrates from T in the direction 260C of TOUT by steps. If you wish, it will return the 261C computed solution and derivative at the next 262C intermediate step (the intermediate-output mode) or 263C TOUT, whichever comes first. This is a good way to 264C proceed if you want to see the behavior of the solution. 265C If you must have solutions at a great many specific 266C TOUT points, this code will compute them efficiently. 267C 268C **** Do you want the solution only at 269C TOUT (and NOT at the next intermediate step) ... 270C YES -- Set INFO(3) = 0 271C NO -- Set INFO(3) = 1 **** 272C 273C INFO(4) -- To handle solutions at a great many specific 274C values TOUT efficiently, this code may integrate past 275C TOUT and interpolate to obtain the result at TOUT. 276C Sometimes it is not possible to integrate beyond some 277C point TSTOP because the equation changes there or it is 278C not defined past TSTOP. Then you must tell the code 279C not to go past. 280C 281C **** Can the integration be carried out without any 282C restrictions on the independent variable T ... 283C YES -- Set INFO(4)=0 284C NO -- Set INFO(4)=1 285C and define the stopping point TSTOP by 286C setting RWORK(1)=TSTOP **** 287C 288C INFO(5) -- To solve stiff problems it is necessary to use the 289C Jacobian matrix of partial derivatives of the system 290C of differential equations. If you do not provide a 291C subroutine to evaluate it analytically (see the 292C description of the item DJAC in the call list), it will 293C be approximated by numerical differencing in this code. 294C Although it is less trouble for you to have the code 295C compute partial derivatives by numerical differencing, 296C the solution will be more reliable if you provide the 297C derivatives via DJAC. Sometimes numerical differencing 298C is cheaper than evaluating derivatives in DJAC and 299C sometimes it is not - this depends on your problem. 300C 301C If your problem is linear, i.e. has the form 302C DU/DX = DF(X,U) = J(X)*U + G(X) for some matrix J(X) 303C and vector G(X), the Jacobian matrix DF/DU = J(X). 304C Since you must provide a subroutine to evaluate DF(X,U) 305C analytically, it is little extra trouble to provide 306C subroutine DJAC for evaluating J(X) analytically. 307C Furthermore, in such cases, numerical differencing is 308C much more expensive than analytic evaluation. 309C 310C **** Do you want the code to evaluate the partial 311C derivatives automatically by numerical differences ... 312C YES -- Set INFO(5)=0 313C NO -- Set INFO(5)=1 314C and provide subroutine DJAC for evaluating the 315C Jacobian matrix **** 316C 317C INFO(6) -- DDEBDF will perform much better if the Jacobian 318C matrix is banded and the code is told this. In this 319C case, the storage needed will be greatly reduced, 320C numerical differencing will be performed more cheaply, 321C and a number of important algorithms will execute much 322C faster. The differential equation is said to have 323C half-bandwidths ML (lower) and MU (upper) if equation I 324C involves only unknowns Y(J) with 325C I-ML .LE. J .LE. I+MU 326C for all I=1,2,...,NEQ. Thus, ML and MU are the widths 327C of the lower and upper parts of the band, respectively, 328C with the main diagonal being excluded. If you do not 329C indicate that the equation has a banded Jacobian, 330C the code works with a full matrix of NEQ**2 elements 331C (stored in the conventional way). Computations with 332C banded matrices cost less time and storage than with 333C full matrices if 2*ML+MU .LT. NEQ. If you tell the 334C code that the Jacobian matrix has a banded structure and 335C you want to provide subroutine DJAC to compute the 336C partial derivatives, then you must be careful to store 337C the elements of the Jacobian matrix in the special form 338C indicated in the description of DJAC. 339C 340C **** Do you want to solve the problem using a full 341C (dense) Jacobian matrix (and not a special banded 342C structure) ... 343C YES -- Set INFO(6)=0 344C NO -- Set INFO(6)=1 345C and provide the lower (ML) and upper (MU) 346C bandwidths by setting 347C IWORK(1)=ML 348C IWORK(2)=MU **** 349C 350C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL) 351C error tolerances to tell the code how accurately you want 352C the solution to be computed. They must be defined as 353C program variables because the code may change them. You 354C have two choices -- 355C Both RTOL and ATOL are scalars. (INFO(2)=0) 356C Both RTOL and ATOL are vectors. (INFO(2)=1) 357C In either case all components must be non-negative. 358C 359C The tolerances are used by the code in a local error test 360C at each step which requires roughly that 361C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL 362C for each vector component. 363C (More specifically, a root-mean-square norm is used to 364C measure the size of vectors, and the error test uses the 365C magnitude of the solution at the beginning of the step.) 366C 367C The true (global) error is the difference between the true 368C solution of the initial value problem and the computed 369C approximation. Practically all present day codes, 370C including this one, control the local error at each step 371C and do not even attempt to control the global error 372C directly. Roughly speaking, they produce a solution Y(T) 373C which satisfies the differential equations with a 374C residual R(T), DY(T)/DT = DF(T,Y(T)) + R(T) , 375C and, almost always, R(T) is bounded by the error 376C tolerances. Usually, but not always, the true accuracy of 377C the computed Y is comparable to the error tolerances. This 378C code will usually, but not always, deliver a more accurate 379C solution if you reduce the tolerances and integrate again. 380C By comparing two such solutions you can get a fairly 381C reliable idea of the true error in the solution at the 382C bigger tolerances. 383C 384C Setting ATOL=0. results in a pure relative error test on 385C that component. Setting RTOL=0. results in a pure abso- 386C lute error test on that component. A mixed test with non- 387C zero RTOL and ATOL corresponds roughly to a relative error 388C test when the solution component is much bigger than ATOL 389C and to an absolute error test when the solution component 390C is smaller than the threshold ATOL. 391C 392C Proper selection of the absolute error control parameters 393C ATOL requires you to have some idea of the scale of the 394C solution components. To acquire this information may mean 395C that you will have to solve the problem more than once. In 396C the absence of scale information, you should ask for some 397C relative accuracy in all the components (by setting RTOL 398C values non-zero) and perhaps impose extremely small 399C absolute error tolerances to protect against the danger of 400C a solution component becoming zero. 401C 402C The code will not attempt to compute a solution at an 403C accuracy unreasonable for the machine being used. It will 404C advise you if you ask for too much accuracy and inform 405C you as to the maximum accuracy it believes possible. 406C 407C RWORK(*) -- Dimension this DOUBLE PRECISION work array of length 408C LRW in your calling program. 409C 410C RWORK(1) -- If you have set INFO(4)=0, you can ignore this 411C optional input parameter. Otherwise you must define a 412C stopping point TSTOP by setting RWORK(1) = TSTOP. 413C (For some problems it may not be permissible to integrate 414C past a point TSTOP because a discontinuity occurs there 415C or the solution or its derivative is not defined beyond 416C TSTOP.) 417C 418C LRW -- Set it to the declared length of the RWORK array. 419C You must have 420C LRW .GE. 250+10*NEQ+NEQ**2 421C for the full (dense) Jacobian case (when INFO(6)=0), or 422C LRW .GE. 250+10*NEQ+(2*ML+MU+1)*NEQ 423C for the banded Jacobian case (when INFO(6)=1). 424C 425C IWORK(*) -- Dimension this INTEGER work array of length LIW in 426C your calling program. 427C 428C IWORK(1), IWORK(2) -- If you have set INFO(6)=0, you can ignore 429C these optional input parameters. Otherwise you must define 430C the half-bandwidths ML (lower) and MU (upper) of the 431C Jacobian matrix by setting IWORK(1) = ML and 432C IWORK(2) = MU. (The code will work with a full matrix 433C of NEQ**2 elements unless it is told that the problem has 434C a banded Jacobian, in which case the code will work with 435C a matrix containing at most (2*ML+MU+1)*NEQ elements.) 436C 437C LIW -- Set it to the declared length of the IWORK array. 438C You must have LIW .GE. 56+NEQ. 439C 440C RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and 441C INTEGER type, respectively. You can use them for 442C communication between your program that calls DDEBDF and 443C the DF subroutine (and the DJAC subroutine). They are not 444C used or altered by DDEBDF. If you do not need RPAR or 445C IPAR, ignore these parameters by treating them as dummy 446C arguments. If you do choose to use them, dimension them in 447C your calling program and in DF (and in DJAC) as arrays of 448C appropriate length. 449C 450C DJAC -- If you have set INFO(5)=0, you can ignore this parameter 451C by treating it as a dummy argument. (For some compilers 452C you may have to write a dummy subroutine named DJAC in 453C order to avoid problems associated with missing external 454C routine names.) Otherwise, you must provide a subroutine 455C of the form 456C DJAC(X,U,PD,NROWPD,RPAR,IPAR) 457C to define the Jacobian matrix of partial derivatives DF/DU 458C of the system of differential equations DU/DX = DF(X,U). 459C For the given values of X and the vector 460C U(*)=(U(1),U(2),...,U(NEQ)), the subroutine must evaluate 461C the non-zero partial derivatives DF(I)/DU(J) for each 462C differential equation I=1,...,NEQ and each solution 463C component J=1,...,NEQ , and store these values in the 464C matrix PD. The elements of PD are set to zero before each 465C call to DJAC so only non-zero elements need to be defined. 466C 467C Subroutine DJAC must not alter X, U(*), or NROWPD. You 468C must declare the name DJAC in an external statement in 469C your program that calls DDEBDF. NROWPD is the row 470C dimension of the PD matrix and is assigned by the code. 471C Therefore you must dimension PD in DJAC according to 472C DIMENSION PD(NROWPD,1) 473C You must also dimension U in DJAC. 474C 475C The way you must store the elements into the PD matrix 476C depends on the structure of the Jacobian which you 477C indicated by INFO(6). 478C *** INFO(6)=0 -- Full (Dense) Jacobian *** 479C When you evaluate the (non-zero) partial derivative 480C of equation I with respect to variable J, you must 481C store it in PD according to 482C PD(I,J) = * DF(I)/DU(J) * 483C *** INFO(6)=1 -- Banded Jacobian with ML Lower and MU 484C Upper Diagonal Bands (refer to INFO(6) description of 485C ML and MU) *** 486C When you evaluate the (non-zero) partial derivative 487C of equation I with respect to variable J, you must 488C store it in PD according to 489C IROW = I - J + ML + MU + 1 490C PD(IROW,J) = * DF(I)/DU(J) * 491C 492C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter 493C arrays which you can use for communication between your 494C calling program and your Jacobian subroutine DJAC. They 495C are not altered by DDEBDF. If you do not need RPAR or 496C IPAR, ignore these parameters by treating them as dummy 497C arguments. If you do choose to use them, dimension them 498C in your calling program and in DJAC as arrays of 499C appropriate length. 500C 501C ********************************************************************** 502C * OUTPUT -- After any return from DDEBDF * 503C ********************************************************************** 504C 505C The principal aim of the code is to return a computed solution at 506C TOUT, although it is also possible to obtain intermediate results 507C along the way. To find out whether the code achieved its goal 508C or if the integration process was interrupted before the task was 509C completed, you must check the IDID parameter. 510C 511C 512C T -- The solution was successfully advanced to the 513C output value of T. 514C 515C Y(*) -- Contains the computed solution approximation at T. 516C You may also be interested in the approximate derivative 517C of the solution at T. It is contained in 518C RWORK(21),...,RWORK(20+NEQ). 519C 520C IDID -- Reports what the code did 521C 522C *** Task Completed *** 523C Reported by positive values of IDID 524C 525C IDID = 1 -- A step was successfully taken in the 526C intermediate-output mode. The code has not 527C yet reached TOUT. 528C 529C IDID = 2 -- The integration to TOUT was successfully 530C completed (T=TOUT) by stepping exactly to TOUT. 531C 532C IDID = 3 -- The integration to TOUT was successfully 533C completed (T=TOUT) by stepping past TOUT. 534C Y(*) is obtained by interpolation. 535C 536C *** Task Interrupted *** 537C Reported by negative values of IDID 538C 539C IDID = -1 -- A large amount of work has been expended. 540C (500 steps attempted) 541C 542C IDID = -2 -- The error tolerances are too stringent. 543C 544C IDID = -3 -- The local error test cannot be satisfied 545C because you specified a zero component in ATOL 546C and the corresponding computed solution 547C component is zero. Thus, a pure relative error 548C test is impossible for this component. 549C 550C IDID = -4,-5 -- Not applicable for this code but used 551C by other members of DEPAC. 552C 553C IDID = -6 -- DDEBDF had repeated convergence test failures 554C on the last attempted step. 555C 556C IDID = -7 -- DDEBDF had repeated error test failures on 557C the last attempted step. 558C 559C IDID = -8,..,-32 -- Not applicable for this code but 560C used by other members of DEPAC or possible 561C future extensions. 562C 563C *** Task Terminated *** 564C Reported by the value of IDID=-33 565C 566C IDID = -33 -- The code has encountered trouble from which 567C it cannot recover. A message is printed 568C explaining the trouble and control is returned 569C to the calling program. For example, this 570C occurs when invalid input is detected. 571C 572C RTOL, ATOL -- These quantities remain unchanged except when 573C IDID = -2. In this case, the error tolerances have been 574C increased by the code to values which are estimated to be 575C appropriate for continuing the integration. However, the 576C reported solution at T was obtained using the input values 577C of RTOL and ATOL. 578C 579C RWORK, IWORK -- Contain information which is usually of no 580C interest to the user but necessary for subsequent calls. 581C However, you may find use for 582C 583C RWORK(11)--which contains the step size H to be 584C attempted on the next step. 585C 586C RWORK(12)--If the tolerances have been increased by the 587C code (IDID = -2) , they were multiplied by the 588C value in RWORK(12). 589C 590C RWORK(13)--which contains the current value of the 591C independent variable, i.e. the farthest point 592C integration has reached. This will be 593C different from T only when interpolation has 594C been performed (IDID=3). 595C 596C RWORK(20+I)--which contains the approximate derivative 597C of the solution component Y(I). In DDEBDF, it 598C is never obtained by calling subroutine DF to 599C evaluate the differential equation using T and 600C Y(*), except at the initial point of 601C integration. 602C 603C ********************************************************************** 604C ** INPUT -- What To Do To Continue The Integration ** 605C ** (calls after the first) ** 606C ********************************************************************** 607C 608C This code is organized so that subsequent calls to continue the 609C integration involve little (if any) additional effort on your 610C part. You must monitor the IDID parameter in order to determine 611C what to do next. 612C 613C Recalling that the principal task of the code is to integrate 614C from T to TOUT (the interval mode), usually all you will need 615C to do is specify a new TOUT upon reaching the current TOUT. 616C 617C Do not alter any quantity not specifically permitted below, 618C in particular do not alter NEQ, T, Y(*), RWORK(*), IWORK(*) or 619C the differential equation in subroutine DF. Any such alteration 620C constitutes a new problem and must be treated as such, i.e. 621C you must start afresh. 622C 623C You cannot change from vector to scalar error control or vice 624C versa (INFO(2)) but you can change the size of the entries of 625C RTOL, ATOL. Increasing a tolerance makes the equation easier 626C to integrate. Decreasing a tolerance will make the equation 627C harder to integrate and should generally be avoided. 628C 629C You can switch from the intermediate-output mode to the 630C interval mode (INFO(3)) or vice versa at any time. 631C 632C If it has been necessary to prevent the integration from going 633C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the 634C code will not integrate to any TOUT beyond the currently 635C specified TSTOP. Once TSTOP has been reached you must change 636C the value of TSTOP or set INFO(4)=0. You may change INFO(4) 637C or TSTOP at any time but you must supply the value of TSTOP in 638C RWORK(1) whenever you set INFO(4)=1. 639C 640C Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2) 641C unless you are going to restart the code. 642C 643C The parameter INFO(1) is used by the code to indicate the 644C beginning of a new problem and to indicate whether integration 645C is to be continued. You must input the value INFO(1) = 0 646C when starting a new problem. You must input the value 647C INFO(1) = 1 if you wish to continue after an interrupted task. 648C Do not set INFO(1) = 0 on a continuation call unless you 649C want the code to restart at the current T. 650C 651C *** Following a Completed Task *** 652C If 653C IDID = 1, call the code again to continue the integration 654C another step in the direction of TOUT. 655C 656C IDID = 2 or 3, define a new TOUT and call the code again. 657C TOUT must be different from T. You cannot change 658C the direction of integration without restarting. 659C 660C *** Following an Interrupted Task *** 661C To show the code that you realize the task was 662C interrupted and that you want to continue, you 663C must take appropriate action and reset INFO(1) = 1 664C If 665C IDID = -1, the code has attempted 500 steps. 666C If you want to continue, set INFO(1) = 1 and 667C call the code again. An additional 500 steps 668C will be allowed. 669C 670C IDID = -2, the error tolerances RTOL, ATOL have been 671C increased to values the code estimates appropriate 672C for continuing. You may want to change them 673C yourself. If you are sure you want to continue 674C with relaxed error tolerances, set INFO(1)=1 and 675C call the code again. 676C 677C IDID = -3, a solution component is zero and you set the 678C corresponding component of ATOL to zero. If you 679C are sure you want to continue, you must first 680C alter the error criterion to use positive values 681C for those components of ATOL corresponding to zero 682C solution components, then set INFO(1)=1 and call 683C the code again. 684C 685C IDID = -4,-5 --- cannot occur with this code but used 686C by other members of DEPAC. 687C 688C IDID = -6, repeated convergence test failures occurred 689C on the last attempted step in DDEBDF. An inaccu- 690C rate Jacobian may be the problem. If you are 691C absolutely certain you want to continue, restart 692C the integration at the current T by setting 693C INFO(1)=0 and call the code again. 694C 695C IDID = -7, repeated error test failures occurred on the 696C last attempted step in DDEBDF. A singularity in 697C the solution may be present. You should re- 698C examine the problem being solved. If you are 699C absolutely certain you want to continue, restart 700C the integration at the current T by setting 701C INFO(1)=0 and call the code again. 702C 703C IDID = -8,..,-32 --- cannot occur with this code but 704C used by other members of DDEPAC or possible future 705C extensions. 706C 707C *** Following a Terminated Task *** 708C If 709C IDID = -33, you cannot continue the solution of this 710C problem. An attempt to do so will result in your 711C run being terminated. 712C 713C ********************************************************************** 714C 715C ***** Warning ***** 716C 717C If DDEBDF is to be used in an overlay situation, you must save and 718C restore certain items used internally by DDEBDF (values in the 719C common block DDEBD1). This can be accomplished as follows. 720C 721C To save the necessary values upon return from DDEBDF, simply call 722C DSVCO(RWORK(22+NEQ),IWORK(21+NEQ)). 723C 724C To restore the necessary values before the next call to DDEBDF, 725C simply call DRSCO(RWORK(22+NEQ),IWORK(21+NEQ)). 726C 727C***REFERENCES L. F. Shampine and H. A. Watts, DEPAC - design of a user 728C oriented package of ODE solvers, Report SAND79-2374, 729C Sandia Laboratories, 1979. 730C***ROUTINES CALLED DLSOD, XERMSG 731C***COMMON BLOCKS DDEBD1 732C***REVISION HISTORY (YYMMDD) 733C 820301 DATE WRITTEN 734C 890831 Modified array declarations. (WRB) 735C 891024 Changed references from DVNORM to DHVNRM. (WRB) 736C 891024 REVISION DATE from Version 3.2 737C 891214 Prologue converted to Version 4.0 format. (BAB) 738C 900326 Removed duplicate information from DESCRIPTION section. 739C (WRB) 740C 900510 Convert XERRWV calls to XERMSG calls, make Prologue comments 741C consistent with DEBDF. (RWC) 742C 920501 Reformatted the REFERENCES section. (WRB) 743C***END PROLOGUE DDEBDF 744 INTEGER IACOR, IBAND, IBEGIN, ICOMI, ICOMR, IDELSN, IDID, IER, 745 1 IEWT, IINOUT, IINTEG, IJAC, ILRW, INFO, INIT, 746 2 IOWNS, IPAR, IQUIT, ISAVF, ITOL, ITSTAR, ITSTOP, IWM, 747 3 IWORK, IYH, IYPOUT, JSTART, KFLAG, KSTEPS, L, LIW, LRW, 748 4 MAXORD, METH, MITER, ML, MU, N, NEQ, NFE, NJE, NQ, NQU, 749 5 NST 750 DOUBLE PRECISION ATOL, EL0, H, HMIN, HMXI, HU, ROWNS, RPAR, 751 1 RTOL, RWORK, T, TN, TOLD, TOUT, UROUND, Y 752 LOGICAL INTOUT 753 CHARACTER*8 XERN1, XERN2 754 CHARACTER*16 XERN3 755C 756 DIMENSION Y(*),INFO(15),RTOL(*),ATOL(*),RWORK(*),IWORK(*), 757 1 RPAR(*),IPAR(*) 758C 759 COMMON /DDEBD1/ TOLD,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND, 760 1 IQUIT,INIT,IYH,IEWT,IACOR,ISAVF,IWM,KSTEPS,IBEGIN, 761 2 ITOL,IINTEG,ITSTOP,IJAC,IBAND,IOWNS(6),IER,JSTART, 762 3 KFLAG,L,METH,MITER,MAXORD,N,NQ,NST,NFE,NJE,NQU 763C 764 EXTERNAL DF, DJAC 765C 766C CHECK FOR AN APPARENT INFINITE LOOP 767C 768C***FIRST EXECUTABLE STATEMENT DDEBDF 769 IF (INFO(1) .EQ. 0) IWORK(LIW) = 0 770C 771 IF (IWORK(LIW).GE. 5) THEN 772 IF (T .EQ. RWORK(21+NEQ)) THEN 773 WRITE (XERN3, '(1PE15.6)') T 774 CALL XERMSG ('SLATEC', 'DDEBDF', 775 * 'AN APPARENT INFINITE LOOP HAS BEEN DETECTED.$$' // 776 * 'YOU HAVE MADE REPEATED CALLS AT T = ' // XERN3 // 777 * ' AND THE INTEGRATION HAS NOT ADVANCED. CHECK THE ' // 778 * 'WAY YOU HAVE SET PARAMETERS FOR THE CALL TO THE ' // 779 * 'CODE, PARTICULARLY INFO(1).', 13, 2) 780 RETURN 781 ENDIF 782 ENDIF 783C 784 IDID = 0 785C 786C CHECK VALIDITY OF INFO PARAMETERS 787C 788 IF (INFO(1) .NE. 0 .AND. INFO(1) .NE. 1) THEN 789 WRITE (XERN1, '(I8)') INFO(1) 790 CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(1) MUST BE SET TO 0 ' // 791 * 'FOR THE START OF A NEW PROBLEM, AND MUST BE SET TO 1 ' // 792 * 'FOLLOWING AN INTERRUPTED TASK. YOU ARE ATTEMPTING TO ' // 793 * 'CONTINUE THE INTEGRATION ILLEGALLY BY CALLING THE ' // 794 * 'CODE WITH INFO(1) = ' // XERN1, 3, 1) 795 IDID = -33 796 ENDIF 797C 798 IF (INFO(2) .NE. 0 .AND. INFO(2) .NE. 1) THEN 799 WRITE (XERN1, '(I8)') INFO(2) 800 CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(2) MUST BE 0 OR 1 ' // 801 * 'INDICATING SCALAR AND VECTOR ERROR TOLERANCES, ' // 802 * 'RESPECTIVELY. YOU HAVE CALLED THE CODE WITH INFO(2) = ' // 803 * XERN1, 4, 1) 804 IDID = -33 805 ENDIF 806C 807 IF (INFO(3) .NE. 0 .AND. INFO(3) .NE. 1) THEN 808 WRITE (XERN1, '(I8)') INFO(3) 809 CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(3) MUST BE 0 OR 1 ' // 810 * 'INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT MODE OF ' // 811 * 'INTEGRATION, RESPECTIVELY. YOU HAVE CALLED THE CODE ' // 812 * 'WITH INFO(3) = ' // XERN1, 5, 1) 813 IDID = -33 814 ENDIF 815C 816 IF (INFO(4) .NE. 0 .AND. INFO(4) .NE. 1) THEN 817 WRITE (XERN1, '(I8)') INFO(4) 818 CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(4) MUST BE 0 OR 1 ' // 819 * 'INDICATING WHETHER OR NOT THE INTEGRATION INTERVAL IS ' // 820 * 'TO BE RESTRICTED BY A POINT TSTOP. YOU HAVE CALLED ' // 821 * 'THE CODE WITH INFO(4) = ' // XERN1, 14, 1) 822 IDID = -33 823 ENDIF 824C 825 IF (INFO(5) .NE. 0 .AND. INFO(5) .NE. 1) THEN 826 WRITE (XERN1, '(I8)') INFO(5) 827 CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(5) MUST BE 0 OR 1 ' // 828 * 'INDICATING WHETHER THE CODE IS TOLD TO FORM THE ' // 829 * 'JACOBIAN MATRIX BY NUMERICAL DIFFERENCING OR YOU ' // 830 * 'PROVIDE A SUBROUTINE TO EVALUATE IT ANALYTICALLY. ' // 831 * 'YOU HAVE CALLED THE CODE WITH INFO(5) = ' // XERN1, 15, 1) 832 IDID = -33 833 ENDIF 834C 835 IF (INFO(6) .NE. 0 .AND. INFO(6) .NE. 1) THEN 836 WRITE (XERN1, '(I8)') INFO(6) 837 CALL XERMSG ('SLATEC', 'DDEBDF', 'INFO(6) MUST BE 0 OR 1 ' // 838 * 'INDICATING WHETHER THE CODE IS TOLD TO TREAT THE ' // 839 * 'JACOBIAN AS A FULL (DENSE) MATRIX OR AS HAVING A ' // 840 * 'SPECIAL BANDED STRUCTURE. YOU HAVE CALLED THE CODE ' // 841 * 'WITH INFO(6) = ' // XERN1, 16, 1) 842 IDID = -33 843 ENDIF 844C 845 ILRW = NEQ 846 IF (INFO(6) .NE. 0) THEN 847C 848C CHECK BANDWIDTH PARAMETERS 849C 850 ML = IWORK(1) 851 MU = IWORK(2) 852 ILRW = 2*ML + MU + 1 853C 854 IF (ML.LT.0 .OR. ML.GE.NEQ .OR. MU.LT.0 .OR. MU.GE.NEQ) THEN 855 WRITE (XERN1, '(I8)') ML 856 WRITE (XERN2, '(I8)') MU 857 CALL XERMSG ('SLATEC', 'DDEBDF', 'YOU HAVE SET INFO(6) ' // 858 * '= 1, TELLING THE CODE THAT THE JACOBIAN MATRIX HAS ' // 859 * 'A SPECIAL BANDED STRUCTURE. HOWEVER, THE LOWER ' // 860 * '(UPPER) BANDWIDTHS ML (MU) VIOLATE THE CONSTRAINTS ' // 861 * 'ML,MU .GE. 0 AND ML,MU .LT. NEQ. YOU HAVE CALLED ' // 862 * 'THE CODE WITH ML = ' // XERN1 // ' AND MU = ' // XERN2, 863 * 17, 1) 864 IDID = -33 865 ENDIF 866 ENDIF 867C 868C CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION 869C 870 IF (LRW .LT. 250 + (10 + ILRW)*NEQ) THEN 871 WRITE (XERN1, '(I8)') LRW 872 IF (INFO(6) .EQ. 0) THEN 873 CALL XERMSG ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY RWORK ' // 874 * 'MUST BE AT LEAST 250 + 10*NEQ + NEQ*NEQ.$$' // 875 * 'YOU HAVE CALLED THE CODE WITH LRW = ' // XERN1, 1, 1) 876 ELSE 877 CALL XERMSG ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY RWORK ' // 878 * 'MUST BE AT LEAST 250 + 10*NEQ + (2*ML+MU+1)*NEQ.$$' // 879 * 'YOU HAVE CALLED THE CODE WITH LRW = ' // XERN1, 18, 1) 880 ENDIF 881 IDID = -33 882 ENDIF 883C 884 IF (LIW .LT. 56 + NEQ) THEN 885 WRITE (XERN1, '(I8)') LIW 886 CALL XERMSG ('SLATEC', 'DDEBDF', 'LENGTH OF ARRAY IWORK ' // 887 * 'BE AT LEAST 56 + NEQ. YOU HAVE CALLED THE CODE WITH ' // 888 * 'LIW = ' // XERN1, 2, 1) 889 IDID = -33 890 ENDIF 891C 892C COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK 893C ARRAY AND RESTORE COMMON BLOCK DATA 894C 895 ICOMI = 21 + NEQ 896 IINOUT = ICOMI + 33 897C 898 IYPOUT = 21 899 ITSTAR = 21 + NEQ 900 ICOMR = 22 + NEQ 901C 902 IF (INFO(1) .NE. 0) INTOUT = IWORK(IINOUT) .NE. (-1) 903C CALL DRSCO(RWORK(ICOMR),IWORK(ICOMI)) 904C 905 IYH = ICOMR + 218 906 IEWT = IYH + 6*NEQ 907 ISAVF = IEWT + NEQ 908 IACOR = ISAVF + NEQ 909 IWM = IACOR + NEQ 910 IDELSN = IWM + 2 + ILRW*NEQ 911C 912 IBEGIN = INFO(1) 913 ITOL = INFO(2) 914 IINTEG = INFO(3) 915 ITSTOP = INFO(4) 916 IJAC = INFO(5) 917 IBAND = INFO(6) 918 RWORK(ITSTAR) = T 919C 920 CALL DLSOD(DF,NEQ,T,Y,TOUT,RTOL,ATOL,IDID,RWORK(IYPOUT), 921 1 RWORK(IYH),RWORK(IYH),RWORK(IEWT),RWORK(ISAVF), 922 2 RWORK(IACOR),RWORK(IWM),IWORK(1),DJAC,INTOUT, 923 3 RWORK(1),RWORK(12),RWORK(IDELSN),RPAR,IPAR) 924C 925 IWORK(IINOUT) = -1 926 IF (INTOUT) IWORK(IINOUT) = 1 927C 928 IF (IDID .NE. (-2)) IWORK(LIW) = IWORK(LIW) + 1 929 IF (T .NE. RWORK(ITSTAR)) IWORK(LIW) = 0 930C CALL DSVCO(RWORK(ICOMR),IWORK(ICOMI)) 931 RWORK(11) = H 932 RWORK(13) = TN 933 INFO(1) = IBEGIN 934C 935 RETURN 936 END 937*DECK DLSOD 938 SUBROUTINE DLSOD (DF, NEQ, T, Y, TOUT, RTOL, ATOL, IDID, YPOUT, 939 + YH, YH1, EWT, SAVF, ACOR, WM, IWM, DJAC, INTOUT, TSTOP, TOLFAC, 940 + DELSGN, RPAR, IPAR) 941C***BEGIN PROLOGUE DLSOD 942C***SUBSIDIARY 943C***PURPOSE Subsidiary to DDEBDF 944C***LIBRARY SLATEC 945C***TYPE DOUBLE PRECISION (LSOD-S, DLSOD-D) 946C***AUTHOR (UNKNOWN) 947C***DESCRIPTION 948C 949C DDEBDF merely allocates storage for DLSOD to relieve the user of 950C the inconvenience of a long call list. Consequently DLSOD is used 951C as described in the comments for DDEBDF . 952C 953C***SEE ALSO DDEBDF 954C***ROUTINES CALLED D1MACH, DHSTRT, DINTYD, DSTOD, DVNRMS, XERMSG 955C***COMMON BLOCKS DDEBD1 956C***REVISION HISTORY (YYMMDD) 957C 820301 DATE WRITTEN 958C 890531 Changed all specific intrinsics to generic. (WRB) 959C 890831 Modified array declarations. (WRB) 960C 891214 Prologue converted to Version 4.0 format. (BAB) 961C 900328 Added TYPE section. (WRB) 962C 900510 Convert XERRWV calls to XERMSG calls. (RWC) 963C***END PROLOGUE DLSOD 964C 965 INTEGER IBAND, IBEGIN, IDID, IER, IINTEG, IJAC, INIT, INTFLG, 966 1 IOWNS, IPAR, IQUIT, ITOL, ITSTOP, IWM, JSTART, K, KFLAG, 967 2 KSTEPS, L, LACOR, LDUM, LEWT, LSAVF, LTOL, LWM, LYH, MAXNUM, 968 3 MAXORD, METH, MITER, N, NATOLP, NEQ, NFE, NJE, NQ, NQU, 969 4 NRTOLP, NST 970 DOUBLE PRECISION ABSDEL, ACOR, ATOL, BIG, D1MACH, DEL, 971 1 DELSGN, DT, DVNRMS, EL0, EWT, 972 2 H, HA, HMIN, HMXI, HU, ROWNS, RPAR, RTOL, SAVF, T, TOL, 973 3 TOLD, TOLFAC, TOUT, TSTOP, U, WM, X, Y, YH, YH1, YPOUT 974 LOGICAL INTOUT 975 CHARACTER*8 XERN1 976 CHARACTER*16 XERN3, XERN4 977C 978 DIMENSION Y(*),YPOUT(*),YH(NEQ,6),YH1(*),EWT(*),SAVF(*), 979 1 ACOR(*),WM(*),IWM(*),RTOL(*),ATOL(*),RPAR(*),IPAR(*) 980C 981C 982 COMMON /DDEBD1/ TOLD,ROWNS(210),EL0,H,HMIN,HMXI,HU,X,U,IQUIT,INIT, 983 1 LYH,LEWT,LACOR,LSAVF,LWM,KSTEPS,IBEGIN,ITOL, 984 2 IINTEG,ITSTOP,IJAC,IBAND,IOWNS(6),IER,JSTART, 985 3 KFLAG,LDUM,METH,MITER,MAXORD,N,NQ,NST,NFE,NJE,NQU 986C 987 EXTERNAL DF, DJAC 988C 989C .................................................................. 990C 991C THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE 992C NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MAXNUM, THE 993C COUNTER IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE 994C EXCESSIVE WORK. 995 SAVE MAXNUM 996C 997 DATA MAXNUM /500/ 998C 999C .................................................................. 1000C 1001C***FIRST EXECUTABLE STATEMENT DLSOD 1002 IF (IBEGIN .EQ. 0) THEN 1003C 1004C ON THE FIRST CALL , PERFORM INITIALIZATION -- 1005C DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING THE 1006C FUNCTION ROUTINE D1MACH. THE USER MUST MAKE SURE THAT THE 1007C VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED. 1008C 1009 U = D1MACH(4) 1010C -- SET ASSOCIATED MACHINE DEPENDENT PARAMETER 1011 WM(1) = SQRT(U) 1012C -- SET TERMINATION FLAG 1013 IQUIT = 0 1014C -- SET INITIALIZATION INDICATOR 1015 INIT = 0 1016C -- SET COUNTER FOR ATTEMPTED STEPS 1017 KSTEPS = 0 1018C -- SET INDICATOR FOR INTERMEDIATE-OUTPUT 1019 INTOUT = .FALSE. 1020C -- SET START INDICATOR FOR DSTOD CODE 1021 JSTART = 0 1022C -- SET BDF METHOD INDICATOR 1023 METH = 2 1024C -- SET MAXIMUM ORDER FOR BDF METHOD 1025 MAXORD = 5 1026C -- SET ITERATION MATRIX INDICATOR 1027C 1028 IF (IJAC .EQ. 0 .AND. IBAND .EQ. 0) MITER = 2 1029 IF (IJAC .EQ. 1 .AND. IBAND .EQ. 0) MITER = 1 1030 IF (IJAC .EQ. 0 .AND. IBAND .EQ. 1) MITER = 5 1031 IF (IJAC .EQ. 1 .AND. IBAND .EQ. 1) MITER = 4 1032C 1033C -- SET OTHER NECESSARY ITEMS IN COMMON BLOCK 1034 N = NEQ 1035 NST = 0 1036 NJE = 0 1037 HMXI = 0.0D0 1038 NQ = 1 1039 H = 1.0D0 1040C -- RESET IBEGIN FOR SUBSEQUENT CALLS 1041 IBEGIN = 1 1042 ENDIF 1043C 1044C .................................................................. 1045C 1046C CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY 1047C 1048 IF (NEQ .LT. 1) THEN 1049 WRITE (XERN1, '(I8)') NEQ 1050 CALL XERMSG ('SLATEC', 'DLSOD', 1051 * 'IN DDEBDF, THE NUMBER OF EQUATIONS MUST BE A ' // 1052 * 'POSITIVE INTEGER.$$YOU HAVE CALLED THE CODE WITH NEQ = ' // 1053 * XERN1, 6, 1) 1054 IDID=-33 1055 ENDIF 1056C 1057 NRTOLP = 0 1058 NATOLP = 0 1059 DO 60 K = 1, NEQ 1060 IF (NRTOLP .LE. 0) THEN 1061 IF (RTOL(K) .LT. 0.) THEN 1062 WRITE (XERN1, '(I8)') K 1063 WRITE (XERN3, '(1PE15.6)') RTOL(K) 1064 CALL XERMSG ('SLATEC', 'DLSOD', 1065 * 'IN DDEBDF, THE RELATIVE ERROR TOLERANCES MUST ' // 1066 * 'BE NON-NEGATIVE.$$YOU HAVE CALLED THE CODE WITH ' // 1067 * 'RTOL(' // XERN1 // ') = ' // XERN3 // '$$IN THE ' // 1068 * 'CASE OF VECTOR ERROR TOLERANCES, NO FURTHER ' // 1069 * 'CHECKING OF RTOL COMPONENTS IS DONE.', 7, 1) 1070 IDID = -33 1071 IF (NATOLP .GT. 0) GO TO 70 1072 NRTOLP = 1 1073 ELSEIF (NATOLP .GT. 0) THEN 1074 GO TO 50 1075 ENDIF 1076 ENDIF 1077C 1078 IF (ATOL(K) .LT. 0.) THEN 1079 WRITE (XERN1, '(I8)') K 1080 WRITE (XERN3, '(1PE15.6)') ATOL(K) 1081 CALL XERMSG ('SLATEC', 'DLSOD', 1082 * 'IN DDEBDF, THE ABSOLUTE ERROR ' // 1083 * 'TOLERANCES MUST BE NON-NEGATIVE.$$YOU HAVE CALLED ' // 1084 * 'THE CODE WITH ATOL(' // XERN1 // ') = ' // XERN3 // 1085 * '$$IN THE CASE OF VECTOR ERROR TOLERANCES, NO FURTHER ' 1086 * // 'CHECKING OF ATOL COMPONENTS IS DONE.', 8, 1) 1087 IDID=-33 1088 IF (NRTOLP .GT. 0) GO TO 70 1089 NATOLP=1 1090 ENDIF 1091 50 IF (ITOL .EQ. 0) GO TO 70 1092 60 CONTINUE 1093C 1094 70 IF (ITSTOP .EQ. 1) THEN 1095 IF (SIGN(1.0D0,TOUT-T) .NE. SIGN(1.0D0,TSTOP-T) .OR. 1096 1 ABS(TOUT-T) .GT. ABS(TSTOP-T)) THEN 1097 WRITE (XERN3, '(1PE15.6)') TOUT 1098 WRITE (XERN4, '(1PE15.6)') TSTOP 1099 CALL XERMSG ('SLATEC', 'DLSOD', 1100 * 'IN DDEBDF, YOU HAVE CALLED THE ' // 1101 * 'CODE WITH TOUT = ' // XERN3 // '$$BUT YOU HAVE ' // 1102 * 'ALSO TOLD THE CODE NOT TO INTEGRATE PAST THE POINT ' // 1103 * 'TSTOP = ' // XERN4 // ' BY SETTING INFO(4) = 1.$$' // 1104 * 'THESE INSTRUCTIONS CONFLICT.', 14, 1) 1105 IDID=-33 1106 ENDIF 1107 ENDIF 1108C 1109C CHECK SOME CONTINUATION POSSIBILITIES 1110C 1111 IF (INIT .NE. 0) THEN 1112 IF (T .EQ. TOUT) THEN 1113 WRITE (XERN3, '(1PE15.6)') T 1114 CALL XERMSG ('SLATEC', 'DLSOD', 1115 * 'IN DDEBDF, YOU HAVE CALLED THE CODE WITH T = TOUT = ' // 1116 * XERN3 // '$$THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 1117 * 9, 1) 1118 IDID=-33 1119 ENDIF 1120C 1121 IF (T .NE. TOLD) THEN 1122 WRITE (XERN3, '(1PE15.6)') TOLD 1123 WRITE (XERN4, '(1PE15.6)') T 1124 CALL XERMSG ('SLATEC', 'DLSOD', 1125 * 'IN DDEBDF, YOU HAVE CHANGED THE VALUE OF T FROM ' // 1126 * XERN3 // ' TO ' // XERN4 // 1127 * ' THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 10, 1) 1128 IDID=-33 1129 ENDIF 1130C 1131 IF (INIT .NE. 1) THEN 1132 IF (DELSGN*(TOUT-T) .LT. 0.0D0) THEN 1133 WRITE (XERN3, '(1PE15.6)') TOUT 1134 CALL XERMSG ('SLATEC', 'DLSOD', 1135 * 'IN DDEBDF, BY CALLING THE CODE WITH TOUT = ' // 1136 * XERN3 // ' YOU ARE ATTEMPTING TO CHANGE THE ' // 1137 * 'DIRECTION OF INTEGRATION.$$THIS IS NOT ALLOWED ' // 1138 * 'WITHOUT RESTARTING.', 11, 1) 1139 IDID=-33 1140 ENDIF 1141 ENDIF 1142 ENDIF 1143C 1144 IF (IDID .EQ. (-33)) THEN 1145 IF (IQUIT .NE. (-33)) THEN 1146C INVALID INPUT DETECTED 1147 IQUIT=-33 1148 IBEGIN=-1 1149 ELSE 1150 CALL XERMSG ('SLATEC', 'DLSOD', 1151 * 'IN DDEBDF, INVALID INPUT WAS DETECTED ON ' // 1152 * 'SUCCESSIVE ENTRIES. IT IS IMPOSSIBLE TO PROCEED ' // 1153 * 'BECAUSE YOU HAVE NOT CORRECTED THE PROBLEM, ' // 1154 * 'SO EXECUTION IS BEING TERMINATED.', 12, 2) 1155 ENDIF 1156 RETURN 1157 ENDIF 1158C 1159C ............................................................... 1160C 1161C RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED 1162C AS ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS 1163C CASE, THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE 1164C SMALLEST VALUE 100*U WHICH IS LIKELY TO BE REASONABLE FOR 1165C THIS METHOD AND MACHINE 1166C 1167 DO 180 K = 1, NEQ 1168 IF (RTOL(K) + ATOL(K) .GT. 0.0D0) GO TO 170 1169 RTOL(K) = 100.0D0*U 1170 IDID = -2 1171 170 CONTINUE 1172C ...EXIT 1173 IF (ITOL .EQ. 0) GO TO 190 1174 180 CONTINUE 1175 190 CONTINUE 1176C 1177 IF (IDID .NE. (-2)) GO TO 200 1178C RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A 1179C SMALL POSITIVE VALUE 1180 IBEGIN = -1 1181 GO TO 460 1182 200 CONTINUE 1183C BEGIN BLOCK PERMITTING ...EXITS TO 450 1184C BEGIN BLOCK PERMITTING ...EXITS TO 430 1185C BEGIN BLOCK PERMITTING ...EXITS TO 260 1186C BEGIN BLOCK PERMITTING ...EXITS TO 230 1187C 1188C BRANCH ON STATUS OF INITIALIZATION INDICATOR 1189C INIT=0 MEANS INITIAL DERIVATIVES AND 1190C NOMINAL STEP SIZE 1191C AND DIRECTION NOT YET SET 1192C INIT=1 MEANS NOMINAL STEP SIZE AND 1193C DIRECTION NOT YET SET INIT=2 MEANS NO 1194C FURTHER INITIALIZATION REQUIRED 1195C 1196 IF (INIT .EQ. 0) GO TO 210 1197C ......EXIT 1198 IF (INIT .EQ. 1) GO TO 230 1199C .........EXIT 1200 GO TO 260 1201 210 CONTINUE 1202C 1203C ................................................ 1204C 1205C MORE INITIALIZATION -- 1206C -- EVALUATE INITIAL 1207C DERIVATIVES 1208C 1209 INIT = 1 1210 CALL DF(T,Y,YH(1,2),RPAR,IPAR) 1211 NFE = 1 1212C ...EXIT 1213 IF (T .NE. TOUT) GO TO 230 1214 IDID = 2 1215 DO 220 L = 1, NEQ 1216 YPOUT(L) = YH(L,2) 1217 220 CONTINUE 1218 TOLD = T 1219C ............EXIT 1220 GO TO 450 1221 230 CONTINUE 1222C 1223C -- COMPUTE INITIAL STEP SIZE 1224C -- SAVE SIGN OF INTEGRATION DIRECTION 1225C -- SET INDEPENDENT AND DEPENDENT VARIABLES 1226C X AND YH(*) FOR DSTOD 1227C 1228 LTOL = 1 1229 DO 240 L = 1, NEQ 1230 IF (ITOL .EQ. 1) LTOL = L 1231 TOL = RTOL(LTOL)*ABS(Y(L)) + ATOL(LTOL) 1232 IF (TOL .EQ. 0.0D0) GO TO 390 1233 EWT(L) = TOL 1234 240 CONTINUE 1235C 1236 BIG = SQRT(D1MACH(2)) 1237 CALL DHSTRT(DF,NEQ,T,TOUT,Y,YH(1,2),EWT,1,U,BIG, 1238 1 YH(1,3),YH(1,4),YH(1,5),YH(1,6),RPAR, 1239 2 IPAR,H) 1240C 1241 DELSGN = SIGN(1.0D0,TOUT-T) 1242 X = T 1243 DO 250 L = 1, NEQ 1244 YH(L,1) = Y(L) 1245 YH(L,2) = H*YH(L,2) 1246 250 CONTINUE 1247 INIT = 2 1248 260 CONTINUE 1249C 1250C ...................................................... 1251C 1252C ON EACH CALL SET INFORMATION WHICH DETERMINES THE 1253C ALLOWED INTERVAL OF INTEGRATION BEFORE RETURNING 1254C WITH AN ANSWER AT TOUT 1255C 1256 DEL = TOUT - T 1257 ABSDEL = ABS(DEL) 1258C 1259C ...................................................... 1260C 1261C IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND 1262C RETURN 1263C 1264 270 CONTINUE 1265C BEGIN BLOCK PERMITTING ...EXITS TO 400 1266C BEGIN BLOCK PERMITTING ...EXITS TO 380 1267 IF (ABS(X-T) .LT. ABSDEL) GO TO 290 1268 CALL DINTYD(TOUT,0,YH,NEQ,Y,INTFLG) 1269 CALL DINTYD(TOUT,1,YH,NEQ,YPOUT,INTFLG) 1270 IDID = 3 1271 IF (X .NE. TOUT) GO TO 280 1272 IDID = 2 1273 INTOUT = .FALSE. 1274 280 CONTINUE 1275 T = TOUT 1276 TOLD = T 1277C ..................EXIT 1278 GO TO 450 1279 290 CONTINUE 1280C 1281C IF CANNOT GO PAST TSTOP AND SUFFICIENTLY 1282C CLOSE, EXTRAPOLATE AND RETURN 1283C 1284 IF (ITSTOP .NE. 1) GO TO 310 1285 IF (ABS(TSTOP-X) .GE. 100.0D0*U*ABS(X)) 1286 1 GO TO 310 1287 DT = TOUT - X 1288 DO 300 L = 1, NEQ 1289 Y(L) = YH(L,1) + (DT/H)*YH(L,2) 1290 300 CONTINUE 1291 CALL DF(TOUT,Y,YPOUT,RPAR,IPAR) 1292 NFE = NFE + 1 1293 IDID = 3 1294 T = TOUT 1295 TOLD = T 1296C ..................EXIT 1297 GO TO 450 1298 310 CONTINUE 1299C 1300 IF (IINTEG .EQ. 0 .OR. .NOT.INTOUT) GO TO 320 1301C 1302C INTERMEDIATE-OUTPUT MODE 1303C 1304 IDID = 1 1305 GO TO 370 1306 320 CONTINUE 1307C 1308C ............................................. 1309C 1310C MONITOR NUMBER OF STEPS ATTEMPTED 1311C 1312 IF (KSTEPS .LE. MAXNUM) GO TO 330 1313C 1314C A SIGNIFICANT AMOUNT OF WORK HAS BEEN 1315C EXPENDED 1316 IDID = -1 1317 KSTEPS = 0 1318 IBEGIN = -1 1319 GO TO 370 1320 330 CONTINUE 1321C 1322C .......................................... 1323C 1324C LIMIT STEP SIZE AND SET WEIGHT VECTOR 1325C 1326 HMIN = 100.0D0*U*ABS(X) 1327 HA = MAX(ABS(H),HMIN) 1328 IF (ITSTOP .EQ. 1) 1329 1 HA = MIN(HA,ABS(TSTOP-X)) 1330 H = SIGN(HA,H) 1331 LTOL = 1 1332 DO 340 L = 1, NEQ 1333 IF (ITOL .EQ. 1) LTOL = L 1334 EWT(L) = RTOL(LTOL)*ABS(YH(L,1)) 1335 1 + ATOL(LTOL) 1336C .........EXIT 1337 IF (EWT(L) .LE. 0.0D0) GO TO 380 1338 340 CONTINUE 1339 TOLFAC = U*DVNRMS(NEQ,YH,EWT) 1340C .........EXIT 1341 IF (TOLFAC .LE. 1.0D0) GO TO 400 1342C 1343C TOLERANCES TOO SMALL 1344 IDID = -2 1345 TOLFAC = 2.0D0*TOLFAC 1346 RTOL(1) = TOLFAC*RTOL(1) 1347 ATOL(1) = TOLFAC*ATOL(1) 1348 IF (ITOL .EQ. 0) GO TO 360 1349 DO 350 L = 2, NEQ 1350 RTOL(L) = TOLFAC*RTOL(L) 1351 ATOL(L) = TOLFAC*ATOL(L) 1352 350 CONTINUE 1353 360 CONTINUE 1354 IBEGIN = -1 1355 370 CONTINUE 1356C ............EXIT 1357 GO TO 430 1358 380 CONTINUE 1359C 1360C RELATIVE ERROR CRITERION INAPPROPRIATE 1361 390 CONTINUE 1362 IDID = -3 1363 IBEGIN = -1 1364C .........EXIT 1365 GO TO 430 1366 400 CONTINUE 1367C 1368C ................................................... 1369C 1370C TAKE A STEP 1371C 1372 CALL DSTOD(NEQ,Y,YH,NEQ,YH1,EWT,SAVF,ACOR,WM,IWM, 1373 1 DF,DJAC,RPAR,IPAR) 1374C 1375 JSTART = -2 1376 INTOUT = .TRUE. 1377 IF (KFLAG .EQ. 0) GO TO 270 1378C 1379C ...................................................... 1380C 1381 IF (KFLAG .EQ. -1) GO TO 410 1382C 1383C REPEATED CORRECTOR CONVERGENCE FAILURES 1384 IDID = -6 1385 IBEGIN = -1 1386 GO TO 420 1387 410 CONTINUE 1388C 1389C REPEATED ERROR TEST FAILURES 1390 IDID = -7 1391 IBEGIN = -1 1392 420 CONTINUE 1393 430 CONTINUE 1394C 1395C ......................................................... 1396C 1397C STORE VALUES BEFORE RETURNING TO 1398C DDEBDF 1399 DO 440 L = 1, NEQ 1400 Y(L) = YH(L,1) 1401 YPOUT(L) = YH(L,2)/H 1402 440 CONTINUE 1403 T = X 1404 TOLD = T 1405 INTOUT = .FALSE. 1406 450 CONTINUE 1407 460 CONTINUE 1408 RETURN 1409 END 1410*DECK DSTOD 1411 SUBROUTINE DSTOD (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, 1412 + DF, DJAC, RPAR, IPAR) 1413C***BEGIN PROLOGUE DSTOD 1414C***SUBSIDIARY 1415C***PURPOSE Subsidiary to DDEBDF 1416C***LIBRARY SLATEC 1417C***TYPE DOUBLE PRECISION (STOD-S, DSTOD-D) 1418C***AUTHOR Watts, H. A., (SNLA) 1419C***DESCRIPTION 1420C 1421C DSTOD integrates a system of first order odes over one step in the 1422C integrator package DDEBDF. 1423C ---------------------------------------------------------------------- 1424C DSTOD performs one step of the integration of an initial value 1425C problem for a system of ordinary differential equations. 1426C Note.. DSTOD is independent of the value of the iteration method 1427C indicator MITER, when this is .NE. 0, and hence is independent 1428C of the type of chord method used, or the Jacobian structure. 1429C Communication with DSTOD is done with the following variables.. 1430C 1431C Y = An array of length .GE. N used as the Y argument in 1432C all calls to DF and DJAC. 1433C NEQ = Integer array containing problem size in NEQ(1), and 1434C passed as the NEQ argument in all calls to DF and DJAC. 1435C YH = An NYH by LMAX array containing the dependent variables 1436C and their approximate scaled derivatives, where 1437C LMAX = MAXORD + 1. YH(I,J+1) contains the approximate 1438C J-th derivative of Y(I), scaled by H**J/FACTORIAL(J) 1439C (J = 0,1,...,NQ). On entry for the first step, the first 1440C two columns of YH must be set from the initial values. 1441C NYH = A constant integer .GE. N, the first dimension of YH. 1442C YH1 = A one-dimensional array occupying the same space as YH. 1443C EWT = An array of N elements with which the estimated local 1444C errors in YH are compared. 1445C SAVF = An array of working storage, of length N. 1446C ACOR = A work array of length N, used for the accumulated 1447C corrections. On a successful return, ACOR(I) contains 1448C the estimated one-step local error in Y(I). 1449C WM,IWM = DOUBLE PRECISION and INTEGER work arrays associated with 1450C matrix operations in chord iteration (MITER .NE. 0). 1451C DPJAC = Name of routine to evaluate and preprocess Jacobian matrix 1452C if a chord method is being used. 1453C DSLVS = Name of routine to solve linear system in chord iteration. 1454C H = The step size to be attempted on the next step. 1455C H is altered by the error control algorithm during the 1456C problem. H can be either positive or negative, but its 1457C sign must remain constant throughout the problem. 1458C HMIN = The minimum absolute value of the step size H to be used. 1459C HMXI = Inverse of the maximum absolute value of H to be used. 1460C HMXI = 0.0 is allowed and corresponds to an infinite HMAX. 1461C HMIN and HMXI may be changed at any time, but will not 1462C take effect until the next change of H is considered. 1463C TN = The independent variable. TN is updated on each step taken. 1464C JSTART = An integer used for input only, with the following 1465C values and meanings.. 1466C 0 Perform the first step. 1467C .GT.0 Take a new step continuing from the last. 1468C -1 Take the next step with a new value of H, MAXORD, 1469C N, METH, MITER, and/or matrix parameters. 1470C -2 Take the next step with a new value of H, 1471C but with other inputs unchanged. 1472C On return, JSTART is set to 1 to facilitate continuation. 1473C KFLAG = a completion code with the following meanings.. 1474C 0 The step was successful. 1475C -1 The requested error could not be achieved. 1476C -2 Corrector convergence could not be achieved. 1477C A return with KFLAG = -1 or -2 means either 1478C ABS(H) = HMIN or 10 consecutive failures occurred. 1479C On a return with KFLAG negative, the values of TN and 1480C the YH array are as of the beginning of the last 1481C step, and H is the last step size attempted. 1482C MAXORD = The maximum order of integration method to be allowed. 1483C METH/MITER = The method flags. See description in driver. 1484C N = The number of first-order differential equations. 1485C ---------------------------------------------------------------------- 1486C 1487C***SEE ALSO DDEBDF 1488C***ROUTINES CALLED DCFOD, DPJAC, DSLVS, DVNRMS 1489C***COMMON BLOCKS DDEBD1 1490C***REVISION HISTORY (YYMMDD) 1491C 820301 DATE WRITTEN 1492C 890531 Changed all specific intrinsics to generic. (WRB) 1493C 890911 Removed unnecessary intrinsics. (WRB) 1494C 891214 Prologue converted to Version 4.0 format. (BAB) 1495C 900328 Added TYPE section. (WRB) 1496C 910722 Updated AUTHOR section. (ALS) 1497C 920422 Changed DIMENSION statement. (WRB) 1498C***END PROLOGUE DSTOD 1499C 1500 INTEGER I, I1, IALTH, IER, IOD, IOWND, IPAR, IPUP, IREDO, IRET, 1501 1 IWM, J, JB, JSTART, KFLAG, KSTEPS, L, LMAX, M, MAXORD, 1502 2 MEO, METH, MITER, N, NCF, NEQ, NEWQ, NFE, NJE, NQ, NQNYH, 1503 3 NQU, NST, NSTEPJ, NYH 1504 DOUBLE PRECISION ACOR, CONIT, CRATE, DCON, DDN, 1505 1 DEL, DELP, DSM, DUP, DVNRMS, EL, EL0, ELCO, 1506 2 EWT, EXDN, EXSM, EXUP, H, HMIN, HMXI, HOLD, HU, R, RC, 1507 3 RH, RHDN, RHSM, RHUP, RMAX, ROWND, RPAR, SAVF, TESCO, 1508 4 TN, TOLD, UROUND, WM, Y, YH, YH1 1509 EXTERNAL DF, DJAC 1510C 1511 DIMENSION Y(*),YH(NYH,*),YH1(*),EWT(*),SAVF(*),ACOR(*),WM(*), 1512 1 IWM(*),RPAR(*),IPAR(*) 1513 COMMON /DDEBD1/ ROWND,CONIT,CRATE,EL(13),ELCO(13,12),HOLD,RC,RMAX, 1514 1 TESCO(3,12),EL0,H,HMIN,HMXI,HU,TN,UROUND,IOWND(7), 1515 2 KSTEPS,IOD(6),IALTH,IPUP,LMAX,MEO,NQNYH,NSTEPJ, 1516 3 IER,JSTART,KFLAG,L,METH,MITER,MAXORD,N,NQ,NST,NFE, 1517 4 NJE,NQU 1518C 1519C 1520C BEGIN BLOCK PERMITTING ...EXITS TO 690 1521C BEGIN BLOCK PERMITTING ...EXITS TO 60 1522C***FIRST EXECUTABLE STATEMENT DSTOD 1523 KFLAG = 0 1524 TOLD = TN 1525 NCF = 0 1526 IF (JSTART .GT. 0) GO TO 160 1527 IF (JSTART .EQ. -1) GO TO 10 1528 IF (JSTART .EQ. -2) GO TO 90 1529C --------------------------------------------------------- 1530C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER 1531C VARIABLES ARE INITIALIZED. RMAX IS THE MAXIMUM RATIO BY 1532C WHICH H CAN BE INCREASED IN A SINGLE STEP. IT IS 1533C INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL INITIAL H, 1534C BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE OCCURS 1535C (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 1536C 2 FOR THE NEXT INCREASE. 1537C --------------------------------------------------------- 1538 LMAX = MAXORD + 1 1539 NQ = 1 1540 L = 2 1541 IALTH = 2 1542 RMAX = 10000.0D0 1543 RC = 0.0D0 1544 EL0 = 1.0D0 1545 CRATE = 0.7D0 1546 DELP = 0.0D0 1547 HOLD = H 1548 MEO = METH 1549 NSTEPJ = 0 1550 IRET = 3 1551 GO TO 50 1552 10 CONTINUE 1553C BEGIN BLOCK PERMITTING ...EXITS TO 30 1554C ------------------------------------------------------ 1555C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN 1556C JSTART = -1. IPUP IS SET TO MITER TO FORCE A MATRIX 1557C UPDATE. IF AN ORDER INCREASE IS ABOUT TO BE 1558C CONSIDERED (IALTH = 1), IALTH IS RESET TO 2 TO 1559C POSTPONE CONSIDERATION ONE MORE STEP. IF THE CALLER 1560C HAS CHANGED METH, DCFOD IS CALLED TO RESET THE 1561C COEFFICIENTS OF THE METHOD. IF THE CALLER HAS 1562C CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT 1563C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN 1564C ACCORDINGLY. IF H IS TO BE CHANGED, YH MUST BE 1565C RESCALED. IF H OR METH IS BEING CHANGED, IALTH IS 1566C RESET TO L = NQ + 1 TO PREVENT FURTHER CHANGES IN H 1567C FOR THAT MANY STEPS. 1568C ------------------------------------------------------ 1569 IPUP = MITER 1570 LMAX = MAXORD + 1 1571 IF (IALTH .EQ. 1) IALTH = 2 1572 IF (METH .EQ. MEO) GO TO 20 1573 CALL DCFOD(METH,ELCO,TESCO) 1574 MEO = METH 1575C ......EXIT 1576 IF (NQ .GT. MAXORD) GO TO 30 1577 IALTH = L 1578 IRET = 1 1579C ............EXIT 1580 GO TO 60 1581 20 CONTINUE 1582 IF (NQ .LE. MAXORD) GO TO 90 1583 30 CONTINUE 1584 NQ = MAXORD 1585 L = LMAX 1586 DO 40 I = 1, L 1587 EL(I) = ELCO(I,NQ) 1588 40 CONTINUE 1589 NQNYH = NQ*NYH 1590 RC = RC*EL(1)/EL0 1591 EL0 = EL(1) 1592 CONIT = 0.5D0/(NQ+2) 1593 DDN = DVNRMS(N,SAVF,EWT)/TESCO(1,L) 1594 EXDN = 1.0D0/L 1595 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 1596 RH = MIN(RHDN,1.0D0) 1597 IREDO = 3 1598 IF (H .EQ. HOLD) GO TO 660 1599 RH = MIN(RH,ABS(H/HOLD)) 1600 H = HOLD 1601 GO TO 100 1602 50 CONTINUE 1603C ------------------------------------------------------------ 1604C DCFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS 1605C FOR THE CURRENT METH. THEN THE EL VECTOR AND RELATED 1606C CONSTANTS ARE RESET WHENEVER THE ORDER NQ IS CHANGED, OR AT 1607C THE START OF THE PROBLEM. 1608C ------------------------------------------------------------ 1609 CALL DCFOD(METH,ELCO,TESCO) 1610 60 CONTINUE 1611 70 CONTINUE 1612C BEGIN BLOCK PERMITTING ...EXITS TO 680 1613 DO 80 I = 1, L 1614 EL(I) = ELCO(I,NQ) 1615 80 CONTINUE 1616 NQNYH = NQ*NYH 1617 RC = RC*EL(1)/EL0 1618 EL0 = EL(1) 1619 CONIT = 0.5D0/(NQ+2) 1620 GO TO (90,660,160), IRET 1621C --------------------------------------------------------- 1622C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST 1623C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH 1624C IS SET TO L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT 1625C MANY STEPS, UNLESS FORCED BY A CONVERGENCE OR ERROR TEST 1626C FAILURE. 1627C --------------------------------------------------------- 1628 90 CONTINUE 1629 IF (H .EQ. HOLD) GO TO 160 1630 RH = H/HOLD 1631 H = HOLD 1632 IREDO = 3 1633 100 CONTINUE 1634 110 CONTINUE 1635 RH = MIN(RH,RMAX) 1636 RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH) 1637 R = 1.0D0 1638 DO 130 J = 2, L 1639 R = R*RH 1640 DO 120 I = 1, N 1641 YH(I,J) = YH(I,J)*R 1642 120 CONTINUE 1643 130 CONTINUE 1644 H = H*RH 1645 RC = RC*RH 1646 IALTH = L 1647 IF (IREDO .NE. 0) GO TO 150 1648 RMAX = 10.0D0 1649 R = 1.0D0/TESCO(2,NQU) 1650 DO 140 I = 1, N 1651 ACOR(I) = ACOR(I)*R 1652 140 CONTINUE 1653C ...............EXIT 1654 GO TO 690 1655 150 CONTINUE 1656C ------------------------------------------------------ 1657C THIS SECTION COMPUTES THE PREDICTED VALUES BY 1658C EFFECTIVELY MULTIPLYING THE YH ARRAY BY THE PASCAL 1659C TRIANGLE MATRIX. RC IS THE RATIO OF NEW TO OLD 1660C VALUES OF THE COEFFICIENT H*EL(1). WHEN RC DIFFERS 1661C FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER 1662C TO FORCE DPJAC TO BE CALLED, IF A JACOBIAN IS 1663C INVOLVED. IN ANY CASE, DPJAC IS CALLED AT LEAST 1664C EVERY 20-TH STEP. 1665C ------------------------------------------------------ 1666 160 CONTINUE 1667 170 CONTINUE 1668C BEGIN BLOCK PERMITTING ...EXITS TO 610 1669C BEGIN BLOCK PERMITTING ...EXITS TO 490 1670 IF (ABS(RC-1.0D0) .GT. 0.3D0) IPUP = MITER 1671 IF (NST .GE. NSTEPJ + 20) IPUP = MITER 1672 TN = TN + H 1673 I1 = NQNYH + 1 1674 DO 190 JB = 1, NQ 1675 I1 = I1 - NYH 1676 DO 180 I = I1, NQNYH 1677 YH1(I) = YH1(I) + YH1(I+NYH) 1678 180 CONTINUE 1679 190 CONTINUE 1680 KSTEPS = KSTEPS + 1 1681C --------------------------------------------- 1682C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A 1683C CONVERGENCE TEST IS MADE ON THE R.M.S. NORM 1684C OF EACH CORRECTION, WEIGHTED BY THE ERROR 1685C WEIGHT VECTOR EWT. THE SUM OF THE 1686C CORRECTIONS IS ACCUMULATED IN THE VECTOR 1687C ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE 1688C CORRECTOR LOOP. 1689C --------------------------------------------- 1690 200 CONTINUE 1691 M = 0 1692 DO 210 I = 1, N 1693 Y(I) = YH(I,1) 1694 210 CONTINUE 1695 CALL DF(TN,Y,SAVF,RPAR,IPAR) 1696 NFE = NFE + 1 1697 IF (IPUP .LE. 0) GO TO 220 1698C --------------------------------------- 1699C IF INDICATED, THE MATRIX P = I - 1700C H*EL(1)*J IS REEVALUATED AND 1701C PREPROCESSED BEFORE STARTING THE 1702C CORRECTOR ITERATION. IPUP IS SET TO 0 1703C AS AN INDICATOR THAT THIS HAS BEEN 1704C DONE. 1705C --------------------------------------- 1706 IPUP = 0 1707 RC = 1.0D0 1708 NSTEPJ = NST 1709 CRATE = 0.7D0 1710 CALL DPJAC(NEQ,Y,YH,NYH,EWT,ACOR,SAVF, 1711 1 WM,IWM,DF,DJAC,RPAR,IPAR) 1712C ......EXIT 1713 IF (IER .NE. 0) GO TO 440 1714 220 CONTINUE 1715 DO 230 I = 1, N 1716 ACOR(I) = 0.0D0 1717 230 CONTINUE 1718 240 CONTINUE 1719 IF (MITER .NE. 0) GO TO 270 1720C ------------------------------------ 1721C IN THE CASE OF FUNCTIONAL 1722C ITERATION, UPDATE Y DIRECTLY FROM 1723C THE RESULT OF THE LAST FUNCTION 1724C EVALUATION. 1725C ------------------------------------ 1726 DO 250 I = 1, N 1727 SAVF(I) = H*SAVF(I) - YH(I,2) 1728 Y(I) = SAVF(I) - ACOR(I) 1729 250 CONTINUE 1730 DEL = DVNRMS(N,Y,EWT) 1731 DO 260 I = 1, N 1732 Y(I) = YH(I,1) + EL(1)*SAVF(I) 1733 ACOR(I) = SAVF(I) 1734 260 CONTINUE 1735 GO TO 300 1736 270 CONTINUE 1737C ------------------------------------ 1738C IN THE CASE OF THE CHORD METHOD, 1739C COMPUTE THE CORRECTOR ERROR, AND 1740C SOLVE THE LINEAR SYSTEM WITH THAT 1741C AS RIGHT-HAND SIDE AND P AS 1742C COEFFICIENT MATRIX. 1743C ------------------------------------ 1744 DO 280 I = 1, N 1745 Y(I) = H*SAVF(I) 1746 1 - (YH(I,2) + ACOR(I)) 1747 280 CONTINUE 1748 CALL DSLVS(WM,IWM,Y,SAVF) 1749C ......EXIT 1750 IF (IER .NE. 0) GO TO 430 1751 DEL = DVNRMS(N,Y,EWT) 1752 DO 290 I = 1, N 1753 ACOR(I) = ACOR(I) + Y(I) 1754 Y(I) = YH(I,1) + EL(1)*ACOR(I) 1755 290 CONTINUE 1756 300 CONTINUE 1757C --------------------------------------- 1758C TEST FOR CONVERGENCE. IF M.GT.0, AN 1759C ESTIMATE OF THE CONVERGENCE RATE 1760C CONSTANT IS STORED IN CRATE, AND THIS 1761C IS USED IN THE TEST. 1762C --------------------------------------- 1763 IF (M .NE. 0) 1764 1 CRATE = MAX(0.2D0*CRATE,DEL/DELP) 1765 DCON = DEL*MIN(1.0D0,1.5D0*CRATE) 1766 1 /(TESCO(2,NQ)*CONIT) 1767 IF (DCON .GT. 1.0D0) GO TO 420 1768C ------------------------------------ 1769C THE CORRECTOR HAS CONVERGED. IPUP 1770C IS SET TO -1 IF MITER .NE. 0, TO 1771C SIGNAL THAT THE JACOBIAN INVOLVED 1772C MAY NEED UPDATING LATER. THE LOCAL 1773C ERROR TEST IS MADE AND CONTROL 1774C PASSES TO STATEMENT 500 IF IT 1775C FAILS. 1776C ------------------------------------ 1777 IF (MITER .NE. 0) IPUP = -1 1778 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 1779 IF (M .GT. 0) 1780 1 DSM = DVNRMS(N,ACOR,EWT) 1781 2 /TESCO(2,NQ) 1782 IF (DSM .GT. 1.0D0) GO TO 380 1783C BEGIN BLOCK 1784C PERMITTING ...EXITS TO 360 1785C ------------------------------ 1786C AFTER A SUCCESSFUL STEP, 1787C UPDATE THE YH ARRAY. 1788C CONSIDER CHANGING H IF IALTH 1789C = 1. OTHERWISE DECREASE 1790C IALTH BY 1. IF IALTH IS THEN 1791C 1 AND NQ .LT. MAXORD, THEN 1792C ACOR IS SAVED FOR USE IN A 1793C POSSIBLE ORDER INCREASE ON 1794C THE NEXT STEP. IF A CHANGE 1795C IN H IS CONSIDERED, AN 1796C INCREASE OR DECREASE IN ORDER 1797C BY ONE IS CONSIDERED ALSO. A 1798C CHANGE IN H IS MADE ONLY IF 1799C IT IS BY A FACTOR OF AT LEAST 1800C 1.1. IF NOT, IALTH IS SET TO 1801C 3 TO PREVENT TESTING FOR THAT 1802C MANY STEPS. 1803C ------------------------------ 1804 KFLAG = 0 1805 IREDO = 0 1806 NST = NST + 1 1807 HU = H 1808 NQU = NQ 1809 DO 320 J = 1, L 1810 DO 310 I = 1, N 1811 YH(I,J) = YH(I,J) 1812 1 + EL(J) 1813 2 *ACOR(I) 1814 310 CONTINUE 1815 320 CONTINUE 1816 IALTH = IALTH - 1 1817 IF (IALTH .NE. 0) GO TO 340 1818C --------------------------- 1819C REGARDLESS OF THE SUCCESS 1820C OR FAILURE OF THE STEP, 1821C FACTORS RHDN, RHSM, AND 1822C RHUP ARE COMPUTED, BY 1823C WHICH H COULD BE 1824C MULTIPLIED AT ORDER NQ - 1825C 1, ORDER NQ, OR ORDER NQ + 1826C 1, RESPECTIVELY. IN THE 1827C CASE OF FAILURE, RHUP = 1828C 0.0 TO AVOID AN ORDER 1829C INCREASE. THE LARGEST OF 1830C THESE IS DETERMINED AND 1831C THE NEW ORDER CHOSEN 1832C ACCORDINGLY. IF THE ORDER 1833C IS TO BE INCREASED, WE 1834C COMPUTE ONE ADDITIONAL 1835C SCALED DERIVATIVE. 1836C --------------------------- 1837 RHUP = 0.0D0 1838C .....................EXIT 1839 IF (L .EQ. LMAX) GO TO 490 1840 DO 330 I = 1, N 1841 SAVF(I) = ACOR(I) 1842 1 - YH(I,LMAX) 1843 330 CONTINUE 1844 DUP = DVNRMS(N,SAVF,EWT) 1845 1 /TESCO(3,NQ) 1846 EXUP = 1.0D0/(L+1) 1847 RHUP = 1.0D0 1848 1 /(1.4D0*DUP**EXUP 1849 2 + 0.0000014D0) 1850C .....................EXIT 1851 GO TO 490 1852 340 CONTINUE 1853C ...EXIT 1854 IF (IALTH .GT. 1) GO TO 360 1855C ...EXIT 1856 IF (L .EQ. LMAX) GO TO 360 1857 DO 350 I = 1, N 1858 YH(I,LMAX) = ACOR(I) 1859 350 CONTINUE 1860 360 CONTINUE 1861 R = 1.0D0/TESCO(2,NQU) 1862 DO 370 I = 1, N 1863 ACOR(I) = ACOR(I)*R 1864 370 CONTINUE 1865C .................................EXIT 1866 GO TO 690 1867 380 CONTINUE 1868C ------------------------------------ 1869C THE ERROR TEST FAILED. KFLAG KEEPS 1870C TRACK OF MULTIPLE FAILURES. 1871C RESTORE TN AND THE YH ARRAY TO 1872C THEIR PREVIOUS VALUES, AND PREPARE 1873C TO TRY THE STEP AGAIN. COMPUTE THE 1874C OPTIMUM STEP SIZE FOR THIS OR ONE 1875C LOWER ORDER. AFTER 2 OR MORE 1876C FAILURES, H IS FORCED TO DECREASE 1877C BY A FACTOR OF 0.2 OR LESS. 1878C ------------------------------------ 1879 KFLAG = KFLAG - 1 1880 TN = TOLD 1881 I1 = NQNYH + 1 1882 DO 400 JB = 1, NQ 1883 I1 = I1 - NYH 1884 DO 390 I = I1, NQNYH 1885 YH1(I) = YH1(I) - YH1(I+NYH) 1886 390 CONTINUE 1887 400 CONTINUE 1888 RMAX = 2.0D0 1889 IF (ABS(H) .GT. HMIN*1.00001D0) 1890 1 GO TO 410 1891C --------------------------------- 1892C ALL RETURNS ARE MADE THROUGH 1893C THIS SECTION. H IS SAVED IN 1894C HOLD TO ALLOW THE CALLER TO 1895C CHANGE H ON THE NEXT STEP. 1896C --------------------------------- 1897 KFLAG = -1 1898C .................................EXIT 1899 GO TO 690 1900 410 CONTINUE 1901C ...............EXIT 1902 IF (KFLAG .LE. -3) GO TO 610 1903 IREDO = 2 1904 RHUP = 0.0D0 1905C ............EXIT 1906 GO TO 490 1907 420 CONTINUE 1908 M = M + 1 1909C ...EXIT 1910 IF (M .EQ. 3) GO TO 430 1911C ...EXIT 1912 IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) 1913 1 GO TO 430 1914 DELP = DEL 1915 CALL DF(TN,Y,SAVF,RPAR,IPAR) 1916 NFE = NFE + 1 1917 GO TO 240 1918 430 CONTINUE 1919C ------------------------------------------ 1920C THE CORRECTOR ITERATION FAILED TO 1921C CONVERGE IN 3 TRIES. IF MITER .NE. 0 AND 1922C THE JACOBIAN IS OUT OF DATE, DPJAC IS 1923C CALLED FOR THE NEXT TRY. OTHERWISE THE 1924C YH ARRAY IS RETRACTED TO ITS VALUES 1925C BEFORE PREDICTION, AND H IS REDUCED, IF 1926C POSSIBLE. IF H CANNOT BE REDUCED OR 10 1927C FAILURES HAVE OCCURRED, EXIT WITH KFLAG = 1928C -2. 1929C ------------------------------------------ 1930C ...EXIT 1931 IF (IPUP .EQ. 0) GO TO 440 1932 IPUP = MITER 1933 GO TO 200 1934 440 CONTINUE 1935 TN = TOLD 1936 NCF = NCF + 1 1937 RMAX = 2.0D0 1938 I1 = NQNYH + 1 1939 DO 460 JB = 1, NQ 1940 I1 = I1 - NYH 1941 DO 450 I = I1, NQNYH 1942 YH1(I) = YH1(I) - YH1(I+NYH) 1943 450 CONTINUE 1944 460 CONTINUE 1945 IF (ABS(H) .GT. HMIN*1.00001D0) GO TO 470 1946 KFLAG = -2 1947C ........................EXIT 1948 GO TO 690 1949 470 CONTINUE 1950 IF (NCF .NE. 10) GO TO 480 1951 KFLAG = -2 1952C ........................EXIT 1953 GO TO 690 1954 480 CONTINUE 1955 RH = 0.25D0 1956 IPUP = MITER 1957 IREDO = 1 1958C .........EXIT 1959 GO TO 650 1960 490 CONTINUE 1961 EXSM = 1.0D0/L 1962 RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0) 1963 RHDN = 0.0D0 1964 IF (NQ .EQ. 1) GO TO 500 1965 DDN = DVNRMS(N,YH(1,L),EWT)/TESCO(1,NQ) 1966 EXDN = 1.0D0/NQ 1967 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 1968 500 CONTINUE 1969 IF (RHSM .GE. RHUP) GO TO 550 1970 IF (RHUP .LE. RHDN) GO TO 540 1971 NEWQ = L 1972 RH = RHUP 1973 IF (RH .GE. 1.1D0) GO TO 520 1974 IALTH = 3 1975 R = 1.0D0/TESCO(2,NQU) 1976 DO 510 I = 1, N 1977 ACOR(I) = ACOR(I)*R 1978 510 CONTINUE 1979C ...........................EXIT 1980 GO TO 690 1981 520 CONTINUE 1982 R = EL(L)/L 1983 DO 530 I = 1, N 1984 YH(I,NEWQ+1) = ACOR(I)*R 1985 530 CONTINUE 1986 NQ = NEWQ 1987 L = NQ + 1 1988 IRET = 2 1989C ..................EXIT 1990 GO TO 680 1991 540 CONTINUE 1992 GO TO 580 1993 550 CONTINUE 1994 IF (RHSM .LT. RHDN) GO TO 580 1995 NEWQ = NQ 1996 RH = RHSM 1997 IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) 1998 1 GO TO 560 1999 IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) 2000C ------------------------------------------ 2001C IF THERE IS A CHANGE OF ORDER, RESET NQ, 2002C L, AND THE COEFFICIENTS. IN ANY CASE H 2003C IS RESET ACCORDING TO RH AND THE YH ARRAY 2004C IS RESCALED. THEN EXIT FROM 680 IF THE 2005C STEP WAS OK, OR REDO THE STEP OTHERWISE. 2006C ------------------------------------------ 2007C ............EXIT 2008 IF (NEWQ .EQ. NQ) GO TO 650 2009 NQ = NEWQ 2010 L = NQ + 1 2011 IRET = 2 2012C ..................EXIT 2013 GO TO 680 2014 560 CONTINUE 2015 IALTH = 3 2016 R = 1.0D0/TESCO(2,NQU) 2017 DO 570 I = 1, N 2018 ACOR(I) = ACOR(I)*R 2019 570 CONTINUE 2020C .....................EXIT 2021 GO TO 690 2022 580 CONTINUE 2023 NEWQ = NQ - 1 2024 RH = RHDN 2025 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0 2026 IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) GO TO 590 2027 IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) 2028C --------------------------------------------- 2029C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, 2030C AND THE COEFFICIENTS. IN ANY CASE H IS 2031C RESET ACCORDING TO RH AND THE YH ARRAY IS 2032C RESCALED. THEN EXIT FROM 680 IF THE STEP 2033C WAS OK, OR REDO THE STEP OTHERWISE. 2034C --------------------------------------------- 2035C .........EXIT 2036 IF (NEWQ .EQ. NQ) GO TO 650 2037 NQ = NEWQ 2038 L = NQ + 1 2039 IRET = 2 2040C ...............EXIT 2041 GO TO 680 2042 590 CONTINUE 2043 IALTH = 3 2044 R = 1.0D0/TESCO(2,NQU) 2045 DO 600 I = 1, N 2046 ACOR(I) = ACOR(I)*R 2047 600 CONTINUE 2048C ..................EXIT 2049 GO TO 690 2050 610 CONTINUE 2051C --------------------------------------------------- 2052C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES 2053C HAVE OCCURRED. IF 10 FAILURES HAVE OCCURRED, EXIT 2054C WITH KFLAG = -1. IT IS ASSUMED THAT THE 2055C DERIVATIVES THAT HAVE ACCUMULATED IN THE YH ARRAY 2056C HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST 2057C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 2058C 1. THEN H IS REDUCED BY A FACTOR OF 10, AND THE 2059C STEP IS RETRIED, UNTIL IT SUCCEEDS OR H REACHES 2060C HMIN. 2061C --------------------------------------------------- 2062 IF (KFLAG .NE. -10) GO TO 620 2063C ------------------------------------------------ 2064C ALL RETURNS ARE MADE THROUGH THIS SECTION. H 2065C IS SAVED IN HOLD TO ALLOW THE CALLER TO CHANGE 2066C H ON THE NEXT STEP. 2067C ------------------------------------------------ 2068 KFLAG = -1 2069C ..................EXIT 2070 GO TO 690 2071 620 CONTINUE 2072 RH = 0.1D0 2073 RH = MAX(HMIN/ABS(H),RH) 2074 H = H*RH 2075 DO 630 I = 1, N 2076 Y(I) = YH(I,1) 2077 630 CONTINUE 2078 CALL DF(TN,Y,SAVF,RPAR,IPAR) 2079 NFE = NFE + 1 2080 DO 640 I = 1, N 2081 YH(I,2) = H*SAVF(I) 2082 640 CONTINUE 2083 IPUP = MITER 2084 IALTH = 5 2085C ......EXIT 2086 IF (NQ .NE. 1) GO TO 670 2087 GO TO 170 2088 650 CONTINUE 2089 660 CONTINUE 2090 RH = MAX(RH,HMIN/ABS(H)) 2091 GO TO 110 2092 670 CONTINUE 2093 NQ = 1 2094 L = 2 2095 IRET = 3 2096 680 CONTINUE 2097 GO TO 70 2098 690 CONTINUE 2099 HOLD = H 2100 JSTART = 1 2101 RETURN 2102C ----------------------- END OF SUBROUTINE DSTOD 2103C ----------------------- 2104 END 2105*DECK DCFOD 2106 SUBROUTINE DCFOD (METH, ELCO, TESCO) 2107C***BEGIN PROLOGUE DCFOD 2108C***SUBSIDIARY 2109C***PURPOSE Subsidiary to DDEBDF 2110C***LIBRARY SLATEC 2111C***TYPE DOUBLE PRECISION (CFOD-S, DCFOD-D) 2112C***AUTHOR (UNKNOWN) 2113C***DESCRIPTION 2114C 2115C DCFOD defines coefficients needed in the integrator package DDEBDF 2116C 2117C***SEE ALSO DDEBDF 2118C***ROUTINES CALLED (NONE) 2119C***REVISION HISTORY (YYMMDD) 2120C 820301 DATE WRITTEN 2121C 890911 Removed unnecessary intrinsics. (WRB) 2122C 891214 Prologue converted to Version 4.0 format. (BAB) 2123C 900328 Added TYPE section. (WRB) 2124C***END PROLOGUE DCFOD 2125C 2126C 2127 INTEGER I, IB, METH, NQ, NQM1, NQP1 2128 DOUBLE PRECISION AGAMQ, ELCO, FNQ, FNQM1, PC, PINT, RAGQ, 2129 1 RQ1FAC, RQFAC, TESCO, TSIGN, XPIN 2130 DIMENSION ELCO(13,12),TESCO(3,12) 2131C ------------------------------------------------------------------ 2132C DCFOD IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS 2133C NEEDED THERE. THE COEFFICIENTS FOR THE CURRENT METHOD, AS 2134C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED. 2135C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2136C 2. (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.) 2137C DCFOD IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM, 2138C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED. 2139C 2140C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS. 2141C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF 2142C ORDER NQ ARE STORED IN ELCO(I,NQ). THEY ARE GIVEN BY A 2143C GENERATING POLYNOMIAL, I.E., 2144C L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ. 2145C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY 2146C DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1), L(-1) = 2147C 0. FOR THE BDF METHODS, L(X) IS GIVEN BY 2148C L(X) = (X+1)*(X+2)* ... *(X+NQ)/K, 2149C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ). 2150C 2151C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE 2152C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER. 2153C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP 2154C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER 2155C NQ + 1 IF K = 3. 2156C ------------------------------------------------------------------ 2157 DIMENSION PC(12) 2158C 2159C***FIRST EXECUTABLE STATEMENT DCFOD 2160 GO TO (10,60), METH 2161C 2162 10 CONTINUE 2163 ELCO(1,1) = 1.0D0 2164 ELCO(2,1) = 1.0D0 2165 TESCO(1,1) = 0.0D0 2166 TESCO(2,1) = 2.0D0 2167 TESCO(1,2) = 1.0D0 2168 TESCO(3,12) = 0.0D0 2169 PC(1) = 1.0D0 2170 RQFAC = 1.0D0 2171 DO 50 NQ = 2, 12 2172C ------------------------------------------------------------ 2173C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE 2174C POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ-1). 2175C INITIALLY, P(X) = 1. 2176C ------------------------------------------------------------ 2177 RQ1FAC = RQFAC 2178 RQFAC = RQFAC/NQ 2179 NQM1 = NQ - 1 2180 FNQM1 = NQM1 2181 NQP1 = NQ + 1 2182C FORM COEFFICIENTS OF P(X)*(X+NQ-1). 2183C ---------------------------------- 2184 PC(NQ) = 0.0D0 2185 DO 20 IB = 1, NQM1 2186 I = NQP1 - IB 2187 PC(I) = PC(I-1) + FNQM1*PC(I) 2188 20 CONTINUE 2189 PC(1) = FNQM1*PC(1) 2190C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). 2191C ----------------------- 2192 PINT = PC(1) 2193 XPIN = PC(1)/2.0D0 2194 TSIGN = 1.0D0 2195 DO 30 I = 2, NQ 2196 TSIGN = -TSIGN 2197 PINT = PINT + TSIGN*PC(I)/I 2198 XPIN = XPIN + TSIGN*PC(I)/(I+1) 2199 30 CONTINUE 2200C STORE COEFFICIENTS IN ELCO AND TESCO. 2201C -------------------------------- 2202 ELCO(1,NQ) = PINT*RQ1FAC 2203 ELCO(2,NQ) = 1.0D0 2204 DO 40 I = 2, NQ 2205 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I 2206 40 CONTINUE 2207 AGAMQ = RQFAC*XPIN 2208 RAGQ = 1.0D0/AGAMQ 2209 TESCO(2,NQ) = RAGQ 2210 IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1 2211 TESCO(3,NQM1) = RAGQ 2212 50 CONTINUE 2213 GO TO 100 2214C 2215 60 CONTINUE 2216 PC(1) = 1.0D0 2217 RQ1FAC = 1.0D0 2218 DO 90 NQ = 1, 5 2219C ------------------------------------------------------------ 2220C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE 2221C POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ). 2222C INITIALLY, P(X) = 1. 2223C ------------------------------------------------------------ 2224 FNQ = NQ 2225 NQP1 = NQ + 1 2226C FORM COEFFICIENTS OF P(X)*(X+NQ). 2227C ------------------------------------ 2228 PC(NQP1) = 0.0D0 2229 DO 70 IB = 1, NQ 2230 I = NQ + 2 - IB 2231 PC(I) = PC(I-1) + FNQ*PC(I) 2232 70 CONTINUE 2233 PC(1) = FNQ*PC(1) 2234C STORE COEFFICIENTS IN ELCO AND TESCO. 2235C -------------------------------- 2236 DO 80 I = 1, NQP1 2237 ELCO(I,NQ) = PC(I)/PC(2) 2238 80 CONTINUE 2239 ELCO(2,NQ) = 1.0D0 2240 TESCO(1,NQ) = RQ1FAC 2241 TESCO(2,NQ) = NQP1/ELCO(1,NQ) 2242 TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ) 2243 RQ1FAC = RQ1FAC/FNQ 2244 90 CONTINUE 2245 100 CONTINUE 2246 RETURN 2247C ----------------------- END OF SUBROUTINE DCFOD 2248C ----------------------- 2249 END 2250*DECK DVNRMS 2251 DOUBLE PRECISION FUNCTION DVNRMS (N, V, W) 2252C***BEGIN PROLOGUE DVNRMS 2253C***SUBSIDIARY 2254C***PURPOSE Subsidiary to DDEBDF 2255C***LIBRARY SLATEC 2256C***TYPE DOUBLE PRECISION (VNWRMS-S, DVNRMS-D) 2257C***AUTHOR (UNKNOWN) 2258C***DESCRIPTION 2259C 2260C DVNRMS computes a weighted root-mean-square vector norm for the 2261C integrator package DDEBDF. 2262C 2263C***SEE ALSO DDEBDF 2264C***ROUTINES CALLED (NONE) 2265C***REVISION HISTORY (YYMMDD) 2266C 820301 DATE WRITTEN 2267C 890531 Changed all specific intrinsics to generic. (WRB) 2268C 890831 Modified array declarations. (WRB) 2269C 890911 Removed unnecessary intrinsics. (WRB) 2270C 891214 Prologue converted to Version 4.0 format. (BAB) 2271C 900328 Added TYPE section. (WRB) 2272C***END PROLOGUE DVNRMS 2273 INTEGER I, N 2274 DOUBLE PRECISION SUM, V, W 2275 DIMENSION V(*),W(*) 2276C***FIRST EXECUTABLE STATEMENT DVNRMS 2277 SUM = 0.0D0 2278 DO 10 I = 1, N 2279 SUM = SUM + (V(I)/W(I))**2 2280 10 CONTINUE 2281 DVNRMS = SQRT(SUM/N) 2282 RETURN 2283C ----------------------- END OF FUNCTION DVNRMS 2284C ------------------------ 2285 END 2286*DECK DINTYD 2287 SUBROUTINE DINTYD (T, K, YH, NYH, DKY, IFLAG) 2288C***BEGIN PROLOGUE DINTYD 2289C***SUBSIDIARY 2290C***PURPOSE Subsidiary to DDEBDF 2291C***LIBRARY SLATEC 2292C***TYPE DOUBLE PRECISION (INTYD-S, DINTYD-D) 2293C***AUTHOR Watts, H. A., (SNLA) 2294C***DESCRIPTION 2295C 2296C DINTYD approximates the solution and derivatives at T by polynomial 2297C interpolation. Must be used in conjunction with the integrator 2298C package DDEBDF. 2299C ---------------------------------------------------------------------- 2300C DINTYD computes interpolated values of the K-th derivative of the 2301C dependent variable vector Y, and stores it in DKY. 2302C This routine is called by DDEBDF with K = 0,1 and T = TOUT, but may 2303C also be called by the user for any K up to the current order. 2304C (see detailed instructions in LSODE usage documentation.) 2305C ---------------------------------------------------------------------- 2306C The computed values in DKY are gotten by interpolation using the 2307C Nordsieck history array YH. This array corresponds uniquely to a 2308C vector-valued polynomial of degree NQCUR or less, and DKY is set 2309C to the K-th derivative of this polynomial at T. 2310C The formula for DKY is.. 2311C Q 2312C DKY(I) = Sum C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1) 2313C J=K 2314C where C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR. 2315C The quantities NQ = NQCUR, L = NQ+1, N = NEQ, TN, and H are 2316C communicated by common. The above sum is done in reverse order. 2317C IFLAG is returned negative if either K or T is out of bounds. 2318C ---------------------------------------------------------------------- 2319C 2320C***SEE ALSO DDEBDF 2321C***ROUTINES CALLED (NONE) 2322C***COMMON BLOCKS DDEBD1 2323C***REVISION HISTORY (YYMMDD) 2324C 820301 DATE WRITTEN 2325C 890911 Removed unnecessary intrinsics. (WRB) 2326C 891214 Prologue converted to Version 4.0 format. (BAB) 2327C 900328 Added TYPE section. (WRB) 2328C 910722 Updated AUTHOR section. (ALS) 2329C***END PROLOGUE DINTYD 2330C 2331 INTEGER I, IC, IER, IFLAG, IOWND, IOWNS, J, JB, JB2, JJ, JJ1, 2332 1 JP1, JSTART, K, KFLAG, L, MAXORD, METH, MITER, N, NFE, 2333 2 NJE, NQ, NQU, NST, NYH 2334 DOUBLE PRECISION C, DKY, EL0, H, HMIN, HMXI, HU, R, ROWND, 2335 1 ROWNS, S, T, TN, TP, UROUND, YH 2336 DIMENSION YH(NYH,*),DKY(*) 2337 COMMON /DDEBD1/ ROWND,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND, 2338 1 IOWND(14),IOWNS(6),IER,JSTART,KFLAG,L,METH,MITER, 2339 2 MAXORD,N,NQ,NST,NFE,NJE,NQU 2340C 2341C BEGIN BLOCK PERMITTING ...EXITS TO 130 2342C***FIRST EXECUTABLE STATEMENT DINTYD 2343 IFLAG = 0 2344 IF (K .LT. 0 .OR. K .GT. NQ) GO TO 110 2345 TP = TN - HU*(1.0D0 + 100.0D0*UROUND) 2346 IF ((T - TP)*(T - TN) .LE. 0.0D0) GO TO 10 2347 IFLAG = -2 2348C .........EXIT 2349 GO TO 130 2350 10 CONTINUE 2351C 2352 S = (T - TN)/H 2353 IC = 1 2354 IF (K .EQ. 0) GO TO 30 2355 JJ1 = L - K 2356 DO 20 JJ = JJ1, NQ 2357 IC = IC*JJ 2358 20 CONTINUE 2359 30 CONTINUE 2360 C = IC 2361 DO 40 I = 1, N 2362 DKY(I) = C*YH(I,L) 2363 40 CONTINUE 2364 IF (K .EQ. NQ) GO TO 90 2365 JB2 = NQ - K 2366 DO 80 JB = 1, JB2 2367 J = NQ - JB 2368 JP1 = J + 1 2369 IC = 1 2370 IF (K .EQ. 0) GO TO 60 2371 JJ1 = JP1 - K 2372 DO 50 JJ = JJ1, J 2373 IC = IC*JJ 2374 50 CONTINUE 2375 60 CONTINUE 2376 C = IC 2377 DO 70 I = 1, N 2378 DKY(I) = C*YH(I,JP1) + S*DKY(I) 2379 70 CONTINUE 2380 80 CONTINUE 2381C .........EXIT 2382 IF (K .EQ. 0) GO TO 130 2383 90 CONTINUE 2384 R = H**(-K) 2385 DO 100 I = 1, N 2386 DKY(I) = R*DKY(I) 2387 100 CONTINUE 2388 GO TO 120 2389 110 CONTINUE 2390C 2391 IFLAG = -1 2392 120 CONTINUE 2393 130 CONTINUE 2394 RETURN 2395C ----------------------- END OF SUBROUTINE DINTYD 2396C ----------------------- 2397 END 2398*DECK DPJAC 2399 SUBROUTINE DPJAC (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, DF, 2400 + DJAC, RPAR, IPAR) 2401C***BEGIN PROLOGUE DPJAC 2402C***SUBSIDIARY 2403C***PURPOSE Subsidiary to DDEBDF 2404C***LIBRARY SLATEC 2405C***TYPE DOUBLE PRECISION (PJAC-S, DPJAC-D) 2406C***AUTHOR Watts, H. A., (SNLA) 2407C***DESCRIPTION 2408C 2409C DPJAC sets up the iteration matrix (involving the Jacobian) for the 2410C integration package DDEBDF. 2411C 2412C***SEE ALSO DDEBDF 2413C***ROUTINES CALLED DGBFA, DGEFA, DVNRMS 2414C***COMMON BLOCKS DDEBD1 2415C***REVISION HISTORY (YYMMDD) 2416C 820301 DATE WRITTEN 2417C 890531 Changed all specific intrinsics to generic. (WRB) 2418C 890911 Removed unnecessary intrinsics. (WRB) 2419C 891214 Prologue converted to Version 4.0 format. (BAB) 2420C 900328 Added TYPE section. (WRB) 2421C 910722 Updated AUTHOR section. (ALS) 2422C 920422 Changed DIMENSION statement. (WRB) 2423C***END PROLOGUE DPJAC 2424C 2425 INTEGER I, I1, I2, IER, II, IOWND, IOWNS, IPAR, IWM, J, J1, 2426 1 JJ, JSTART, KFLAG, L, LENP, MAXORD, MBA, MBAND, 2427 2 MEB1, MEBAND, METH, MITER, ML, ML3, MU, N, NEQ, 2428 3 NFE, NJE, NQ, NQU, NST, NYH 2429 DOUBLE PRECISION CON, DI, DVNRMS, EL0, EWT, 2430 1 FAC, FTEM, H, HL0, HMIN, HMXI, HU, R, R0, ROWND, ROWNS, 2431 2 RPAR, SAVF, SRUR, TN, UROUND, WM, Y, YH, YI, YJ, YJJ 2432 EXTERNAL DF, DJAC 2433 DIMENSION Y(*),YH(NYH,*),EWT(*),FTEM(*),SAVF(*),WM(*),IWM(*), 2434 1 RPAR(*),IPAR(*) 2435 COMMON /DDEBD1/ ROWND,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND, 2436 1 IOWND(14),IOWNS(6),IER,JSTART,KFLAG,L,METH,MITER, 2437 2 MAXORD,N,NQ,NST,NFE,NJE,NQU 2438C ------------------------------------------------------------------ 2439C DPJAC IS CALLED BY DSTOD TO COMPUTE AND PROCESS THE MATRIX 2440C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN. 2441C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE DJAC IF 2442C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5. 2443C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED. 2444C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN 2445C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION 2446C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE 2447C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5. 2448C 2449C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION 2450C WITH DPJAC USES THE FOLLOWING.. 2451C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY. 2452C FTEM = WORK ARRAY OF LENGTH N (ACOR IN DSTOD ). 2453C SAVF = ARRAY CONTAINING DF EVALUATED AT PREDICTED Y. 2454C WM = DOUBLE PRECISION WORK SPACE FOR MATRICES. ON OUTPUT IT 2455C CONTAINS THE 2456C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU 2457C DECOMPOSITION OF P IF MITER IS 1, 2 , 4, OR 5. 2458C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). 2459C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. 2460C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN 2461C INCREMENTS. WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 2462C 3. 2463C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING 2464C AT IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS 2465C THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER 2466C IS 4 OR 5. 2467C EL0 = EL(1) (INPUT). 2468C IER = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .NE. 0 IF 2469C P MATRIX FOUND TO BE SINGULAR. 2470C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND, 2471C MITER, N, NFE, AND NJE. 2472C----------------------------------------------------------------------- 2473C BEGIN BLOCK PERMITTING ...EXITS TO 240 2474C BEGIN BLOCK PERMITTING ...EXITS TO 220 2475C BEGIN BLOCK PERMITTING ...EXITS TO 130 2476C BEGIN BLOCK PERMITTING ...EXITS TO 70 2477C***FIRST EXECUTABLE STATEMENT DPJAC 2478 NJE = NJE + 1 2479 HL0 = H*EL0 2480 GO TO (10,40,90,140,170), MITER 2481C IF MITER = 1, CALL DJAC AND MULTIPLY BY SCALAR. 2482C ----------------------- 2483 10 CONTINUE 2484 LENP = N*N 2485 DO 20 I = 1, LENP 2486 WM(I+2) = 0.0D0 2487 20 CONTINUE 2488 CALL DJAC(TN,Y,WM(3),N,RPAR,IPAR) 2489 CON = -HL0 2490 DO 30 I = 1, LENP 2491 WM(I+2) = WM(I+2)*CON 2492 30 CONTINUE 2493C ...EXIT 2494 GO TO 70 2495C IF MITER = 2, MAKE N CALLS TO DF TO APPROXIMATE J. 2496C -------------------- 2497 40 CONTINUE 2498 FAC = DVNRMS(N,SAVF,EWT) 2499 R0 = 1000.0D0*ABS(H)*UROUND*N*FAC 2500 IF (R0 .EQ. 0.0D0) R0 = 1.0D0 2501 SRUR = WM(1) 2502 J1 = 2 2503 DO 60 J = 1, N 2504 YJ = Y(J) 2505 R = MAX(SRUR*ABS(YJ),R0*EWT(J)) 2506 Y(J) = Y(J) + R 2507 FAC = -HL0/R 2508 CALL DF(TN,Y,FTEM,RPAR,IPAR) 2509 DO 50 I = 1, N 2510 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC 2511 50 CONTINUE 2512 Y(J) = YJ 2513 J1 = J1 + N 2514 60 CONTINUE 2515 NFE = NFE + N 2516 70 CONTINUE 2517C ADD IDENTITY MATRIX. 2518C ------------------------------------------------- 2519 J = 3 2520 DO 80 I = 1, N 2521 WM(J) = WM(J) + 1.0D0 2522 J = J + (N + 1) 2523 80 CONTINUE 2524C DO LU DECOMPOSITION ON P. 2525C -------------------------------------------- 2526 CALL DGEFA(WM(3),N,N,IWM(21),IER) 2527C .........EXIT 2528 GO TO 240 2529C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND 2530C P. --------- 2531 90 CONTINUE 2532 WM(2) = HL0 2533 IER = 0 2534 R = EL0*0.1D0 2535 DO 100 I = 1, N 2536 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2)) 2537 100 CONTINUE 2538 CALL DF(TN,Y,WM(3),RPAR,IPAR) 2539 NFE = NFE + 1 2540 DO 120 I = 1, N 2541 R0 = H*SAVF(I) - YH(I,2) 2542 DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I)) 2543 WM(I+2) = 1.0D0 2544 IF (ABS(R0) .LT. UROUND*EWT(I)) GO TO 110 2545C .........EXIT 2546 IF (ABS(DI) .EQ. 0.0D0) GO TO 130 2547 WM(I+2) = 0.1D0*R0/DI 2548 110 CONTINUE 2549 120 CONTINUE 2550C .........EXIT 2551 GO TO 240 2552 130 CONTINUE 2553 IER = -1 2554C ......EXIT 2555 GO TO 240 2556C IF MITER = 4, CALL DJAC AND MULTIPLY BY SCALAR. 2557C ----------------------- 2558 140 CONTINUE 2559 ML = IWM(1) 2560 MU = IWM(2) 2561 ML3 = 3 2562 MBAND = ML + MU + 1 2563 MEBAND = MBAND + ML 2564 LENP = MEBAND*N 2565 DO 150 I = 1, LENP 2566 WM(I+2) = 0.0D0 2567 150 CONTINUE 2568 CALL DJAC(TN,Y,WM(ML3),MEBAND,RPAR,IPAR) 2569 CON = -HL0 2570 DO 160 I = 1, LENP 2571 WM(I+2) = WM(I+2)*CON 2572 160 CONTINUE 2573C ...EXIT 2574 GO TO 220 2575C IF MITER = 5, MAKE MBAND CALLS TO DF TO APPROXIMATE J. 2576C ---------------- 2577 170 CONTINUE 2578 ML = IWM(1) 2579 MU = IWM(2) 2580 MBAND = ML + MU + 1 2581 MBA = MIN(MBAND,N) 2582 MEBAND = MBAND + ML 2583 MEB1 = MEBAND - 1 2584 SRUR = WM(1) 2585 FAC = DVNRMS(N,SAVF,EWT) 2586 R0 = 1000.0D0*ABS(H)*UROUND*N*FAC 2587 IF (R0 .EQ. 0.0D0) R0 = 1.0D0 2588 DO 210 J = 1, MBA 2589 DO 180 I = J, N, MBAND 2590 YI = Y(I) 2591 R = MAX(SRUR*ABS(YI),R0*EWT(I)) 2592 Y(I) = Y(I) + R 2593 180 CONTINUE 2594 CALL DF(TN,Y,FTEM,RPAR,IPAR) 2595 DO 200 JJ = J, N, MBAND 2596 Y(JJ) = YH(JJ,1) 2597 YJJ = Y(JJ) 2598 R = MAX(SRUR*ABS(YJJ),R0*EWT(JJ)) 2599 FAC = -HL0/R 2600 I1 = MAX(JJ-MU,1) 2601 I2 = MIN(JJ+ML,N) 2602 II = JJ*MEB1 - ML + 2 2603 DO 190 I = I1, I2 2604 WM(II+I) = (FTEM(I) - SAVF(I))*FAC 2605 190 CONTINUE 2606 200 CONTINUE 2607 210 CONTINUE 2608 NFE = NFE + MBA 2609 220 CONTINUE 2610C ADD IDENTITY MATRIX. 2611C ------------------------------------------------- 2612 II = MBAND + 2 2613 DO 230 I = 1, N 2614 WM(II) = WM(II) + 1.0D0 2615 II = II + MEBAND 2616 230 CONTINUE 2617C DO LU DECOMPOSITION OF P. 2618C -------------------------------------------- 2619 CALL DGBFA(WM(3),MEBAND,N,ML,MU,IWM(21),IER) 2620 240 CONTINUE 2621 RETURN 2622C ----------------------- END OF SUBROUTINE DPJAC 2623C ----------------------- 2624 END 2625*DECK DSLVS 2626 SUBROUTINE DSLVS (WM, IWM, X, TEM) 2627C***BEGIN PROLOGUE DSLVS 2628C***SUBSIDIARY 2629C***PURPOSE Subsidiary to DDEBDF 2630C***LIBRARY SLATEC 2631C***TYPE DOUBLE PRECISION (SLVS-S, DSLVS-D) 2632C***AUTHOR Watts, H. A., (SNLA) 2633C***DESCRIPTION 2634C 2635C DSLVS solves the linear system in the iteration scheme for the 2636C integrator package DDEBDF. 2637C 2638C***SEE ALSO DDEBDF 2639C***ROUTINES CALLED DGBSL, DGESL 2640C***COMMON BLOCKS DDEBD1 2641C***REVISION HISTORY (YYMMDD) 2642C 820301 DATE WRITTEN 2643C 890531 Changed all specific intrinsics to generic. (WRB) 2644C 891214 Prologue converted to Version 4.0 format. (BAB) 2645C 900328 Added TYPE section. (WRB) 2646C 910722 Updated AUTHOR section. (ALS) 2647C 920422 Changed DIMENSION statement. (WRB) 2648C***END PROLOGUE DSLVS 2649C 2650 INTEGER I, IER, IOWND, IOWNS, IWM, JSTART, KFLAG, L, MAXORD, 2651 1 MEBAND, METH, MITER, ML, MU, N, NFE, NJE, NQ, NQU, NST 2652 DOUBLE PRECISION DI, EL0, H, HL0, HMIN, HMXI, HU, PHL0, 2653 1 R, ROWND, ROWNS, TEM, TN, UROUND, WM, X 2654 DIMENSION WM(*), IWM(*), X(*), TEM(*) 2655 COMMON /DDEBD1/ ROWND,ROWNS(210),EL0,H,HMIN,HMXI,HU,TN,UROUND, 2656 1 IOWND(14),IOWNS(6),IER,JSTART,KFLAG,L,METH,MITER, 2657 2 MAXORD,N,NQ,NST,NFE,NJE,NQU 2658C ------------------------------------------------------------------ 2659C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING 2660C FROM A CHORD ITERATION. IT IS CALLED BY DSTOD IF MITER .NE. 0. 2661C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS. 2662C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL 2663C MATRIX, AND THEN COMPUTES THE SOLUTION. 2664C IF MITER IS 4 OR 5, IT CALLS DGBSL. 2665C COMMUNICATION WITH DSLVS USES THE FOLLOWING VARIABLES.. 2666C WM = DOUBLE PRECISION WORK SPACE CONTAINING THE INVERSE DIAGONAL 2667C MATRIX IF MITER 2668C IS 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. 2669C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3). 2670C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA.. 2671C WM(1) = SQRT(UROUND) (NOT USED HERE), 2672C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 2673C 3. 2674C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING 2675C AT IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS 2676C THE BAND PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 2677C 4 OR 5. 2678C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION 2679C VECTOR ON OUTPUT, OF LENGTH N. 2680C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION. 2681C IER = OUTPUT FLAG (IN COMMON). IER = 0 IF NO TROUBLE OCCURRED. 2682C IER = -1 IF A SINGULAR MATRIX AROSE WITH MITER = 3. 2683C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N. 2684C----------------------------------------------------------------------- 2685C BEGIN BLOCK PERMITTING ...EXITS TO 80 2686C BEGIN BLOCK PERMITTING ...EXITS TO 60 2687C***FIRST EXECUTABLE STATEMENT DSLVS 2688 IER = 0 2689 GO TO (10,10,20,70,70), MITER 2690 10 CONTINUE 2691 CALL DGESL(WM(3),N,N,IWM(21),X,0) 2692C ......EXIT 2693 GO TO 80 2694C 2695 20 CONTINUE 2696 PHL0 = WM(2) 2697 HL0 = H*EL0 2698 WM(2) = HL0 2699 IF (HL0 .EQ. PHL0) GO TO 40 2700 R = HL0/PHL0 2701 DO 30 I = 1, N 2702 DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2)) 2703C .........EXIT 2704 IF (ABS(DI) .EQ. 0.0D0) GO TO 60 2705 WM(I+2) = 1.0D0/DI 2706 30 CONTINUE 2707 40 CONTINUE 2708 DO 50 I = 1, N 2709 X(I) = WM(I+2)*X(I) 2710 50 CONTINUE 2711C ......EXIT 2712 GO TO 80 2713 60 CONTINUE 2714 IER = -1 2715C ...EXIT 2716 GO TO 80 2717C 2718 70 CONTINUE 2719 ML = IWM(1) 2720 MU = IWM(2) 2721 MEBAND = 2*ML + MU + 1 2722 CALL DGBSL(WM(3),MEBAND,N,ML,MU,IWM(21),X,0) 2723 80 CONTINUE 2724 RETURN 2725C ----------------------- END OF SUBROUTINE DSLVS 2726C ----------------------- 2727 END 2728*DECK DGBFA 2729 SUBROUTINE DGBFA (ABD, LDA, N, ML, MU, IPVT, INFO) 2730C***BEGIN PROLOGUE DGBFA 2731C***PURPOSE Factor a band matrix using Gaussian elimination. 2732C***LIBRARY SLATEC (LINPACK) 2733C***CATEGORY D2A2 2734C***TYPE DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C) 2735C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION 2736C***AUTHOR Moler, C. B., (U. of New Mexico) 2737C***DESCRIPTION 2738C 2739C DGBFA factors a double precision band matrix by elimination. 2740C 2741C DGBFA is usually called by DGBCO, but it can be called 2742C directly with a saving in time if RCOND is not needed. 2743C 2744C On Entry 2745C 2746C ABD DOUBLE PRECISION(LDA, N) 2747C contains the matrix in band storage. The columns 2748C of the matrix are stored in the columns of ABD and 2749C the diagonals of the matrix are stored in rows 2750C ML+1 through 2*ML+MU+1 of ABD . 2751C See the comments below for details. 2752C 2753C LDA INTEGER 2754C the leading dimension of the array ABD . 2755C LDA must be .GE. 2*ML + MU + 1 . 2756C 2757C N INTEGER 2758C the order of the original matrix. 2759C 2760C ML INTEGER 2761C number of diagonals below the main diagonal. 2762C 0 .LE. ML .LT. N . 2763C 2764C MU INTEGER 2765C number of diagonals above the main diagonal. 2766C 0 .LE. MU .LT. N . 2767C More efficient if ML .LE. MU . 2768C On Return 2769C 2770C ABD an upper triangular matrix in band storage and 2771C the multipliers which were used to obtain it. 2772C The factorization can be written A = L*U where 2773C L is a product of permutation and unit lower 2774C triangular matrices and U is upper triangular. 2775C 2776C IPVT INTEGER(N) 2777C an integer vector of pivot indices. 2778C 2779C INFO INTEGER 2780C = 0 normal value. 2781C = K if U(K,K) .EQ. 0.0 . This is not an error 2782C condition for this subroutine, but it does 2783C indicate that DGBSL will divide by zero if 2784C called. Use RCOND in DGBCO for a reliable 2785C indication of singularity. 2786C 2787C Band Storage 2788C 2789C If A is a band matrix, the following program segment 2790C will set up the input. 2791C 2792C ML = (band width below the diagonal) 2793C MU = (band width above the diagonal) 2794C M = ML + MU + 1 2795C DO 20 J = 1, N 2796C I1 = MAX(1, J-MU) 2797C I2 = MIN(N, J+ML) 2798C DO 10 I = I1, I2 2799C K = I - J + M 2800C ABD(K,J) = A(I,J) 2801C 10 CONTINUE 2802C 20 CONTINUE 2803C 2804C This uses rows ML+1 through 2*ML+MU+1 of ABD . 2805C In addition, the first ML rows in ABD are used for 2806C elements generated during the triangularization. 2807C The total number of rows needed in ABD is 2*ML+MU+1 . 2808C The ML+MU by ML+MU upper left triangle and the 2809C ML by ML lower right triangle are not referenced. 2810C 2811C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. 2812C Stewart, LINPACK Users' Guide, SIAM, 1979. 2813C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX 2814C***REVISION HISTORY (YYMMDD) 2815C 780814 DATE WRITTEN 2816C 890531 Changed all specific intrinsics to generic. (WRB) 2817C 890831 Modified array declarations. (WRB) 2818C 890831 REVISION DATE from Version 3.2 2819C 891214 Prologue converted to Version 4.0 format. (BAB) 2820C 900326 Removed duplicate information from DESCRIPTION section. 2821C (WRB) 2822C 920501 Reformatted the REFERENCES section. (WRB) 2823C***END PROLOGUE DGBFA 2824 INTEGER LDA,N,ML,MU,IPVT(*),INFO 2825 DOUBLE PRECISION ABD(LDA,*) 2826C 2827 DOUBLE PRECISION T 2828 INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1 2829C 2830C***FIRST EXECUTABLE STATEMENT DGBFA 2831 M = ML + MU + 1 2832 INFO = 0 2833C 2834C ZERO INITIAL FILL-IN COLUMNS 2835C 2836 J0 = MU + 2 2837 J1 = MIN(N,M) - 1 2838 IF (J1 .LT. J0) GO TO 30 2839 DO 20 JZ = J0, J1 2840 I0 = M + 1 - JZ 2841 DO 10 I = I0, ML 2842 ABD(I,JZ) = 0.0D0 2843 10 CONTINUE 2844 20 CONTINUE 2845 30 CONTINUE 2846 JZ = J1 2847 JU = 0 2848C 2849C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 2850C 2851 NM1 = N - 1 2852 IF (NM1 .LT. 1) GO TO 130 2853 DO 120 K = 1, NM1 2854 KP1 = K + 1 2855C 2856C ZERO NEXT FILL-IN COLUMN 2857C 2858 JZ = JZ + 1 2859 IF (JZ .GT. N) GO TO 50 2860 IF (ML .LT. 1) GO TO 50 2861 DO 40 I = 1, ML 2862 ABD(I,JZ) = 0.0D0 2863 40 CONTINUE 2864 50 CONTINUE 2865C 2866C FIND L = PIVOT INDEX 2867C 2868 LM = MIN(ML,N-K) 2869 L = IDAMAX(LM+1,ABD(M,K),1) + M - 1 2870 IPVT(K) = L + K - M 2871C 2872C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED 2873C 2874 IF (ABD(L,K) .EQ. 0.0D0) GO TO 100 2875C 2876C INTERCHANGE IF NECESSARY 2877C 2878 IF (L .EQ. M) GO TO 60 2879 T = ABD(L,K) 2880 ABD(L,K) = ABD(M,K) 2881 ABD(M,K) = T 2882 60 CONTINUE 2883C 2884C COMPUTE MULTIPLIERS 2885C 2886 T = -1.0D0/ABD(M,K) 2887 CALL DSCAL(LM,T,ABD(M+1,K),1) 2888C 2889C ROW ELIMINATION WITH COLUMN INDEXING 2890C 2891 JU = MIN(MAX(JU,MU+IPVT(K)),N) 2892 MM = M 2893 IF (JU .LT. KP1) GO TO 90 2894 DO 80 J = KP1, JU 2895 L = L - 1 2896 MM = MM - 1 2897 T = ABD(L,J) 2898 IF (L .EQ. MM) GO TO 70 2899 ABD(L,J) = ABD(MM,J) 2900 ABD(MM,J) = T 2901 70 CONTINUE 2902 CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1) 2903 80 CONTINUE 2904 90 CONTINUE 2905 GO TO 110 2906 100 CONTINUE 2907 INFO = K 2908 110 CONTINUE 2909 120 CONTINUE 2910 130 CONTINUE 2911 IPVT(N) = N 2912 IF (ABD(M,N) .EQ. 0.0D0) INFO = N 2913 RETURN 2914 END 2915*DECK DGBSL 2916 SUBROUTINE DGBSL (ABD, LDA, N, ML, MU, IPVT, B, JOB) 2917C***BEGIN PROLOGUE DGBSL 2918C***PURPOSE Solve the real band system A*X=B or TRANS(A)*X=B using 2919C the factors computed by DGBCO or DGBFA. 2920C***LIBRARY SLATEC (LINPACK) 2921C***CATEGORY D2A2 2922C***TYPE DOUBLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C) 2923C***KEYWORDS BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE 2924C***AUTHOR Moler, C. B., (U. of New Mexico) 2925C***DESCRIPTION 2926C 2927C DGBSL solves the double precision band system 2928C A * X = B or TRANS(A) * X = B 2929C using the factors computed by DGBCO or DGBFA. 2930C 2931C On Entry 2932C 2933C ABD DOUBLE PRECISION(LDA, N) 2934C the output from DGBCO or DGBFA. 2935C 2936C LDA INTEGER 2937C the leading dimension of the array ABD . 2938C 2939C N INTEGER 2940C the order of the original matrix. 2941C 2942C ML INTEGER 2943C number of diagonals below the main diagonal. 2944C 2945C MU INTEGER 2946C number of diagonals above the main diagonal. 2947C 2948C IPVT INTEGER(N) 2949C the pivot vector from DGBCO or DGBFA. 2950C 2951C B DOUBLE PRECISION(N) 2952C the right hand side vector. 2953C 2954C JOB INTEGER 2955C = 0 to solve A*X = B , 2956C = nonzero to solve TRANS(A)*X = B , where 2957C TRANS(A) is the transpose. 2958C 2959C On Return 2960C 2961C B the solution vector X . 2962C 2963C Error Condition 2964C 2965C A division by zero will occur if the input factor contains a 2966C zero on the diagonal. Technically this indicates singularity 2967C but it is often caused by improper arguments or improper 2968C setting of LDA . It will not occur if the subroutines are 2969C called correctly and if DGBCO has set RCOND .GT. 0.0 2970C or DGBFA has set INFO .EQ. 0 . 2971C 2972C To compute INVERSE(A) * C where C is a matrix 2973C with P columns 2974C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z) 2975C IF (RCOND is too small) GO TO ... 2976C DO 10 J = 1, P 2977C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0) 2978C 10 CONTINUE 2979C 2980C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. 2981C Stewart, LINPACK Users' Guide, SIAM, 1979. 2982C***ROUTINES CALLED DAXPY, DDOT 2983C***REVISION HISTORY (YYMMDD) 2984C 780814 DATE WRITTEN 2985C 890531 Changed all specific intrinsics to generic. (WRB) 2986C 890831 Modified array declarations. (WRB) 2987C 890831 REVISION DATE from Version 3.2 2988C 891214 Prologue converted to Version 4.0 format. (BAB) 2989C 900326 Removed duplicate information from DESCRIPTION section. 2990C (WRB) 2991C 920501 Reformatted the REFERENCES section. (WRB) 2992C***END PROLOGUE DGBSL 2993 INTEGER LDA,N,ML,MU,IPVT(*),JOB 2994 DOUBLE PRECISION ABD(LDA,*),B(*) 2995C 2996 DOUBLE PRECISION DDOT,T 2997 INTEGER K,KB,L,LA,LB,LM,M,NM1 2998C***FIRST EXECUTABLE STATEMENT DGBSL 2999 M = MU + ML + 1 3000 NM1 = N - 1 3001 IF (JOB .NE. 0) GO TO 50 3002C 3003C JOB = 0 , SOLVE A * X = B 3004C FIRST SOLVE L*Y = B 3005C 3006 IF (ML .EQ. 0) GO TO 30 3007 IF (NM1 .LT. 1) GO TO 30 3008 DO 20 K = 1, NM1 3009 LM = MIN(ML,N-K) 3010 L = IPVT(K) 3011 T = B(L) 3012 IF (L .EQ. K) GO TO 10 3013 B(L) = B(K) 3014 B(K) = T 3015 10 CONTINUE 3016 CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1) 3017 20 CONTINUE 3018 30 CONTINUE 3019C 3020C NOW SOLVE U*X = Y 3021C 3022 DO 40 KB = 1, N 3023 K = N + 1 - KB 3024 B(K) = B(K)/ABD(M,K) 3025 LM = MIN(K,M) - 1 3026 LA = M - LM 3027 LB = K - LM 3028 T = -B(K) 3029 CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1) 3030 40 CONTINUE 3031 GO TO 100 3032 50 CONTINUE 3033C 3034C JOB = NONZERO, SOLVE TRANS(A) * X = B 3035C FIRST SOLVE TRANS(U)*Y = B 3036C 3037 DO 60 K = 1, N 3038 LM = MIN(K,M) - 1 3039 LA = M - LM 3040 LB = K - LM 3041 T = DDOT(LM,ABD(LA,K),1,B(LB),1) 3042 B(K) = (B(K) - T)/ABD(M,K) 3043 60 CONTINUE 3044C 3045C NOW SOLVE TRANS(L)*X = Y 3046C 3047 IF (ML .EQ. 0) GO TO 90 3048 IF (NM1 .LT. 1) GO TO 90 3049 DO 80 KB = 1, NM1 3050 K = N - KB 3051 LM = MIN(ML,N-K) 3052 B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1) 3053 L = IPVT(K) 3054 IF (L .EQ. K) GO TO 70 3055 T = B(L) 3056 B(L) = B(K) 3057 B(K) = T 3058 70 CONTINUE 3059 80 CONTINUE 3060 90 CONTINUE 3061 100 CONTINUE 3062 RETURN 3063 END 3064*DECK DGEFA 3065 SUBROUTINE DGEFA (A, LDA, N, IPVT, INFO) 3066C***BEGIN PROLOGUE DGEFA 3067C***PURPOSE Factor a matrix using Gaussian elimination. 3068C***LIBRARY SLATEC (LINPACK) 3069C***CATEGORY D2A1 3070C***TYPE DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C) 3071C***KEYWORDS GENERAL MATRIX, LINEAR ALGEBRA, LINPACK, 3072C MATRIX FACTORIZATION 3073C***AUTHOR Moler, C. B., (U. of New Mexico) 3074C***DESCRIPTION 3075C 3076C DGEFA factors a double precision matrix by Gaussian elimination. 3077C 3078C DGEFA is usually called by DGECO, but it can be called 3079C directly with a saving in time if RCOND is not needed. 3080C (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) . 3081C 3082C On Entry 3083C 3084C A DOUBLE PRECISION(LDA, N) 3085C the matrix to be factored. 3086C 3087C LDA INTEGER 3088C the leading dimension of the array A . 3089C 3090C N INTEGER 3091C the order of the matrix A . 3092C 3093C On Return 3094C 3095C A an upper triangular matrix and the multipliers 3096C which were used to obtain it. 3097C The factorization can be written A = L*U where 3098C L is a product of permutation and unit lower 3099C triangular matrices and U is upper triangular. 3100C 3101C IPVT INTEGER(N) 3102C an integer vector of pivot indices. 3103C 3104C INFO INTEGER 3105C = 0 normal value. 3106C = K if U(K,K) .EQ. 0.0 . This is not an error 3107C condition for this subroutine, but it does 3108C indicate that DGESL or DGEDI will divide by zero 3109C if called. Use RCOND in DGECO for a reliable 3110C indication of singularity. 3111C 3112C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. 3113C Stewart, LINPACK Users' Guide, SIAM, 1979. 3114C***ROUTINES CALLED DAXPY, DSCAL, IDAMAX 3115C***REVISION HISTORY (YYMMDD) 3116C 780814 DATE WRITTEN 3117C 890831 Modified array declarations. (WRB) 3118C 890831 REVISION DATE from Version 3.2 3119C 891214 Prologue converted to Version 4.0 format. (BAB) 3120C 900326 Removed duplicate information from DESCRIPTION section. 3121C (WRB) 3122C 920501 Reformatted the REFERENCES section. (WRB) 3123C***END PROLOGUE DGEFA 3124 INTEGER LDA,N,IPVT(*),INFO 3125 DOUBLE PRECISION A(LDA,*) 3126C 3127 DOUBLE PRECISION T 3128 INTEGER IDAMAX,J,K,KP1,L,NM1 3129C 3130C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 3131C 3132C***FIRST EXECUTABLE STATEMENT DGEFA 3133 INFO = 0 3134 NM1 = N - 1 3135 IF (NM1 .LT. 1) GO TO 70 3136 DO 60 K = 1, NM1 3137 KP1 = K + 1 3138C 3139C FIND L = PIVOT INDEX 3140C 3141 L = IDAMAX(N-K+1,A(K,K),1) + K - 1 3142 IPVT(K) = L 3143C 3144C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED 3145C 3146 IF (A(L,K) .EQ. 0.0D0) GO TO 40 3147C 3148C INTERCHANGE IF NECESSARY 3149C 3150 IF (L .EQ. K) GO TO 10 3151 T = A(L,K) 3152 A(L,K) = A(K,K) 3153 A(K,K) = T 3154 10 CONTINUE 3155C 3156C COMPUTE MULTIPLIERS 3157C 3158 T = -1.0D0/A(K,K) 3159 CALL DSCAL(N-K,T,A(K+1,K),1) 3160C 3161C ROW ELIMINATION WITH COLUMN INDEXING 3162C 3163 DO 30 J = KP1, N 3164 T = A(L,J) 3165 IF (L .EQ. K) GO TO 20 3166 A(L,J) = A(K,J) 3167 A(K,J) = T 3168 20 CONTINUE 3169 CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 3170 30 CONTINUE 3171 GO TO 50 3172 40 CONTINUE 3173 INFO = K 3174 50 CONTINUE 3175 60 CONTINUE 3176 70 CONTINUE 3177 IPVT(N) = N 3178 IF (A(N,N) .EQ. 0.0D0) INFO = N 3179 RETURN 3180 END 3181*DECK DGESL 3182 SUBROUTINE DGESL (A, LDA, N, IPVT, B, JOB) 3183C***BEGIN PROLOGUE DGESL 3184C***PURPOSE Solve the real system A*X=B or TRANS(A)*X=B using the 3185C factors computed by DGECO or DGEFA. 3186C***LIBRARY SLATEC (LINPACK) 3187C***CATEGORY D2A1 3188C***TYPE DOUBLE PRECISION (SGESL-S, DGESL-D, CGESL-C) 3189C***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE 3190C***AUTHOR Moler, C. B., (U. of New Mexico) 3191C***DESCRIPTION 3192C 3193C DGESL solves the double precision system 3194C A * X = B or TRANS(A) * X = B 3195C using the factors computed by DGECO or DGEFA. 3196C 3197C On Entry 3198C 3199C A DOUBLE PRECISION(LDA, N) 3200C the output from DGECO or DGEFA. 3201C 3202C LDA INTEGER 3203C the leading dimension of the array A . 3204C 3205C N INTEGER 3206C the order of the matrix A . 3207C 3208C IPVT INTEGER(N) 3209C the pivot vector from DGECO or DGEFA. 3210C 3211C B DOUBLE PRECISION(N) 3212C the right hand side vector. 3213C 3214C JOB INTEGER 3215C = 0 to solve A*X = B , 3216C = nonzero to solve TRANS(A)*X = B where 3217C TRANS(A) is the transpose. 3218C 3219C On Return 3220C 3221C B the solution vector X . 3222C 3223C Error Condition 3224C 3225C A division by zero will occur if the input factor contains a 3226C zero on the diagonal. Technically this indicates singularity 3227C but it is often caused by improper arguments or improper 3228C setting of LDA . It will not occur if the subroutines are 3229C called correctly and if DGECO has set RCOND .GT. 0.0 3230C or DGEFA has set INFO .EQ. 0 . 3231C 3232C To compute INVERSE(A) * C where C is a matrix 3233C with P columns 3234C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) 3235C IF (RCOND is too small) GO TO ... 3236C DO 10 J = 1, P 3237C CALL DGESL(A,LDA,N,IPVT,C(1,J),0) 3238C 10 CONTINUE 3239C 3240C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. 3241C Stewart, LINPACK Users' Guide, SIAM, 1979. 3242C***ROUTINES CALLED DAXPY, DDOT 3243C***REVISION HISTORY (YYMMDD) 3244C 780814 DATE WRITTEN 3245C 890831 Modified array declarations. (WRB) 3246C 890831 REVISION DATE from Version 3.2 3247C 891214 Prologue converted to Version 4.0 format. (BAB) 3248C 900326 Removed duplicate information from DESCRIPTION section. 3249C (WRB) 3250C 920501 Reformatted the REFERENCES section. (WRB) 3251C***END PROLOGUE DGESL 3252 INTEGER LDA,N,IPVT(*),JOB 3253 DOUBLE PRECISION A(LDA,*),B(*) 3254C 3255 DOUBLE PRECISION DDOT,T 3256 INTEGER K,KB,L,NM1 3257C***FIRST EXECUTABLE STATEMENT DGESL 3258 NM1 = N - 1 3259 IF (JOB .NE. 0) GO TO 50 3260C 3261C JOB = 0 , SOLVE A * X = B 3262C FIRST SOLVE L*Y = B 3263C 3264 IF (NM1 .LT. 1) GO TO 30 3265 DO 20 K = 1, NM1 3266 L = IPVT(K) 3267 T = B(L) 3268 IF (L .EQ. K) GO TO 10 3269 B(L) = B(K) 3270 B(K) = T 3271 10 CONTINUE 3272 CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 3273 20 CONTINUE 3274 30 CONTINUE 3275C 3276C NOW SOLVE U*X = Y 3277C 3278 DO 40 KB = 1, N 3279 K = N + 1 - KB 3280 B(K) = B(K)/A(K,K) 3281 T = -B(K) 3282 CALL DAXPY(K-1,T,A(1,K),1,B(1),1) 3283 40 CONTINUE 3284 GO TO 100 3285 50 CONTINUE 3286C 3287C JOB = NONZERO, SOLVE TRANS(A) * X = B 3288C FIRST SOLVE TRANS(U)*Y = B 3289C 3290 DO 60 K = 1, N 3291 T = DDOT(K-1,A(1,K),1,B(1),1) 3292 B(K) = (B(K) - T)/A(K,K) 3293 60 CONTINUE 3294C 3295C NOW SOLVE TRANS(L)*X = Y 3296C 3297 IF (NM1 .LT. 1) GO TO 90 3298 DO 80 KB = 1, NM1 3299 K = N - KB 3300 B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) 3301 L = IPVT(K) 3302 IF (L .EQ. K) GO TO 70 3303 T = B(L) 3304 B(L) = B(K) 3305 B(K) = T 3306 70 CONTINUE 3307 80 CONTINUE 3308 90 CONTINUE 3309 100 CONTINUE 3310 RETURN 3311 END 3312