1! 2! SLATEC: public domain 3! 4*DECK DDEABM 5 SUBROUTINE DDEABM (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, 6 + RWORK, LRW, IWORK, LIW, RPAR, IPAR) 7C***BEGIN PROLOGUE DDEABM 8C***PURPOSE Solve an initial value problem in ordinary differential 9C equations using an Adams-Bashforth method. 10C***LIBRARY SLATEC (DEPAC) 11C***CATEGORY I1A1B 12C***TYPE DOUBLE PRECISION (DEABM-S, DDEABM-D) 13C***KEYWORDS ADAMS-BASHFORTH METHOD, DEPAC, INITIAL VALUE PROBLEMS, 14C ODE, ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR 15C***AUTHOR Shampine, L. F., (SNLA) 16C Watts, H. A., (SNLA) 17C***DESCRIPTION 18C 19C This is the Adams code in the package of differential equation 20C solvers DEPAC, consisting of the codes DDERKF, DDEABM, and DDEBDF. 21C Design of the package was by L. F. Shampine and H. A. Watts. 22C It is documented in 23C SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE 24C Solvers. 25C DDEABM is a driver for a modification of the code ODE written by 26C L. F. Shampine and M. K. Gordon 27C Sandia Laboratories 28C Albuquerque, New Mexico 87185 29C 30C ********************************************************************** 31C * ABSTRACT * 32C ************ 33C 34C Subroutine DDEABM uses the Adams-Bashforth-Moulton 35C Predictor-Corrector formulas of orders one through twelve to 36C integrate a system of NEQ first order ordinary differential 37C equations of the form 38C DU/DX = DF(X,U) 39C when the vector Y(*) of initial values for U(*) at X=T is given. 40C The subroutine integrates from T to TOUT. It is easy to continue the 41C integration to get results at additional TOUT. This is the interval 42C mode of operation. It is also easy for the routine to return with 43C the solution at each intermediate step on the way to TOUT. This is 44C the intermediate-output mode of operation. 45C 46C DDEABM uses subprograms DDES, DSTEPS, DINTP, DHSTRT, DHVNRM, 47C D1MACH, and the error handling routine XERMSG. The only machine 48C dependent parameters to be assigned appear in D1MACH. 49C 50C ********************************************************************** 51C * Description of The Arguments To DDEABM (An Overview) * 52C ********************************************************************** 53C 54C The Parameters are 55C 56C DF -- This is the name of a subroutine which you provide to 57C define the differential equations. 58C 59C NEQ -- This is the number of (first order) differential 60C equations to be integrated. 61C 62C T -- This is a DOUBLE PRECISION value of the independent 63C variable. 64C 65C Y(*) -- This DOUBLE PRECISION array contains the solution 66C components at T. 67C 68C TOUT -- This is a DOUBLE PRECISION point at which a solution is 69C desired. 70C 71C INFO(*) -- The basic task of the code is to integrate the 72C differential equations from T to TOUT and return an 73C answer at TOUT. INFO(*) is an INTEGER array which is used 74C to communicate exactly how you want this task to be 75C carried out. 76C 77C RTOL, ATOL -- These DOUBLE PRECISION quantities represent 78C relative and absolute error tolerances which you 79C provide to indicate how accurately you wish the 80C solution to be computed. You may choose them to be 81C both scalars or else both vectors. 82C 83C IDID -- This scalar quantity is an indicator reporting what 84C the code did. You must monitor this INTEGER variable to 85C decide what action to take next. 86C 87C RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of 88C length LRW which provides the code with needed storage 89C space. 90C 91C IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW 92C which provides the code with needed storage space and an 93C across call flag. 94C 95C RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter 96C arrays which you can use for communication between your 97C calling program and the DF subroutine. 98C 99C Quantities which are used as input items are 100C NEQ, T, Y(*), TOUT, INFO(*), 101C RTOL, ATOL, RWORK(1), LRW and LIW. 102C 103C Quantities which may be altered by the code are 104C T, Y(*), INFO(1), RTOL, ATOL, 105C IDID, RWORK(*) and IWORK(*). 106C 107C ********************************************************************** 108C * INPUT -- What To Do On The First Call To DDEABM * 109C ********************************************************************** 110C 111C The first call of the code is defined to be the start of each new 112C problem. Read through the descriptions of all the following items, 113C provide sufficient storage space for designated arrays, set 114C appropriate variables for the initialization of the problem, and 115C give information about how you want the problem to be solved. 116C 117C 118C DF -- Provide a subroutine of the form 119C DF(X,U,UPRIME,RPAR,IPAR) 120C to define the system of first order differential equations 121C which is to be solved. For the given values of X and the 122C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must 123C evaluate the NEQ components of the system of differential 124C equations DU/DX=DF(X,U) and store the derivatives in the 125C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for 126C equations I=1,...,NEQ. 127C 128C Subroutine DF must NOT alter X or U(*). You must declare 129C the name df in an external statement in your program that 130C calls DDEABM. You must dimension U and UPRIME in DF. 131C 132C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter 133C arrays which you can use for communication between your 134C calling program and subroutine DF. They are not used or 135C altered by DDEABM. If you do not need RPAR or IPAR, 136C ignore these parameters by treating them as dummy 137C arguments. If you do choose to use them, dimension them in 138C your calling program and in DF as arrays of appropriate 139C length. 140C 141C NEQ -- Set it to the number of differential equations. 142C (NEQ .GE. 1) 143C 144C T -- Set it to the initial point of the integration. 145C You must use a program variable for T because the code 146C changes its value. 147C 148C Y(*) -- Set this vector to the initial values of the NEQ solution 149C components at the initial point. You must dimension Y at 150C least NEQ in your calling program. 151C 152C TOUT -- Set it to the first point at which a solution 153C is desired. You can take TOUT = T, in which case the code 154C will evaluate the derivative of the solution at T and 155C return. Integration either forward in T (TOUT .GT. T) or 156C backward in T (TOUT .LT. T) is permitted. 157C 158C The code advances the solution from T to TOUT using 159C step sizes which are automatically selected so as to 160C achieve the desired accuracy. If you wish, the code will 161C return with the solution and its derivative following 162C each intermediate step (intermediate-output mode) so that 163C you can monitor them, but you still must provide TOUT in 164C accord with the basic aim of the code. 165C 166C The first step taken by the code is a critical one 167C because it must reflect how fast the solution changes near 168C the initial point. The code automatically selects an 169C initial step size which is practically always suitable for 170C the problem. By using the fact that the code will not step 171C past TOUT in the first step, you could, if necessary, 172C restrict the length of the initial step size. 173C 174C For some problems it may not be permissible to integrate 175C past a point TSTOP because a discontinuity occurs there 176C or the solution or its derivative is not defined beyond 177C TSTOP. When you have declared a TSTOP point (see INFO(4) 178C and RWORK(1)), you have told the code not to integrate 179C past TSTOP. In this case any TOUT beyond TSTOP is invalid 180C input. 181C 182C INFO(*) -- Use the INFO array to give the code more details about 183C how you want your problem solved. This array should be 184C dimensioned of length 15 to accommodate other members of 185C DEPAC or possible future extensions, though DDEABM uses 186C only the first four entries. You must respond to all of 187C the following items which are arranged as questions. The 188C simplest use of the code corresponds to answering all 189C questions as YES ,i.e. setting ALL entries of INFO to 0. 190C 191C INFO(1) -- This parameter enables the code to initialize 192C itself. You must set it to indicate the start of every 193C new problem. 194C 195C **** Is this the first call for this problem ... 196C YES -- set INFO(1) = 0 197C NO -- not applicable here. 198C See below for continuation calls. **** 199C 200C INFO(2) -- How much accuracy you want of your solution 201C is specified by the error tolerances RTOL and ATOL. 202C The simplest use is to take them both to be scalars. 203C To obtain more flexibility, they can both be vectors. 204C The code must be told your choice. 205C 206C **** Are both error tolerances RTOL, ATOL scalars ... 207C YES -- set INFO(2) = 0 208C and input scalars for both RTOL and ATOL 209C NO -- set INFO(2) = 1 210C and input arrays for both RTOL and ATOL **** 211C 212C INFO(3) -- The code integrates from T in the direction 213C of TOUT by steps. If you wish, it will return the 214C computed solution and derivative at the next 215C intermediate step (the intermediate-output mode) or 216C TOUT, whichever comes first. This is a good way to 217C proceed if you want to see the behavior of the solution. 218C If you must have solutions at a great many specific 219C TOUT points, this code will compute them efficiently. 220C 221C **** Do you want the solution only at 222C TOUT (and not at the next intermediate step) ... 223C YES -- set INFO(3) = 0 224C NO -- set INFO(3) = 1 **** 225C 226C INFO(4) -- To handle solutions at a great many specific 227C values TOUT efficiently, this code may integrate past 228C TOUT and interpolate to obtain the result at TOUT. 229C Sometimes it is not possible to integrate beyond some 230C point TSTOP because the equation changes there or it is 231C not defined past TSTOP. Then you must tell the code 232C not to go past. 233C 234C **** Can the integration be carried out without any 235C Restrictions on the independent variable T ... 236C YES -- set INFO(4)=0 237C NO -- set INFO(4)=1 238C and define the stopping point TSTOP by 239C setting RWORK(1)=TSTOP **** 240C 241C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL) 242C error tolerances to tell the code how accurately you want 243C the solution to be computed. They must be defined as 244C program variables because the code may change them. You 245C have two choices -- 246C Both RTOL and ATOL are scalars. (INFO(2)=0) 247C Both RTOL and ATOL are vectors. (INFO(2)=1) 248C In either case all components must be non-negative. 249C 250C The tolerances are used by the code in a local error test 251C at each step which requires roughly that 252C ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL 253C for each vector component. 254C (More specifically, a Euclidean norm is used to measure 255C the size of vectors, and the error test uses the magnitude 256C of the solution at the beginning of the step.) 257C 258C The true (global) error is the difference between the true 259C solution of the initial value problem and the computed 260C approximation. Practically all present day codes, 261C including this one, control the local error at each step 262C and do not even attempt to control the global error 263C directly. Roughly speaking, they produce a solution Y(T) 264C which satisfies the differential equations with a 265C residual R(T), DY(T)/DT = DF(T,Y(T)) + R(T) , 266C and, almost always, R(T) is bounded by the error 267C tolerances. Usually, but not always, the true accuracy of 268C the computed Y is comparable to the error tolerances. This 269C code will usually, but not always, deliver a more accurate 270C solution if you reduce the tolerances and integrate again. 271C By comparing two such solutions you can get a fairly 272C reliable idea of the true error in the solution at the 273C bigger tolerances. 274C 275C Setting ATOL=0.D0 results in a pure relative error test on 276C that component. Setting RTOL=0. results in a pure absolute 277C error test on that component. A mixed test with non-zero 278C RTOL and ATOL corresponds roughly to a relative error 279C test when the solution component is much bigger than ATOL 280C and to an absolute error test when the solution component 281C is smaller than the threshold ATOL. 282C 283C Proper selection of the absolute error control parameters 284C ATOL requires you to have some idea of the scale of the 285C solution components. To acquire this information may mean 286C that you will have to solve the problem more than once. In 287C the absence of scale information, you should ask for some 288C relative accuracy in all the components (by setting RTOL 289C values non-zero) and perhaps impose extremely small 290C absolute error tolerances to protect against the danger of 291C a solution component becoming zero. 292C 293C The code will not attempt to compute a solution at an 294C accuracy unreasonable for the machine being used. It will 295C advise you if you ask for too much accuracy and inform 296C you as to the maximum accuracy it believes possible. 297C 298C RWORK(*) -- Dimension this DOUBLE PRECISION work array of length 299C LRW in your calling program. 300C 301C RWORK(1) -- If you have set INFO(4)=0, you can ignore this 302C optional input parameter. Otherwise you must define a 303C stopping point TSTOP by setting RWORK(1) = TSTOP. 304C (for some problems it may not be permissible to integrate 305C past a point TSTOP because a discontinuity occurs there 306C or the solution or its derivative is not defined beyond 307C TSTOP.) 308C 309C LRW -- Set it to the declared length of the RWORK array. 310C You must have LRW .GE. 130+21*NEQ 311C 312C IWORK(*) -- Dimension this INTEGER work array of length LIW in 313C your calling program. 314C 315C LIW -- Set it to the declared length of the IWORK array. 316C You must have LIW .GE. 51 317C 318C RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and 319C INTEGER type, respectively. You can use them for 320C communication between your program that calls DDEABM and 321C the DF subroutine. They are not used or altered by 322C DDEABM. If you do not need RPAR or IPAR, ignore these 323C parameters by treating them as dummy arguments. If you do 324C choose to use them, dimension them in your calling program 325C and in DF as arrays of appropriate length. 326C 327C ********************************************************************** 328C * OUTPUT -- After Any Return From DDEABM * 329C ********************************************************************** 330C 331C The principal aim of the code is to return a computed solution at 332C TOUT, although it is also possible to obtain intermediate results 333C along the way. To find out whether the code achieved its goal 334C or if the integration process was interrupted before the task was 335C completed, you must check the IDID parameter. 336C 337C 338C T -- The solution was successfully advanced to the 339C output value of T. 340C 341C Y(*) -- Contains the computed solution approximation at T. 342C You may also be interested in the approximate derivative 343C of the solution at T. It is contained in 344C RWORK(21),...,RWORK(20+NEQ). 345C 346C IDID -- Reports what the code did 347C 348C *** Task Completed *** 349C Reported by positive values of IDID 350C 351C IDID = 1 -- A step was successfully taken in the 352C intermediate-output mode. The code has not 353C yet reached TOUT. 354C 355C IDID = 2 -- The integration to TOUT was successfully 356C completed (T=TOUT) by stepping exactly to TOUT. 357C 358C IDID = 3 -- The integration to TOUT was successfully 359C completed (T=TOUT) by stepping past TOUT. 360C Y(*) is obtained by interpolation. 361C 362C *** Task Interrupted *** 363C Reported by negative values of IDID 364C 365C IDID = -1 -- A large amount of work has been expended. 366C (500 steps attempted) 367C 368C IDID = -2 -- The error tolerances are too stringent. 369C 370C IDID = -3 -- The local error test cannot be satisfied 371C because you specified a zero component in ATOL 372C and the corresponding computed solution 373C component is zero. Thus, a pure relative error 374C test is impossible for this component. 375C 376C IDID = -4 -- The problem appears to be stiff. 377C 378C IDID = -5,-6,-7,..,-32 -- Not applicable for this code 379C but used by other members of DEPAC or possible 380C future extensions. 381C 382C *** Task Terminated *** 383C Reported by the value of IDID=-33 384C 385C IDID = -33 -- The code has encountered trouble from which 386C it cannot recover. A message is printed 387C explaining the trouble and control is returned 388C to the calling program. For example, this occurs 389C when invalid input is detected. 390C 391C RTOL, ATOL -- These quantities remain unchanged except when 392C IDID = -2. In this case, the error tolerances have been 393C increased by the code to values which are estimated to be 394C appropriate for continuing the integration. However, the 395C reported solution at T was obtained using the input values 396C of RTOL and ATOL. 397C 398C RWORK, IWORK -- Contain information which is usually of no 399C interest to the user but necessary for subsequent calls. 400C However, you may find use for 401C 402C RWORK(11)--which contains the step size H to be 403C attempted on the next step. 404C 405C RWORK(12)--if the tolerances have been increased by the 406C code (IDID = -2) , they were multiplied by the 407C value in RWORK(12). 408C 409C RWORK(13)--Which contains the current value of the 410C independent variable, i.e. the farthest point 411C integration has reached. This will be different 412C from T only when interpolation has been 413C performed (IDID=3). 414C 415C RWORK(20+I)--Which contains the approximate derivative 416C of the solution component Y(I). In DDEABM, it 417C is obtained by calling subroutine DF to 418C evaluate the differential equation using T and 419C Y(*) when IDID=1 or 2, and by interpolation 420C when IDID=3. 421C 422C ********************************************************************** 423C * INPUT -- What To Do To Continue The Integration * 424C * (calls after the first) * 425C ********************************************************************** 426C 427C This code is organized so that subsequent calls to continue the 428C integration involve little (if any) additional effort on your 429C part. You must monitor the IDID parameter in order to determine 430C what to do next. 431C 432C Recalling that the principal task of the code is to integrate 433C from T to TOUT (the interval mode), usually all you will need 434C to do is specify a new TOUT upon reaching the current TOUT. 435C 436C Do not alter any quantity not specifically permitted below, 437C in particular do not alter NEQ, T, Y(*), RWORK(*), IWORK(*) or 438C the differential equation in subroutine DF. Any such alteration 439C constitutes a new problem and must be treated as such, i.e. 440C you must start afresh. 441C 442C You cannot change from vector to scalar error control or vice 443C versa (INFO(2)) but you can change the size of the entries of 444C RTOL, ATOL. Increasing a tolerance makes the equation easier 445C to integrate. Decreasing a tolerance will make the equation 446C harder to integrate and should generally be avoided. 447C 448C You can switch from the intermediate-output mode to the 449C interval mode (INFO(3)) or vice versa at any time. 450C 451C If it has been necessary to prevent the integration from going 452C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the 453C code will not integrate to any TOUT beyond the currently 454C specified TSTOP. Once TSTOP has been reached you must change 455C the value of TSTOP or set INFO(4)=0. You may change INFO(4) 456C or TSTOP at any time but you must supply the value of TSTOP in 457C RWORK(1) whenever you set INFO(4)=1. 458C 459C The parameter INFO(1) is used by the code to indicate the 460C beginning of a new problem and to indicate whether integration 461C is to be continued. You must input the value INFO(1) = 0 462C when starting a new problem. You must input the value 463C INFO(1) = 1 if you wish to continue after an interrupted task. 464C Do not set INFO(1) = 0 on a continuation call unless you 465C want the code to restart at the current T. 466C 467C *** Following A Completed Task *** 468C If 469C IDID = 1, call the code again to continue the integration 470C another step in the direction of TOUT. 471C 472C IDID = 2 or 3, define a new TOUT and call the code again. 473C TOUT must be different from T. You cannot change 474C the direction of integration without restarting. 475C 476C *** Following An Interrupted Task *** 477C To show the code that you realize the task was 478C interrupted and that you want to continue, you 479C must take appropriate action and reset INFO(1) = 1 480C If 481C IDID = -1, the code has attempted 500 steps. 482C If you want to continue, set INFO(1) = 1 and 483C call the code again. An additional 500 steps 484C will be allowed. 485C 486C IDID = -2, the error tolerances RTOL, ATOL have been 487C increased to values the code estimates appropriate 488C for continuing. You may want to change them 489C yourself. If you are sure you want to continue 490C with relaxed error tolerances, set INFO(1)=1 and 491C call the code again. 492C 493C IDID = -3, a solution component is zero and you set the 494C corresponding component of ATOL to zero. If you 495C are sure you want to continue, you must first 496C alter the error criterion to use positive values 497C for those components of ATOL corresponding to zero 498C solution components, then set INFO(1)=1 and call 499C the code again. 500C 501C IDID = -4, the problem appears to be stiff. It is very 502C inefficient to solve such problems with DDEABM. 503C The code DDEBDF in DEPAC handles this task 504C efficiently. If you are absolutely sure you want 505C to continue with DDEABM, set INFO(1)=1 and call 506C the code again. 507C 508C IDID = -5,-6,-7,..,-32 --- cannot occur with this code 509C but used by other members of DEPAC or possible 510C future extensions. 511C 512C *** Following A Terminated Task *** 513C If 514C IDID = -33, you cannot continue the solution of this 515C problem. An attempt to do so will result in your 516C run being terminated. 517C 518C ********************************************************************** 519C *Long Description: 520C 521C ********************************************************************** 522C * DEPAC Package Overview * 523C ********************************************************************** 524C 525C .... You have a choice of three differential equation solvers from 526C .... DEPAC. The following brief descriptions are meant to aid you in 527C .... choosing the most appropriate code for your problem. 528C 529C .... DDERKF is a fifth order Runge-Kutta code. It is the simplest of 530C .... the three choices, both algorithmically and in the use of the 531C .... code. DDERKF is primarily designed to solve non-stiff and 532C .... mildly stiff differential equations when derivative evaluations 533C .... are not expensive. It should generally not be used to get high 534C .... accuracy results nor answers at a great many specific points. 535C .... Because DDERKF has very low overhead costs, it will usually 536C .... result in the least expensive integration when solving 537C .... problems requiring a modest amount of accuracy and having 538C .... equations that are not costly to evaluate. DDERKF attempts to 539C .... discover when it is not suitable for the task posed. 540C 541C .... DDEABM is a variable order (one through twelve) Adams code. 542C .... Its complexity lies somewhere between that of DDERKF and 543C .... DDEBDF. DDEABM is primarily designed to solve non-stiff and 544C .... mildly stiff differential equations when derivative evaluations 545C .... are expensive, high accuracy results are needed or answers at 546C .... many specific points are required. DDEABM attempts to discover 547C .... when it is not suitable for the task posed. 548C 549C .... DDEBDF is a variable order (one through five) backward 550C .... differentiation formula code. it is the most complicated of 551C .... the three choices. DDEBDF is primarily designed to solve stiff 552C .... differential equations at crude to moderate tolerances. 553C .... If the problem is very stiff at all, DDERKF and DDEABM will be 554C .... quite inefficient compared to DDEBDF. However, DDEBDF will be 555C .... inefficient compared to DDERKF and DDEABM on non-stiff problems 556C .... because it uses much more storage, has a much larger overhead, 557C .... and the low order formulas will not give high accuracies 558C .... efficiently. 559C 560C .... The concept of stiffness cannot be described in a few words. 561C .... If you do not know the problem to be stiff, try either DDERKF 562C .... or DDEABM. Both of these codes will inform you of stiffness 563C .... when the cost of solving such problems becomes important. 564C 565C ********************************************************************* 566C 567C***REFERENCES L. F. Shampine and H. A. Watts, DEPAC - design of a user 568C oriented package of ODE solvers, Report SAND79-2374, 569C Sandia Laboratories, 1979. 570C***ROUTINES CALLED DDES, XERMSG 571C***REVISION HISTORY (YYMMDD) 572C 820301 DATE WRITTEN 573C 890531 Changed all specific intrinsics to generic. (WRB) 574C 890831 Modified array declarations. (WRB) 575C 891006 Cosmetic changes to prologue. (WRB) 576C 891024 Changed references from DVNORM to DHVNRM. (WRB) 577C 891024 REVISION DATE from Version 3.2 578C 891214 Prologue converted to Version 4.0 format. (BAB) 579C 900510 Convert XERRWV calls to XERMSG calls. (RWC) 580C 920501 Reformatted the REFERENCES section. (WRB) 581C***END PROLOGUE DDEABM 582C 583 INTEGER IALPHA, IBETA, IDELSN, IDID, IFOURU, IG, IHOLD, 584 1 INFO, IP, IPAR, IPHI, IPSI, ISIG, ITOLD, ITSTAR, ITWOU, 585 2 IV, IW, IWORK, IWT, IYP, IYPOUT, IYY, LIW, LRW, NEQ 586 DOUBLE PRECISION ATOL, RPAR, RTOL, RWORK, T, TOUT, Y 587 LOGICAL START,PHASE1,NORND,STIFF,INTOUT 588C 589 DIMENSION Y(*),INFO(15),RTOL(*),ATOL(*),RWORK(*),IWORK(*), 590 1 RPAR(*),IPAR(*) 591C 592 CHARACTER*8 XERN1 593 CHARACTER*16 XERN3 594C 595 EXTERNAL DF 596C 597C CHECK FOR AN APPARENT INFINITE LOOP 598C 599C***FIRST EXECUTABLE STATEMENT DDEABM 600 IF ( INFO(1) .EQ. 0 ) IWORK(LIW) = 0 601 IF (IWORK(LIW) .GE. 5) THEN 602 IF (T .EQ. RWORK(21 + NEQ)) THEN 603 WRITE (XERN3, '(1PE15.6)') T 604 CALL XERMSG ('SLATEC', 'DDEABM', 605 * 'AN APPARENT INFINITE LOOP HAS BEEN DETECTED.$$' // 606 * 'YOU HAVE MADE REPEATED CALLS AT T = ' // XERN3 // 607 * ' AND THE INTEGRATION HAS NOT ADVANCED. CHECK THE ' // 608 * 'WAY YOU HAVE SET PARAMETERS FOR THE CALL TO THE ' // 609 * 'CODE, PARTICULARLY INFO(1).', 13, 2) 610 RETURN 611 ENDIF 612 ENDIF 613C 614C CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION 615C 616 IDID=0 617 IF (LRW .LT. 130+21*NEQ) THEN 618 WRITE (XERN1, '(I8)') LRW 619 CALL XERMSG ('SLATEC', 'DDEABM', 'THE LENGTH OF THE RWORK ' // 620 * 'ARRAY MUST BE AT LEAST 130 + 21*NEQ.$$' // 621 * 'YOU HAVE CALLED THE CODE WITH LRW = ' // XERN1, 1, 1) 622 IDID=-33 623 ENDIF 624C 625 IF (LIW .LT. 51) THEN 626 WRITE (XERN1, '(I8)') LIW 627 CALL XERMSG ('SLATEC', 'DDEABM', 'THE LENGTH OF THE IWORK ' // 628 * 'ARRAY MUST BE AT LEAST 51.$$YOU HAVE CALLED THE CODE ' // 629 * 'WITH LIW = ' // XERN1, 2, 1) 630 IDID=-33 631 ENDIF 632C 633C COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK ARRAY 634C 635 IYPOUT = 21 636 ITSTAR = NEQ + 21 637 IYP = 1 + ITSTAR 638 IYY = NEQ + IYP 639 IWT = NEQ + IYY 640 IP = NEQ + IWT 641 IPHI = NEQ + IP 642 IALPHA = (NEQ*16) + IPHI 643 IBETA = 12 + IALPHA 644 IPSI = 12 + IBETA 645 IV = 12 + IPSI 646 IW = 12 + IV 647 ISIG = 12 + IW 648 IG = 13 + ISIG 649 IGI = 13 + IG 650 IXOLD = 11 + IGI 651 IHOLD = 1 + IXOLD 652 ITOLD = 1 + IHOLD 653 IDELSN = 1 + ITOLD 654 ITWOU = 1 + IDELSN 655 IFOURU = 1 + ITWOU 656C 657 RWORK(ITSTAR) = T 658 IF (INFO(1) .EQ. 0) GO TO 50 659 START = IWORK(21) .NE. (-1) 660 PHASE1 = IWORK(22) .NE. (-1) 661 NORND = IWORK(23) .NE. (-1) 662 STIFF = IWORK(24) .NE. (-1) 663 INTOUT = IWORK(25) .NE. (-1) 664C 665 50 CALL DDES(DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,RWORK(IYPOUT), 666 1 RWORK(IYP),RWORK(IYY),RWORK(IWT),RWORK(IP),RWORK(IPHI), 667 2 RWORK(IALPHA),RWORK(IBETA),RWORK(IPSI),RWORK(IV), 668 3 RWORK(IW),RWORK(ISIG),RWORK(IG),RWORK(IGI),RWORK(11), 669 4 RWORK(12),RWORK(13),RWORK(IXOLD),RWORK(IHOLD), 670 5 RWORK(ITOLD),RWORK(IDELSN),RWORK(1),RWORK(ITWOU), 671 5 RWORK(IFOURU),START,PHASE1,NORND,STIFF,INTOUT,IWORK(26), 672 6 IWORK(27),IWORK(28),IWORK(29),IWORK(30),IWORK(31), 673 7 IWORK(32),IWORK(33),IWORK(34),IWORK(35),IWORK(45), 674 8 RPAR,IPAR) 675C 676 IWORK(21) = -1 677 IF (START) IWORK(21) = 1 678 IWORK(22) = -1 679 IF (PHASE1) IWORK(22) = 1 680 IWORK(23) = -1 681 IF (NORND) IWORK(23) = 1 682 IWORK(24) = -1 683 IF (STIFF) IWORK(24) = 1 684 IWORK(25) = -1 685 IF (INTOUT) IWORK(25) = 1 686C 687 IF (IDID .NE. (-2)) IWORK(LIW) = IWORK(LIW) + 1 688 IF (T .NE. RWORK(ITSTAR)) IWORK(LIW) = 0 689C 690 RETURN 691 END 692*DECK DDES 693 SUBROUTINE DDES (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID, 694 + YPOUT, YP, YY, WT, P, PHI, ALPHA, BETA, PSI, V, W, SIG, G, GI, 695 + H, EPS, X, XOLD, HOLD, TOLD, DELSGN, TSTOP, TWOU, FOURU, START, 696 + PHASE1, NORND, STIFF, INTOUT, NS, KORD, KOLD, INIT, KSTEPS, 697 + KLE4, IQUIT, KPREV, IVC, IV, KGI, RPAR, IPAR) 698C***BEGIN PROLOGUE DDES 699C***SUBSIDIARY 700C***PURPOSE Subsidiary to DDEABM 701C***LIBRARY SLATEC 702C***TYPE DOUBLE PRECISION (DES-S, DDES-D) 703C***AUTHOR Watts, H. A., (SNLA) 704C***DESCRIPTION 705C 706C DDEABM merely allocates storage for DDES to relieve the user of the 707C inconvenience of a long call list. Consequently DDES is used as 708C described in the comments for DDEABM . 709C 710C***SEE ALSO DDEABM 711C***ROUTINES CALLED D1MACH, DINTP, DSTEPS, XERMSG 712C***REVISION HISTORY (YYMMDD) 713C 820301 DATE WRITTEN 714C 890531 Changed all specific intrinsics to generic. (WRB) 715C 890831 Modified array declarations. (WRB) 716C 891214 Prologue converted to Version 4.0 format. (BAB) 717C 900328 Added TYPE section. (WRB) 718C 900510 Convert XERRWV calls to XERMSG calls, cvt GOTOs to 719C IF-THEN-ELSE. (RWC) 720C 910722 Updated AUTHOR section. (ALS) 721C***END PROLOGUE DDES 722C 723 INTEGER IDID, INFO, INIT, IPAR, IQUIT, IV, IVC, K, KGI, KLE4, 724 1 KOLD, KORD, KPREV, KSTEPS, L, LTOL, MAXNUM, NATOLP, NEQ, 725 2 NRTOLP, NS 726 DOUBLE PRECISION A, ABSDEL, ALPHA, ATOL, BETA, D1MACH, 727 1 DEL, DELSGN, DT, EPS, FOURU, G, GI, H, 728 2 HA, HOLD, P, PHI, PSI, RPAR, RTOL, SIG, T, TOLD, TOUT, 729 3 TSTOP, TWOU, U, V, W, WT, X, XOLD, Y, YP, YPOUT, YY 730 LOGICAL STIFF,CRASH,START,PHASE1,NORND,INTOUT 731C 732 DIMENSION Y(*),YY(*),WT(*),PHI(NEQ,16),P(*),YP(*), 733 1 YPOUT(*),PSI(12),ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13), 734 2 GI(11),IV(10),INFO(15),RTOL(*),ATOL(*),RPAR(*),IPAR(*) 735 CHARACTER*8 XERN1 736 CHARACTER*16 XERN3, XERN4 737C 738 EXTERNAL DF 739C 740C....................................................................... 741C 742C THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE 743C NUMBER OF STEPS ATTEMPTED. WHEN THIS EXCEEDS MAXNUM, THE COUNTER 744C IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE EXCESSIVE 745C WORK. 746C 747 SAVE MAXNUM 748 DATA MAXNUM/500/ 749C 750C....................................................................... 751C 752C***FIRST EXECUTABLE STATEMENT DDES 753 IF (INFO(1) .EQ. 0) THEN 754C 755C ON THE FIRST CALL , PERFORM INITIALIZATION -- 756C DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY U BY CALLING THE 757C FUNCTION ROUTINE D1MACH. THE USER MUST MAKE SURE THAT THE 758C VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED. 759C 760 U=D1MACH(4) 761C -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS 762 TWOU=2.D0*U 763 FOURU=4.D0*U 764C -- SET TERMINATION FLAG 765 IQUIT=0 766C -- SET INITIALIZATION INDICATOR 767 INIT=0 768C -- SET COUNTER FOR ATTEMPTED STEPS 769 KSTEPS=0 770C -- SET INDICATOR FOR INTERMEDIATE-OUTPUT 771 INTOUT= .FALSE. 772C -- SET INDICATOR FOR STIFFNESS DETECTION 773 STIFF= .FALSE. 774C -- SET STEP COUNTER FOR STIFFNESS DETECTION 775 KLE4=0 776C -- SET INDICATORS FOR STEPS CODE 777 START= .TRUE. 778 PHASE1= .TRUE. 779 NORND= .TRUE. 780C -- RESET INFO(1) FOR SUBSEQUENT CALLS 781 INFO(1)=1 782 ENDIF 783C 784C....................................................................... 785C 786C CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY 787C 788 IF (INFO(1) .NE. 0 .AND. INFO(1) .NE. 1) THEN 789 WRITE (XERN1, '(I8)') INFO(1) 790 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INFO(1) MUST BE ' // 791 * 'SET TO 0 FOR THE START OF A NEW PROBLEM, AND MUST BE ' // 792 * 'SET TO 1 FOLLOWING AN INTERRUPTED TASK. YOU ARE ' // 793 * 'ATTEMPTING TO CONTINUE THE INTEGRATION ILLEGALLY BY ' // 794 * 'CALLING THE 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', 'DDES', 'IN DDEABM, INFO(2) MUST BE ' // 801 * '0 OR 1 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', 'DDES', 'IN DDEABM, INFO(3) MUST BE ' // 810 * '0 OR 1 INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT ' // 811 * 'MODE OF INTEGRATION, RESPECTIVELY. YOU HAVE CALLED ' // 812 * 'THE CODE 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', 'DDES', 'IN DDEABM, INFO(4) MUST BE ' // 819 * '0 OR 1 INDICATING WHETHER OR NOT THE INTEGRATION ' // 820 * 'INTERVAL IS TO BE RESTRICTED BY A POINT TSTOP. YOU ' // 821 * 'HAVE CALLED THE CODE WITH INFO(4) = ' // XERN1, 14, 1) 822 IDID=-33 823 ENDIF 824C 825 IF (NEQ .LT. 1) THEN 826 WRITE (XERN1, '(I8)') NEQ 827 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, THE NUMBER OF ' // 828 * 'EQUATIONS NEQ MUST BE A POSITIVE INTEGER. YOU HAVE ' // 829 * 'CALLED THE CODE WITH NEQ = ' // XERN1, 6, 1) 830 IDID=-33 831 ENDIF 832C 833 NRTOLP = 0 834 NATOLP = 0 835 DO 90 K=1,NEQ 836 IF (NRTOLP .EQ. 0 .AND. RTOL(K) .LT. 0.D0) THEN 837 WRITE (XERN1, '(I8)') K 838 WRITE (XERN3, '(1PE15.6)') RTOL(K) 839 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, THE RELATIVE ' // 840 * 'ERROR TOLERANCES RTOL MUST BE NON-NEGATIVE. YOU ' // 841 * 'HAVE CALLED THE CODE WITH RTOL(' // XERN1 // ') = ' // 842 * XERN3 // '. IN THE CASE OF VECTOR ERROR TOLERANCES, ' // 843 * 'NO FURTHER CHECKING OF RTOL COMPONENTS IS DONE.', 7, 1) 844 IDID = -33 845 NRTOLP = 1 846 ENDIF 847C 848 IF (NATOLP .EQ. 0 .AND. ATOL(K) .LT. 0.D0) THEN 849 WRITE (XERN1, '(I8)') K 850 WRITE (XERN3, '(1PE15.6)') ATOL(K) 851 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, THE ABSOLUTE ' // 852 * 'ERROR TOLERANCES ATOL MUST BE NON-NEGATIVE. YOU ' // 853 * 'HAVE CALLED THE CODE WITH ATOL(' // XERN1 // ') = ' // 854 * XERN3 // '. IN THE CASE OF VECTOR ERROR TOLERANCES, ' // 855 * 'NO FURTHER CHECKING OF ATOL COMPONENTS IS DONE.', 8, 1) 856 IDID = -33 857 NATOLP = 1 858 ENDIF 859C 860 IF (INFO(2) .EQ. 0) GO TO 100 861 IF (NATOLP.GT.0 .AND. NRTOLP.GT.0) GO TO 100 862 90 CONTINUE 863C 864 100 IF (INFO(4) .EQ. 1) THEN 865 IF (SIGN(1.D0,TOUT-T) .NE. SIGN(1.D0,TSTOP-T) 866 1 .OR. ABS(TOUT-T) .GT. ABS(TSTOP-T)) THEN 867 WRITE (XERN3, '(1PE15.6)') TOUT 868 WRITE (XERN4, '(1PE15.6)') TSTOP 869 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, YOU HAVE ' // 870 * 'CALLED THE CODE WITH TOUT = ' // XERN3 // ' BUT ' // 871 * 'YOU HAVE ALSO TOLD THE CODE (INFO(4) = 1) NOT TO ' // 872 * 'INTEGRATE PAST THE POINT TSTOP = ' // XERN4 // 873 * ' THESE INSTRUCTIONS CONFLICT.', 14, 1) 874 IDID=-33 875 ENDIF 876 ENDIF 877C 878C CHECK SOME CONTINUATION POSSIBILITIES 879C 880 IF (INIT .NE. 0) THEN 881 IF (T .EQ. TOUT) THEN 882 WRITE (XERN3, '(1PE15.6)') T 883 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, YOU HAVE ' // 884 * 'CALLED THE CODE WITH T = TOUT = ' // XERN3 // 885 * '$$THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 9, 1) 886 IDID=-33 887 ENDIF 888C 889 IF (T .NE. TOLD) THEN 890 WRITE (XERN3, '(1PE15.6)') TOLD 891 WRITE (XERN4, '(1PE15.6)') T 892 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, YOU HAVE ' // 893 * 'CHANGED THE VALUE OF T FROM ' // XERN3 // ' TO ' // 894 * XERN4 //' THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 895 * 10, 1) 896 IDID=-33 897 ENDIF 898C 899 IF (INIT .NE. 1) THEN 900 IF (DELSGN*(TOUT-T) .LT. 0.D0) THEN 901 WRITE (XERN3, '(1PE15.6)') TOUT 902 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, BY ' // 903 * 'CALLING THE CODE WITH TOUT = ' // XERN3 // 904 * ' YOU ARE ATTEMPTING TO CHANGE THE DIRECTION OF ' // 905 * 'INTEGRATION.$$THIS IS NOT ALLOWED WITHOUT ' // 906 * 'RESTARTING.', 11, 1) 907 IDID=-33 908 ENDIF 909 ENDIF 910 ENDIF 911C 912C INVALID INPUT DETECTED 913C 914 IF (IDID .EQ. (-33)) THEN 915 IF (IQUIT .NE. (-33)) THEN 916 IQUIT = -33 917 INFO(1) = -1 918 ELSE 919 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INVALID ' // 920 * 'INPUT WAS DETECTED ON SUCCESSIVE ENTRIES. IT IS ' // 921 * 'IMPOSSIBLE TO PROCEED BECAUSE YOU HAVE NOT ' // 922 * 'CORRECTED THE PROBLEM, SO EXECUTION IS BEING ' // 923 * 'TERMINATED.', 12, 2) 924 ENDIF 925 RETURN 926 ENDIF 927C 928C....................................................................... 929C 930C RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED AS 931C ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS CASE, 932C THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE SMALLEST VALUE 933C FOURU WHICH IS LIKELY TO BE REASONABLE FOR THIS METHOD AND MACHINE 934C 935 DO 180 K=1,NEQ 936 IF (RTOL(K)+ATOL(K) .GT. 0.D0) GO TO 170 937 RTOL(K)=FOURU 938 IDID=-2 939 170 IF (INFO(2) .EQ. 0) GO TO 190 940 180 CONTINUE 941C 942 190 IF (IDID .NE. (-2)) GO TO 200 943C RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A 944C SMALL POSITIVE VALUE 945 INFO(1)=-1 946 RETURN 947C 948C BRANCH ON STATUS OF INITIALIZATION INDICATOR 949C INIT=0 MEANS INITIAL DERIVATIVES AND NOMINAL STEP SIZE 950C AND DIRECTION NOT YET SET 951C INIT=1 MEANS NOMINAL STEP SIZE AND DIRECTION NOT YET SET 952C INIT=2 MEANS NO FURTHER INITIALIZATION REQUIRED 953C 954 200 IF (INIT .EQ. 0) GO TO 210 955 IF (INIT .EQ. 1) GO TO 220 956 GO TO 240 957C 958C....................................................................... 959C 960C MORE INITIALIZATION -- 961C -- EVALUATE INITIAL DERIVATIVES 962C 963 210 INIT=1 964 A=T 965 CALL DF(A,Y,YP,RPAR,IPAR) 966 IF (T .NE. TOUT) GO TO 220 967 IDID=2 968 DO 215 L = 1,NEQ 969 215 YPOUT(L) = YP(L) 970 TOLD=T 971 RETURN 972C 973C -- SET INDEPENDENT AND DEPENDENT VARIABLES 974C X AND YY(*) FOR STEPS 975C -- SET SIGN OF INTEGRATION DIRECTION 976C -- INITIALIZE THE STEP SIZE 977C 978 220 INIT = 2 979 X = T 980 DO 230 L = 1,NEQ 981 230 YY(L) = Y(L) 982 DELSGN = SIGN(1.0D0,TOUT-T) 983 H = SIGN(MAX(FOURU*ABS(X),ABS(TOUT-X)),TOUT-X) 984C 985C....................................................................... 986C 987C ON EACH CALL SET INFORMATION WHICH DETERMINES THE ALLOWED INTERVAL 988C OF INTEGRATION BEFORE RETURNING WITH AN ANSWER AT TOUT 989C 990 240 DEL = TOUT - T 991 ABSDEL = ABS(DEL) 992C 993C....................................................................... 994C 995C IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN 996C 997 250 IF(ABS(X-T) .LT. ABSDEL) GO TO 260 998 CALL DINTP(X,YY,TOUT,Y,YPOUT,NEQ,KOLD,PHI,IVC,IV,KGI,GI, 999 1 ALPHA,G,W,XOLD,P) 1000 IDID = 3 1001 IF (X .NE. TOUT) GO TO 255 1002 IDID = 2 1003 INTOUT = .FALSE. 1004 255 T = TOUT 1005 TOLD = T 1006 RETURN 1007C 1008C IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE, 1009C EXTRAPOLATE AND RETURN 1010C 1011 260 IF (INFO(4) .NE. 1) GO TO 280 1012 IF (ABS(TSTOP-X) .GE. FOURU*ABS(X)) GO TO 280 1013 DT = TOUT - X 1014 DO 270 L = 1,NEQ 1015 270 Y(L) = YY(L) + DT*YP(L) 1016 CALL DF(TOUT,Y,YPOUT,RPAR,IPAR) 1017 IDID = 3 1018 T = TOUT 1019 TOLD = T 1020 RETURN 1021C 1022 280 IF (INFO(3) .EQ. 0 .OR. .NOT.INTOUT) GO TO 300 1023C 1024C INTERMEDIATE-OUTPUT MODE 1025C 1026 IDID = 1 1027 DO 290 L = 1,NEQ 1028 Y(L)=YY(L) 1029 290 YPOUT(L) = YP(L) 1030 T = X 1031 TOLD = T 1032 INTOUT = .FALSE. 1033 RETURN 1034C 1035C....................................................................... 1036C 1037C MONITOR NUMBER OF STEPS ATTEMPTED 1038C 1039 300 IF (KSTEPS .LE. MAXNUM) GO TO 330 1040C 1041C A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED 1042 IDID=-1 1043 KSTEPS=0 1044 IF (.NOT. STIFF) GO TO 310 1045C 1046C PROBLEM APPEARS TO BE STIFF 1047 IDID=-4 1048 STIFF= .FALSE. 1049 KLE4=0 1050C 1051 310 DO 320 L = 1,NEQ 1052 Y(L) = YY(L) 1053 320 YPOUT(L) = YP(L) 1054 T = X 1055 TOLD = T 1056 INFO(1) = -1 1057 INTOUT = .FALSE. 1058 RETURN 1059C 1060C....................................................................... 1061C 1062C LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP 1063C 1064 330 HA = ABS(H) 1065 IF (INFO(4) .NE. 1) GO TO 340 1066 HA = MIN(HA,ABS(TSTOP-X)) 1067 340 H = SIGN(HA,H) 1068 EPS = 1.0D0 1069 LTOL = 1 1070 DO 350 L = 1,NEQ 1071 IF (INFO(2) .EQ. 1) LTOL = L 1072 WT(L) = RTOL(LTOL)*ABS(YY(L)) + ATOL(LTOL) 1073 IF (WT(L) .LE. 0.0D0) GO TO 360 1074 350 CONTINUE 1075 GO TO 380 1076C 1077C RELATIVE ERROR CRITERION INAPPROPRIATE 1078 360 IDID = -3 1079 DO 370 L = 1,NEQ 1080 Y(L) = YY(L) 1081 370 YPOUT(L) = YP(L) 1082 T = X 1083 TOLD = T 1084 INFO(1) = -1 1085 INTOUT = .FALSE. 1086 RETURN 1087C 1088 380 CALL DSTEPS(DF,NEQ,YY,X,H,EPS,WT,START,HOLD,KORD,KOLD,CRASH,PHI,P, 1089 1 YP,PSI,ALPHA,BETA,SIG,V,W,G,PHASE1,NS,NORND,KSTEPS, 1090 2 TWOU,FOURU,XOLD,KPREV,IVC,IV,KGI,GI,RPAR,IPAR) 1091C 1092C....................................................................... 1093C 1094 IF(.NOT.CRASH) GO TO 420 1095C 1096C TOLERANCES TOO SMALL 1097 IDID = -2 1098 RTOL(1) = EPS*RTOL(1) 1099 ATOL(1) = EPS*ATOL(1) 1100 IF (INFO(2) .EQ. 0) GO TO 400 1101 DO 390 L = 2,NEQ 1102 RTOL(L) = EPS*RTOL(L) 1103 390 ATOL(L) = EPS*ATOL(L) 1104 400 DO 410 L = 1,NEQ 1105 Y(L) = YY(L) 1106 410 YPOUT(L) = YP(L) 1107 T = X 1108 TOLD = T 1109 INFO(1) = -1 1110 INTOUT = .FALSE. 1111 RETURN 1112C 1113C (STIFFNESS TEST) COUNT NUMBER OF CONSECUTIVE STEPS TAKEN WITH THE 1114C ORDER OF THE METHOD BEING LESS OR EQUAL TO FOUR 1115C 1116 420 KLE4 = KLE4 + 1 1117 IF(KOLD .GT. 4) KLE4 = 0 1118 IF(KLE4 .GE. 50) STIFF = .TRUE. 1119 INTOUT = .TRUE. 1120 GO TO 250 1121 END 1122*DECK DINTP 1123 SUBROUTINE DINTP (X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, IVC, 1124 + IV, KGI, GI, ALPHA, OG, OW, OX, OY) 1125C***BEGIN PROLOGUE DINTP 1126C***PURPOSE Approximate the solution at XOUT by evaluating the 1127C polynomial computed in DSTEPS at XOUT. Must be used in 1128C conjunction with DSTEPS. 1129C***LIBRARY SLATEC (DEPAC) 1130C***CATEGORY I1A1B 1131C***TYPE DOUBLE PRECISION (SINTRP-S, DINTP-D) 1132C***KEYWORDS ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE, 1133C ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR, 1134C SMOOTH INTERPOLANT 1135C***AUTHOR Watts, H. A., (SNLA) 1136C***DESCRIPTION 1137C 1138C The methods in subroutine DSTEPS approximate the solution near X 1139C by a polynomial. Subroutine DINTP approximates the solution at 1140C XOUT by evaluating the polynomial there. Information defining this 1141C polynomial is passed from DSTEPS so DINTP cannot be used alone. 1142C 1143C Subroutine DSTEPS is completely explained and documented in the text 1144C "Computer Solution of Ordinary Differential Equations, the Initial 1145C Value Problem" by L. F. Shampine and M. K. Gordon. 1146C 1147C Input to DINTP -- 1148C 1149C The user provides storage in the calling program for the arrays in 1150C the call list 1151C DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN) 1152C AND ALPHA(12),OG(13),OW(12),GI(11),IV(10) 1153C and defines 1154C XOUT -- point at which solution is desired. 1155C The remaining parameters are defined in DSTEPS and passed to 1156C DINTP from that subroutine 1157C 1158C Output from DINTP -- 1159C 1160C YOUT(*) -- solution at XOUT 1161C YPOUT(*) -- derivative of solution at XOUT 1162C The remaining parameters are returned unaltered from their input 1163C values. Integration with DSTEPS may be continued. 1164C 1165C***REFERENCES H. A. Watts, A smoother interpolant for DE/STEP, INTRP 1166C II, Report SAND84-0293, Sandia Laboratories, 1984. 1167C***ROUTINES CALLED (NONE) 1168C***REVISION HISTORY (YYMMDD) 1169C 840201 DATE WRITTEN 1170C 890831 Modified array declarations. (WRB) 1171C 890831 REVISION DATE from Version 3.2 1172C 891214 Prologue converted to Version 4.0 format. (BAB) 1173C 920501 Reformatted the REFERENCES section. (WRB) 1174C***END PROLOGUE DINTP 1175C 1176 INTEGER I, IQ, IV, IVC, IW, J, JQ, KGI, KOLD, KP1, KP2, 1177 1 L, M, NEQN 1178 DOUBLE PRECISION ALP, ALPHA, C, G, GDI, GDIF, GI, GAMMA, H, HI, 1179 1 HMU, OG, OW, OX, OY, PHI, RMU, SIGMA, TEMP1, TEMP2, TEMP3, 1180 2 W, X, XI, XIM1, XIQ, XOUT, Y, YOUT, YPOUT 1181C 1182 DIMENSION Y(*),YOUT(*),YPOUT(*),PHI(NEQN,16),OY(*) 1183 DIMENSION G(13),C(13),W(13),OG(13),OW(12),ALPHA(12),GI(11),IV(10) 1184C 1185C***FIRST EXECUTABLE STATEMENT DINTP 1186 KP1 = KOLD + 1 1187 KP2 = KOLD + 2 1188C 1189 HI = XOUT - OX 1190 H = X - OX 1191 XI = HI/H 1192 XIM1 = XI - 1.D0 1193C 1194C INITIALIZE W(*) FOR COMPUTING G(*) 1195C 1196 XIQ = XI 1197 DO 10 IQ = 1,KP1 1198 XIQ = XI*XIQ 1199 TEMP1 = IQ*(IQ+1) 1200 10 W(IQ) = XIQ/TEMP1 1201C 1202C COMPUTE THE DOUBLE INTEGRAL TERM GDI 1203C 1204 IF (KOLD .LE. KGI) GO TO 50 1205 IF (IVC .GT. 0) GO TO 20 1206 GDI = 1.0D0/TEMP1 1207 M = 2 1208 GO TO 30 1209 20 IW = IV(IVC) 1210 GDI = OW(IW) 1211 M = KOLD - IW + 3 1212 30 IF (M .GT. KOLD) GO TO 60 1213 DO 40 I = M,KOLD 1214 40 GDI = OW(KP2-I) - ALPHA(I)*GDI 1215 GO TO 60 1216 50 GDI = GI(KOLD) 1217C 1218C COMPUTE G(*) AND C(*) 1219C 1220 60 G(1) = XI 1221 G(2) = 0.5D0*XI*XI 1222 C(1) = 1.0D0 1223 C(2) = XI 1224 IF (KOLD .LT. 2) GO TO 90 1225 DO 80 I = 2,KOLD 1226 ALP = ALPHA(I) 1227 GAMMA = 1.0D0 + XIM1*ALP 1228 L = KP2 - I 1229 DO 70 JQ = 1,L 1230 70 W(JQ) = GAMMA*W(JQ) - ALP*W(JQ+1) 1231 G(I+1) = W(1) 1232 80 C(I+1) = GAMMA*C(I) 1233C 1234C DEFINE INTERPOLATION PARAMETERS 1235C 1236 90 SIGMA = (W(2) - XIM1*W(1))/GDI 1237 RMU = XIM1*C(KP1)/GDI 1238 HMU = RMU/H 1239C 1240C INTERPOLATE FOR THE SOLUTION -- YOUT 1241C AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT 1242C 1243 DO 100 L = 1,NEQN 1244 YOUT(L) = 0.0D0 1245 100 YPOUT(L) = 0.0D0 1246 DO 120 J = 1,KOLD 1247 I = KP2 - J 1248 GDIF = OG(I) - OG(I-1) 1249 TEMP2 = (G(I) - G(I-1)) - SIGMA*GDIF 1250 TEMP3 = (C(I) - C(I-1)) + RMU*GDIF 1251 DO 110 L = 1,NEQN 1252 YOUT(L) = YOUT(L) + TEMP2*PHI(L,I) 1253 110 YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I) 1254 120 CONTINUE 1255 DO 130 L = 1,NEQN 1256 YOUT(L) = ((1.0D0 - SIGMA)*OY(L) + SIGMA*Y(L)) + 1257 1 H*(YOUT(L) + (G(1) - SIGMA*OG(1))*PHI(L,1)) 1258 130 YPOUT(L) = HMU*(OY(L) - Y(L)) + 1259 1 (YPOUT(L) + (C(1) + RMU*OG(1))*PHI(L,1)) 1260C 1261 RETURN 1262 END 1263*DECK XERMSG 1264 SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL) 1265C***BEGIN PROLOGUE XERMSG 1266C***PURPOSE Process error messages for SLATEC and other libraries. 1267C***LIBRARY SLATEC (XERROR) 1268C***CATEGORY R3C 1269C***TYPE ALL (XERMSG-A) 1270C***KEYWORDS ERROR MESSAGE, XERROR 1271C***AUTHOR Fong, Kirby, (NMFECC at LLNL) 1272C***DESCRIPTION 1273C 1274C XERMSG processes a diagnostic message in a manner determined by the 1275C value of LEVEL and the current value of the library error control 1276C flag, KONTRL. See subroutine XSETF for details. 1277C 1278C LIBRAR A character constant (or character variable) with the name 1279C of the library. This will be 'SLATEC' for the SLATEC 1280C Common Math Library. The error handling package is 1281C general enough to be used by many libraries 1282C simultaneously, so it is desirable for the routine that 1283C detects and reports an error to identify the library name 1284C as well as the routine name. 1285C 1286C SUBROU A character constant (or character variable) with the name 1287C of the routine that detected the error. Usually it is the 1288C name of the routine that is calling XERMSG. There are 1289C some instances where a user callable library routine calls 1290C lower level subsidiary routines where the error is 1291C detected. In such cases it may be more informative to 1292C supply the name of the routine the user called rather than 1293C the name of the subsidiary routine that detected the 1294C error. 1295C 1296C MESSG A character constant (or character variable) with the text 1297C of the error or warning message. In the example below, 1298C the message is a character constant that contains a 1299C generic message. 1300C 1301C CALL XERMSG ('SLATEC', 'MMPY', 1302C *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION', 1303C *3, 1) 1304C 1305C It is possible (and is sometimes desirable) to generate a 1306C specific message--e.g., one that contains actual numeric 1307C values. Specific numeric values can be converted into 1308C character strings using formatted WRITE statements into 1309C character variables. This is called standard Fortran 1310C internal file I/O and is exemplified in the first three 1311C lines of the following example. You can also catenate 1312C substrings of characters to construct the error message. 1313C Here is an example showing the use of both writing to 1314C an internal file and catenating character strings. 1315C 1316C CHARACTER*5 CHARN, CHARL 1317C WRITE (CHARN,10) N 1318C WRITE (CHARL,10) LDA 1319C 10 FORMAT(I5) 1320C CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN// 1321C * ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'// 1322C * CHARL, 3, 1) 1323C 1324C There are two subtleties worth mentioning. One is that 1325C the // for character catenation is used to construct the 1326C error message so that no single character constant is 1327C continued to the next line. This avoids confusion as to 1328C whether there are trailing blanks at the end of the line. 1329C The second is that by catenating the parts of the message 1330C as an actual argument rather than encoding the entire 1331C message into one large character variable, we avoid 1332C having to know how long the message will be in order to 1333C declare an adequate length for that large character 1334C variable. XERMSG calls XERPRN to print the message using 1335C multiple lines if necessary. If the message is very long, 1336C XERPRN will break it into pieces of 72 characters (as 1337C requested by XERMSG) for printing on multiple lines. 1338C Also, XERMSG asks XERPRN to prefix each line with ' * ' 1339C so that the total line length could be 76 characters. 1340C Note also that XERPRN scans the error message backwards 1341C to ignore trailing blanks. Another feature is that 1342C the substring '$$' is treated as a new line sentinel 1343C by XERPRN. If you want to construct a multiline 1344C message without having to count out multiples of 72 1345C characters, just use '$$' as a separator. '$$' 1346C obviously must occur within 72 characters of the 1347C start of each line to have its intended effect since 1348C XERPRN is asked to wrap around at 72 characters in 1349C addition to looking for '$$'. 1350C 1351C NERR An integer value that is chosen by the library routine's 1352C author. It must be in the range -99 to 999 (three 1353C printable digits). Each distinct error should have its 1354C own error number. These error numbers should be described 1355C in the machine readable documentation for the routine. 1356C The error numbers need be unique only within each routine, 1357C so it is reasonable for each routine to start enumerating 1358C errors from 1 and proceeding to the next integer. 1359C 1360C LEVEL An integer value in the range 0 to 2 that indicates the 1361C level (severity) of the error. Their meanings are 1362C 1363C -1 A warning message. This is used if it is not clear 1364C that there really is an error, but the user's attention 1365C may be needed. An attempt is made to only print this 1366C message once. 1367C 1368C 0 A warning message. This is used if it is not clear 1369C that there really is an error, but the user's attention 1370C may be needed. 1371C 1372C 1 A recoverable error. This is used even if the error is 1373C so serious that the routine cannot return any useful 1374C answer. If the user has told the error package to 1375C return after recoverable errors, then XERMSG will 1376C return to the Library routine which can then return to 1377C the user's routine. The user may also permit the error 1378C package to terminate the program upon encountering a 1379C recoverable error. 1380C 1381C 2 A fatal error. XERMSG will not return to its caller 1382C after it receives a fatal error. This level should 1383C hardly ever be used; it is much better to allow the 1384C user a chance to recover. An example of one of the few 1385C cases in which it is permissible to declare a level 2 1386C error is a reverse communication Library routine that 1387C is likely to be called repeatedly until it integrates 1388C across some interval. If there is a serious error in 1389C the input such that another step cannot be taken and 1390C the Library routine is called again without the input 1391C error having been corrected by the caller, the Library 1392C routine will probably be called forever with improper 1393C input. In this case, it is reasonable to declare the 1394C error to be fatal. 1395C 1396C Each of the arguments to XERMSG is input; none will be modified by 1397C XERMSG. A routine may make multiple calls to XERMSG with warning 1398C level messages; however, after a call to XERMSG with a recoverable 1399C error, the routine should return to the user. Do not try to call 1400C XERMSG with a second recoverable error after the first recoverable 1401C error because the error package saves the error number. The user 1402C can retrieve this error number by calling another entry point in 1403C the error handling package and then clear the error number when 1404C recovering from the error. Calling XERMSG in succession causes the 1405C old error number to be overwritten by the latest error number. 1406C This is considered harmless for error numbers associated with 1407C warning messages but must not be done for error numbers of serious 1408C errors. After a call to XERMSG with a recoverable error, the user 1409C must be given a chance to call NUMXER or XERCLR to retrieve or 1410C clear the error number. 1411C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 1412C Error-handling Package, SAND82-0800, Sandia 1413C Laboratories, 1982. 1414C***ROUTINES CALLED FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE 1415C***REVISION HISTORY (YYMMDD) 1416C 880101 DATE WRITTEN 1417C 880621 REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988. 1418C THERE ARE TWO BASIC CHANGES. 1419C 1. A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO 1420C PRINT MESSAGES. THIS ROUTINE WILL BREAK LONG MESSAGES 1421C INTO PIECES FOR PRINTING ON MULTIPLE LINES. '$$' IS 1422C ACCEPTED AS A NEW LINE SENTINEL. A PREFIX CAN BE 1423C ADDED TO EACH LINE TO BE PRINTED. XERMSG USES EITHER 1424C ' ***' OR ' * ' AND LONG MESSAGES ARE BROKEN EVERY 1425C 72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE 1426C LENGTH OUTPUT CAN NOW BE AS GREAT AS 76. 1427C 2. THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE 1428C FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE 1429C OF LOWER CASE. 1430C 880708 REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30. 1431C THE PRINCIPAL CHANGES ARE 1432C 1. CLARIFY COMMENTS IN THE PROLOGUES 1433C 2. RENAME XRPRNT TO XERPRN 1434C 3. REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES 1435C SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE / 1436C CHARACTER FOR NEW RECORDS. 1437C 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO 1438C CLEAN UP THE CODING. 1439C 890721 REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN 1440C PREFIX. 1441C 891013 REVISED TO CORRECT COMMENTS. 1442C 891214 Prologue converted to Version 4.0 format. (WRB) 1443C 900510 Changed test on NERR to be -9999999 < NERR < 99999999, but 1444C NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3. Added 1445C LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and 1446C XERCTL to XERCNT. (RWC) 1447C 920501 Reformatted the REFERENCES section. (WRB) 1448C***END PROLOGUE XERMSG 1449 CHARACTER*(*) LIBRAR, SUBROU, MESSG 1450 CHARACTER*8 XLIBR, XSUBR 1451 CHARACTER*72 TEMP 1452 CHARACTER*20 LFIRST 1453C***FIRST EXECUTABLE STATEMENT XERMSG 1454 LKNTRL = J4SAVE (2, 0, .FALSE.) 1455 MAXMES = J4SAVE (4, 0, .FALSE.) 1456C 1457C LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL. 1458C MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE 1459C SHOULD BE PRINTED. 1460C 1461C WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN 1462C CALLING XERMSG. THE ERROR NUMBER SHOULD BE POSITIVE, 1463C AND THE LEVEL SHOULD BE BETWEEN 0 AND 2. 1464C 1465 IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR. 1466 * LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN 1467 CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' // 1468 * 'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '// 1469 * 'JOB ABORT DUE TO FATAL ERROR.', 72) 1470 CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY) 1471 CALL XERHLT (' ***XERMSG -- INVALID INPUT') 1472 RETURN 1473 ENDIF 1474C 1475C RECORD THE MESSAGE. 1476C 1477 I = J4SAVE (1, NERR, .TRUE.) 1478 CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT) 1479C 1480C HANDLE PRINT-ONCE WARNING MESSAGES. 1481C 1482 IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN 1483C 1484C ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG. 1485C 1486 XLIBR = LIBRAR 1487 XSUBR = SUBROU 1488 LFIRST = MESSG 1489 LERR = NERR 1490 LLEVEL = LEVEL 1491 CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL) 1492C 1493 LKNTRL = MAX(-2, MIN(2,LKNTRL)) 1494 MKNTRL = ABS(LKNTRL) 1495C 1496C SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS 1497C ZERO AND THE ERROR IS NOT FATAL. 1498C 1499 IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30 1500 IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30 1501 IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30 1502 IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30 1503C 1504C ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A 1505C MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS) 1506C AND SENDING IT OUT VIA XERPRN. PRINT ONLY IF CONTROL FLAG 1507C IS NOT ZERO. 1508C 1509 IF (LKNTRL .NE. 0) THEN 1510 TEMP(1:21) = 'MESSAGE FROM ROUTINE ' 1511 I = MIN(LEN(SUBROU), 16) 1512 TEMP(22:21+I) = SUBROU(1:I) 1513 TEMP(22+I:33+I) = ' IN LIBRARY ' 1514 LTEMP = 33 + I 1515 I = MIN(LEN(LIBRAR), 16) 1516 TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I) 1517 TEMP(LTEMP+I+1:LTEMP+I+1) = '.' 1518 LTEMP = LTEMP + I + 1 1519 CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72) 1520 ENDIF 1521C 1522C IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE 1523C PRINTING THE MESSAGE. THE INTRODUCTORY LINE TELLS THE CHOICE 1524C FROM EACH OF THE FOLLOWING THREE OPTIONS. 1525C 1. LEVEL OF THE MESSAGE 1526C 'INFORMATIVE MESSAGE' 1527C 'POTENTIALLY RECOVERABLE ERROR' 1528C 'FATAL ERROR' 1529C 2. WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE 1530C 'PROG CONTINUES' 1531C 'PROG ABORTED' 1532C 3. WHETHER OR NOT A TRACEBACK WAS REQUESTED. (THE TRACEBACK 1533C MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS 1534C WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.) 1535C 'TRACEBACK REQUESTED' 1536C 'TRACEBACK NOT REQUESTED' 1537C NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT 1538C EXCEED 74 CHARACTERS. 1539C WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED. 1540C 1541 IF (LKNTRL .GT. 0) THEN 1542C 1543C THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL. 1544C 1545 IF (LEVEL .LE. 0) THEN 1546 TEMP(1:20) = 'INFORMATIVE MESSAGE,' 1547 LTEMP = 20 1548 ELSEIF (LEVEL .EQ. 1) THEN 1549 TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,' 1550 LTEMP = 30 1551 ELSE 1552 TEMP(1:12) = 'FATAL ERROR,' 1553 LTEMP = 12 1554 ENDIF 1555C 1556C THEN WHETHER THE PROGRAM WILL CONTINUE. 1557C 1558 IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR. 1559 * (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN 1560 TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,' 1561 LTEMP = LTEMP + 14 1562 ELSE 1563 TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,' 1564 LTEMP = LTEMP + 16 1565 ENDIF 1566C 1567C FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK. 1568C 1569 IF (LKNTRL .GT. 0) THEN 1570 TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED' 1571 LTEMP = LTEMP + 20 1572 ELSE 1573 TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED' 1574 LTEMP = LTEMP + 24 1575 ENDIF 1576 CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72) 1577 ENDIF 1578C 1579C NOW SEND OUT THE MESSAGE. 1580C 1581 CALL XERPRN (' * ', -1, MESSG, 72) 1582C 1583C IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A 1584C TRACEBACK. 1585C 1586 IF (LKNTRL .GT. 0) THEN 1587 WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR 1588 DO 10 I=16,22 1589 IF (TEMP(I:I) .NE. ' ') GO TO 20 1590 10 CONTINUE 1591C 1592 20 CALL XERPRN (' * ', -1, TEMP(1:15) // TEMP(I:23), 72) 1593 CALL FDUMP 1594 ENDIF 1595C 1596C IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE. 1597C 1598 IF (LKNTRL .NE. 0) THEN 1599 CALL XERPRN (' * ', -1, ' ', 72) 1600 CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72) 1601 CALL XERPRN (' ', 0, ' ', 72) 1602 ENDIF 1603C 1604C IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE 1605C CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN. 1606C 1607 30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN 1608C 1609C THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A 1610C FATAL ERROR. PRINT THE REASON FOR THE ABORT AND THE ERROR 1611C SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT. 1612C 1613 IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN 1614 IF (LEVEL .EQ. 1) THEN 1615 CALL XERPRN 1616 * (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72) 1617 ELSE 1618 CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72) 1619 ENDIF 1620 CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY) 1621 CALL XERHLT (' ') 1622 ELSE 1623 CALL XERHLT (MESSG) 1624 ENDIF 1625 RETURN 1626 END 1627*DECK XERPRN 1628 SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP) 1629C***BEGIN PROLOGUE XERPRN 1630C***SUBSIDIARY 1631C***PURPOSE Print error messages processed by XERMSG. 1632C***LIBRARY SLATEC (XERROR) 1633C***CATEGORY R3C 1634C***TYPE ALL (XERPRN-A) 1635C***KEYWORDS ERROR MESSAGES, PRINTING, XERROR 1636C***AUTHOR Fong, Kirby, (NMFECC at LLNL) 1637C***DESCRIPTION 1638C 1639C This routine sends one or more lines to each of the (up to five) 1640C logical units to which error messages are to be sent. This routine 1641C is called several times by XERMSG, sometimes with a single line to 1642C print and sometimes with a (potentially very long) message that may 1643C wrap around into multiple lines. 1644C 1645C PREFIX Input argument of type CHARACTER. This argument contains 1646C characters to be put at the beginning of each line before 1647C the body of the message. No more than 16 characters of 1648C PREFIX will be used. 1649C 1650C NPREF Input argument of type INTEGER. This argument is the number 1651C of characters to use from PREFIX. If it is negative, the 1652C intrinsic function LEN is used to determine its length. If 1653C it is zero, PREFIX is not used. If it exceeds 16 or if 1654C LEN(PREFIX) exceeds 16, only the first 16 characters will be 1655C used. If NPREF is positive and the length of PREFIX is less 1656C than NPREF, a copy of PREFIX extended with blanks to length 1657C NPREF will be used. 1658C 1659C MESSG Input argument of type CHARACTER. This is the text of a 1660C message to be printed. If it is a long message, it will be 1661C broken into pieces for printing on multiple lines. Each line 1662C will start with the appropriate prefix and be followed by a 1663C piece of the message. NWRAP is the number of characters per 1664C piece; that is, after each NWRAP characters, we break and 1665C start a new line. In addition the characters '$$' embedded 1666C in MESSG are a sentinel for a new line. The counting of 1667C characters up to NWRAP starts over for each new line. The 1668C value of NWRAP typically used by XERMSG is 72 since many 1669C older error messages in the SLATEC Library are laid out to 1670C rely on wrap-around every 72 characters. 1671C 1672C NWRAP Input argument of type INTEGER. This gives the maximum size 1673C piece into which to break MESSG for printing on multiple 1674C lines. An embedded '$$' ends a line, and the count restarts 1675C at the following character. If a line break does not occur 1676C on a blank (it would split a word) that word is moved to the 1677C next line. Values of NWRAP less than 16 will be treated as 1678C 16. Values of NWRAP greater than 132 will be treated as 132. 1679C The actual line length will be NPREF + NWRAP after NPREF has 1680C been adjusted to fall between 0 and 16 and NWRAP has been 1681C adjusted to fall between 16 and 132. 1682C 1683C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 1684C Error-handling Package, SAND82-0800, Sandia 1685C Laboratories, 1982. 1686C***ROUTINES CALLED I1MACH, XGETUA 1687C***REVISION HISTORY (YYMMDD) 1688C 880621 DATE WRITTEN 1689C 880708 REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF 1690C JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK 1691C THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE 1692C SLASH CHARACTER IN FORMAT STATEMENTS. 1693C 890706 REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO 1694C STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK 1695C LINES TO BE PRINTED. 1696C 890721 REVISED TO ADD A NEW FEATURE. A NEGATIVE VALUE OF NPREF 1697C CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH. 1698C 891013 REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH. 1699C 891214 Prologue converted to Version 4.0 format. (WRB) 1700C 900510 Added code to break messages between words. (RWC) 1701C 920501 Reformatted the REFERENCES section. (WRB) 1702C***END PROLOGUE XERPRN 1703 CHARACTER*(*) PREFIX, MESSG 1704 INTEGER NPREF, NWRAP 1705 CHARACTER*148 CBUFF 1706 INTEGER IU(5), NUNIT 1707 CHARACTER*2 NEWLIN 1708 PARAMETER (NEWLIN = '$$') 1709C***FIRST EXECUTABLE STATEMENT XERPRN 1710 CALL XGETUA(IU,NUNIT) 1711C 1712C A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD 1713C ERROR MESSAGE UNIT INSTEAD. I1MACH(4) RETRIEVES THE STANDARD 1714C ERROR MESSAGE UNIT. 1715C 1716 N = I1MACH(4) 1717 DO 10 I=1,NUNIT 1718 IF (IU(I) .EQ. 0) IU(I) = N 1719 10 CONTINUE 1720C 1721C LPREF IS THE LENGTH OF THE PREFIX. THE PREFIX IS PLACED AT THE 1722C BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING 1723C THE REST OF THIS ROUTINE. 1724C 1725 IF ( NPREF .LT. 0 ) THEN 1726 LPREF = LEN(PREFIX) 1727 ELSE 1728 LPREF = NPREF 1729 ENDIF 1730 LPREF = MIN(16, LPREF) 1731 IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX 1732C 1733C LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE 1734C TIME FROM MESSG TO PRINT ON ONE LINE. 1735C 1736 LWRAP = MAX(16, MIN(132, NWRAP)) 1737C 1738C SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS. 1739C 1740 LENMSG = LEN(MESSG) 1741 N = LENMSG 1742 DO 20 I=1,N 1743 IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30 1744 LENMSG = LENMSG - 1 1745 20 CONTINUE 1746 30 CONTINUE 1747C 1748C IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE. 1749C 1750 IF (LENMSG .EQ. 0) THEN 1751 CBUFF(LPREF+1:LPREF+1) = ' ' 1752 DO 40 I=1,NUNIT 1753 WRITE(IU(I), '(A)') CBUFF(1:LPREF+1) 1754 40 CONTINUE 1755 RETURN 1756 ENDIF 1757C 1758C SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING 1759C STARTS. FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL. 1760C WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT. 1761C WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED. 1762C 1763C WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL. THE 1764C INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE 1765C OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH 1766C OF THE SECOND ARGUMENT. 1767C 1768C THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE 1769C FOLLOWING ORDER. WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER 1770C OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT 1771C POSITION NEXTC. 1772C 1773C LPIECE .EQ. 0 THE NEW LINE SENTINEL DOES NOT OCCUR IN THE 1774C REMAINDER OF THE CHARACTER STRING. LPIECE 1775C SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC, 1776C WHICHEVER IS LESS. 1777C 1778C LPIECE .EQ. 1 THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC: 1779C NEXTC). LPIECE IS EFFECTIVELY ZERO, AND WE 1780C PRINT NOTHING TO AVOID PRODUCING UNNECESSARY 1781C BLANK LINES. THIS TAKES CARE OF THE SITUATION 1782C WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF 1783C EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE 1784C SENTINEL FOLLOWED BY MORE CHARACTERS. NEXTC 1785C SHOULD BE INCREMENTED BY 2. 1786C 1787C LPIECE .GT. LWRAP+1 REDUCE LPIECE TO LWRAP. 1788C 1789C ELSE THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1 1790C RESET LPIECE = LPIECE-1. NOTE THAT THIS 1791C PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ. 1792C LWRAP+1. THAT IS, THE SENTINEL FALLS EXACTLY 1793C AT THE END OF A LINE. 1794C 1795 NEXTC = 1 1796 50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN) 1797 IF (LPIECE .EQ. 0) THEN 1798C 1799C THERE WAS NO NEW LINE SENTINEL FOUND. 1800C 1801 IDELTA = 0 1802 LPIECE = MIN(LWRAP, LENMSG+1-NEXTC) 1803 IF (LPIECE .LT. LENMSG+1-NEXTC) THEN 1804 DO 52 I=LPIECE+1,2,-1 1805 IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN 1806 LPIECE = I-1 1807 IDELTA = 1 1808 GOTO 54 1809 ENDIF 1810 52 CONTINUE 1811 ENDIF 1812 54 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1) 1813 NEXTC = NEXTC + LPIECE + IDELTA 1814 ELSEIF (LPIECE .EQ. 1) THEN 1815C 1816C WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1). 1817C DON'T PRINT A BLANK LINE. 1818C 1819 NEXTC = NEXTC + 2 1820 GO TO 50 1821 ELSEIF (LPIECE .GT. LWRAP+1) THEN 1822C 1823C LPIECE SHOULD BE SET DOWN TO LWRAP. 1824C 1825 IDELTA = 0 1826 LPIECE = LWRAP 1827 DO 56 I=LPIECE+1,2,-1 1828 IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN 1829 LPIECE = I-1 1830 IDELTA = 1 1831 GOTO 58 1832 ENDIF 1833 56 CONTINUE 1834 58 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1) 1835 NEXTC = NEXTC + LPIECE + IDELTA 1836 ELSE 1837C 1838C IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1. 1839C WE SHOULD DECREMENT LPIECE BY ONE. 1840C 1841 LPIECE = LPIECE - 1 1842 CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1) 1843 NEXTC = NEXTC + LPIECE + 2 1844 ENDIF 1845C 1846C PRINT 1847C 1848 DO 60 I=1,NUNIT 1849 WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE) 1850 60 CONTINUE 1851C 1852 IF (NEXTC .LE. LENMSG) GO TO 50 1853 RETURN 1854 END 1855*DECK XERSVE 1856 SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, 1857 + ICOUNT) 1858C***BEGIN PROLOGUE XERSVE 1859C***SUBSIDIARY 1860C***PURPOSE Record that an error has occurred. 1861C***LIBRARY SLATEC (XERROR) 1862C***CATEGORY R3 1863C***TYPE ALL (XERSVE-A) 1864C***KEYWORDS ERROR, XERROR 1865C***AUTHOR Jones, R. E., (SNLA) 1866C***DESCRIPTION 1867C 1868C *Usage: 1869C 1870C INTEGER KFLAG, NERR, LEVEL, ICOUNT 1871C CHARACTER * (len) LIBRAR, SUBROU, MESSG 1872C 1873C CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT) 1874C 1875C *Arguments: 1876C 1877C LIBRAR :IN is the library that the message is from. 1878C SUBROU :IN is the subroutine that the message is from. 1879C MESSG :IN is the message to be saved. 1880C KFLAG :IN indicates the action to be performed. 1881C when KFLAG > 0, the message in MESSG is saved. 1882C when KFLAG=0 the tables will be dumped and 1883C cleared. 1884C when KFLAG < 0, the tables will be dumped and 1885C not cleared. 1886C NERR :IN is the error number. 1887C LEVEL :IN is the error severity. 1888C ICOUNT :OUT the number of times this message has been seen, 1889C or zero if the table has overflowed and does not 1890C contain this message specifically. When KFLAG=0, 1891C ICOUNT will not be altered. 1892C 1893C *Description: 1894C 1895C Record that this error occurred and possibly dump and clear the 1896C tables. 1897C 1898C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 1899C Error-handling Package, SAND82-0800, Sandia 1900C Laboratories, 1982. 1901C***ROUTINES CALLED I1MACH, XGETUA 1902C***REVISION HISTORY (YYMMDD) 1903C 800319 DATE WRITTEN 1904C 861211 REVISION DATE from Version 3.2 1905C 891214 Prologue converted to Version 4.0 format. (BAB) 1906C 900413 Routine modified to remove reference to KFLAG. (WRB) 1907C 900510 Changed to add LIBRARY NAME and SUBROUTINE to calling 1908C sequence, use IF-THEN-ELSE, make number of saved entries 1909C easily changeable, changed routine name from XERSAV to 1910C XERSVE. (RWC) 1911C 910626 Added LIBTAB and SUBTAB to SAVE statement. (BKS) 1912C 920501 Reformatted the REFERENCES section. (WRB) 1913C***END PROLOGUE XERSVE 1914 PARAMETER (LENTAB=10) 1915 INTEGER LUN(5) 1916 CHARACTER*(*) LIBRAR, SUBROU, MESSG 1917 CHARACTER*8 LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB 1918 CHARACTER*20 MESTAB(LENTAB), MES 1919 DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB) 1920 SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG 1921 DATA KOUNTX/0/, NMSG/0/ 1922C***FIRST EXECUTABLE STATEMENT XERSVE 1923C 1924 IF (KFLAG.LE.0) THEN 1925C 1926C Dump the table. 1927C 1928 IF (NMSG.EQ.0) RETURN 1929C 1930C Print to each unit. 1931C 1932 CALL XGETUA (LUN, NUNIT) 1933 DO 20 KUNIT = 1,NUNIT 1934 IUNIT = LUN(KUNIT) 1935 IF (IUNIT.EQ.0) IUNIT = I1MACH(4) 1936C 1937C Print the table header. 1938C 1939 WRITE (IUNIT,9000) 1940C 1941C Print body of table. 1942C 1943 DO 10 I = 1,NMSG 1944 WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I), 1945 * NERTAB(I),LEVTAB(I),KOUNT(I) 1946 10 CONTINUE 1947C 1948C Print number of other errors. 1949C 1950 IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX 1951 WRITE (IUNIT,9030) 1952 20 CONTINUE 1953C 1954C Clear the error tables. 1955C 1956 IF (KFLAG.EQ.0) THEN 1957 NMSG = 0 1958 KOUNTX = 0 1959 ENDIF 1960 ELSE 1961C 1962C PROCESS A MESSAGE... 1963C SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, 1964C OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. 1965C 1966 LIB = LIBRAR 1967 SUB = SUBROU 1968 MES = MESSG 1969 DO 30 I = 1,NMSG 1970 IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND. 1971 * MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND. 1972 * LEVEL.EQ.LEVTAB(I)) THEN 1973 KOUNT(I) = KOUNT(I) + 1 1974 ICOUNT = KOUNT(I) 1975 RETURN 1976 ENDIF 1977 30 CONTINUE 1978C 1979 IF (NMSG.LT.LENTAB) THEN 1980C 1981C Empty slot found for new message. 1982C 1983 NMSG = NMSG + 1 1984 LIBTAB(I) = LIB 1985 SUBTAB(I) = SUB 1986 MESTAB(I) = MES 1987 NERTAB(I) = NERR 1988 LEVTAB(I) = LEVEL 1989 KOUNT (I) = 1 1990 ICOUNT = 1 1991 ELSE 1992C 1993C Table is full. 1994C 1995 KOUNTX = KOUNTX+1 1996 ICOUNT = 0 1997 ENDIF 1998 ENDIF 1999 RETURN 2000C 2001C Formats. 2002C 2003 9000 FORMAT ('0 ERROR MESSAGE SUMMARY' / 2004 + ' LIBRARY SUBROUTINE MESSAGE START NERR', 2005 + ' LEVEL COUNT') 2006 9010 FORMAT (1X,A,3X,A,3X,A,3I10) 2007 9020 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10) 2008 9030 FORMAT (1X) 2009 END 2010*DECK D1MACH 2011 DOUBLE PRECISION FUNCTION D1MACH (I) 2012C***BEGIN PROLOGUE D1MACH 2013C***PURPOSE Return floating point machine dependent constants. 2014C***LIBRARY SLATEC 2015C***CATEGORY R1 2016C***TYPE DOUBLE PRECISION (R1MACH-S, D1MACH-D) 2017C***KEYWORDS MACHINE CONSTANTS 2018C***AUTHOR Fox, P. A., (Bell Labs) 2019C Hall, A. D., (Bell Labs) 2020C Schryer, N. L., (Bell Labs) 2021C***DESCRIPTION 2022C 2023C D1MACH can be used to obtain machine-dependent parameters for the 2024C local machine environment. It is a function subprogram with one 2025C (input) argument, and can be referenced as follows: 2026C 2027C D = D1MACH(I) 2028C 2029C where I=1,...,5. The (output) value of D above is determined by 2030C the (input) value of I. The results for various values of I are 2031C discussed below. 2032C 2033C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. 2034C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. 2035C D1MACH( 3) = B**(-T), the smallest relative spacing. 2036C D1MACH( 4) = B**(1-T), the largest relative spacing. 2037C D1MACH( 5) = LOG10(B) 2038C 2039C Assume double precision numbers are represented in the T-digit, 2040C base-B form 2041C 2042C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) 2043C 2044C where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and 2045C EMIN .LE. E .LE. EMAX. 2046C 2047C The values of B, T, EMIN and EMAX are provided in I1MACH as 2048C follows: 2049C I1MACH(10) = B, the base. 2050C I1MACH(14) = T, the number of base-B digits. 2051C I1MACH(15) = EMIN, the smallest exponent E. 2052C I1MACH(16) = EMAX, the largest exponent E. 2053C 2054C To alter this function for a particular environment, the desired 2055C set of DATA statements should be activated by removing the C from 2056C column 1. Also, the values of D1MACH(1) - D1MACH(4) should be 2057C checked for consistency with the local operating system. 2058C 2059C***REFERENCES P. A. Fox, A. D. Hall and N. L. Schryer, Framework for 2060C a portable library, ACM Transactions on Mathematical 2061C Software 4, 2 (June 1978), pp. 177-188. 2062C***ROUTINES CALLED XERMSG 2063C***REVISION HISTORY (YYMMDD) 2064C 750101 DATE WRITTEN 2065C 890213 REVISION DATE from Version 3.2 2066C 891214 Prologue converted to Version 4.0 format. (BAB) 2067C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) 2068C 900618 Added DEC RISC constants. (WRB) 2069C 900723 Added IBM RS 6000 constants. (WRB) 2070C 900911 Added SUN 386i constants. (WRB) 2071C 910710 Added HP 730 constants. (SMR) 2072C 911114 Added Convex IEEE constants. (WRB) 2073C 920121 Added SUN -r8 compiler option constants. (WRB) 2074C 920229 Added Touchstone Delta i860 constants. (WRB) 2075C 920501 Reformatted the REFERENCES section. (WRB) 2076C 920625 Added CONVEX -p8 and -pd8 compiler option constants. 2077C (BKS, WRB) 2078C 930201 Added DEC Alpha and SGI constants. (RWC and WRB) 2079C 010817 Elevated IEEE to highest importance; see next set of 2080C comments below. (DWL) 2081C***END PROLOGUE D1MACH 2082C 2083 INTEGER SMALL(4) 2084 INTEGER LARGE(4) 2085 INTEGER RIGHT(4) 2086 INTEGER DIVER(4) 2087 INTEGER LOG10(4) 2088C 2089C Initial data here correspond to the IEEE standard. The values for 2090C DMACH(1), DMACH(3) and DMACH(4) are slight upper bounds. The value 2091C for DMACH(2) is a slight lower bound. The value for DMACH(5) is 2092C a 20-digit approximation. If one of the sets of initial data below 2093C is preferred, do the necessary commenting and uncommenting. (DWL) 2094 DOUBLE PRECISION DMACH(5) 2095 DATA DMACH / 2.23D-308, 1.79D+308, 1.111D-16, 2.222D-16, 2096 1 0.30102999566398119521D0 / 2097 SAVE DMACH 2098C 2099cc EQUIVALENCE (DMACH(1),SMALL(1)) 2100cc EQUIVALENCE (DMACH(2),LARGE(1)) 2101cc EQUIVALENCE (DMACH(3),RIGHT(1)) 2102cc EQUIVALENCE (DMACH(4),DIVER(1)) 2103cc EQUIVALENCE (DMACH(5),LOG10(1)) 2104C 2105C MACHINE CONSTANTS FOR THE AMIGA 2106C ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION 2107C 2108C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / 2109C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / 2110C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / 2111C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / 2112C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / 2113C 2114C MACHINE CONSTANTS FOR THE AMIGA 2115C ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT 2116C 2117C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / 2118C DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' / 2119C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / 2120C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / 2121C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / 2122C 2123C MACHINE CONSTANTS FOR THE APOLLO 2124C 2125C DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 / 2126C DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF / 2127C DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 / 2128C DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 / 2129C DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF / 2130C 2131C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM 2132C 2133C DATA SMALL(1) / ZC00800000 / 2134C DATA SMALL(2) / Z000000000 / 2135C DATA LARGE(1) / ZDFFFFFFFF / 2136C DATA LARGE(2) / ZFFFFFFFFF / 2137C DATA RIGHT(1) / ZCC5800000 / 2138C DATA RIGHT(2) / Z000000000 / 2139C DATA DIVER(1) / ZCC6800000 / 2140C DATA DIVER(2) / Z000000000 / 2141C DATA LOG10(1) / ZD00E730E7 / 2142C DATA LOG10(2) / ZC77800DC0 / 2143C 2144C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM 2145C 2146C DATA SMALL(1) / O1771000000000000 / 2147C DATA SMALL(2) / O0000000000000000 / 2148C DATA LARGE(1) / O0777777777777777 / 2149C DATA LARGE(2) / O0007777777777777 / 2150C DATA RIGHT(1) / O1461000000000000 / 2151C DATA RIGHT(2) / O0000000000000000 / 2152C DATA DIVER(1) / O1451000000000000 / 2153C DATA DIVER(2) / O0000000000000000 / 2154C DATA LOG10(1) / O1157163034761674 / 2155C DATA LOG10(2) / O0006677466732724 / 2156C 2157C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS 2158C 2159C DATA SMALL(1) / O1771000000000000 / 2160C DATA SMALL(2) / O7770000000000000 / 2161C DATA LARGE(1) / O0777777777777777 / 2162C DATA LARGE(2) / O7777777777777777 / 2163C DATA RIGHT(1) / O1461000000000000 / 2164C DATA RIGHT(2) / O0000000000000000 / 2165C DATA DIVER(1) / O1451000000000000 / 2166C DATA DIVER(2) / O0000000000000000 / 2167C DATA LOG10(1) / O1157163034761674 / 2168C DATA LOG10(2) / O0006677466732724 / 2169C 2170C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE 2171C 2172C DATA SMALL(1) / Z"3001800000000000" / 2173C DATA SMALL(2) / Z"3001000000000000" / 2174C DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" / 2175C DATA LARGE(2) / Z"4FFE000000000000" / 2176C DATA RIGHT(1) / Z"3FD2800000000000" / 2177C DATA RIGHT(2) / Z"3FD2000000000000" / 2178C DATA DIVER(1) / Z"3FD3800000000000" / 2179C DATA DIVER(2) / Z"3FD3000000000000" / 2180C DATA LOG10(1) / Z"3FFF9A209A84FBCF" / 2181C DATA LOG10(2) / Z"3FFFF7988F8959AC" / 2182C 2183C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES 2184C 2185C DATA SMALL(1) / 00564000000000000000B / 2186C DATA SMALL(2) / 00000000000000000000B / 2187C DATA LARGE(1) / 37757777777777777777B / 2188C DATA LARGE(2) / 37157777777777777777B / 2189C DATA RIGHT(1) / 15624000000000000000B / 2190C DATA RIGHT(2) / 00000000000000000000B / 2191C DATA DIVER(1) / 15634000000000000000B / 2192C DATA DIVER(2) / 00000000000000000000B / 2193C DATA LOG10(1) / 17164642023241175717B / 2194C DATA LOG10(2) / 16367571421742254654B / 2195C 2196C MACHINE CONSTANTS FOR THE CELERITY C1260 2197C 2198C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / 2199C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / 2200C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / 2201C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / 2202C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / 2203C 2204C MACHINE CONSTANTS FOR THE CONVEX 2205C USING THE -fn OR -pd8 COMPILER OPTION 2206C 2207C DATA DMACH(1) / Z'0010000000000000' / 2208C DATA DMACH(2) / Z'7FFFFFFFFFFFFFFF' / 2209C DATA DMACH(3) / Z'3CC0000000000000' / 2210C DATA DMACH(4) / Z'3CD0000000000000' / 2211C DATA DMACH(5) / Z'3FF34413509F79FF' / 2212C 2213C MACHINE CONSTANTS FOR THE CONVEX 2214C USING THE -fi COMPILER OPTION 2215C 2216C DATA DMACH(1) / Z'0010000000000000' / 2217C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / 2218C DATA DMACH(3) / Z'3CA0000000000000' / 2219C DATA DMACH(4) / Z'3CB0000000000000' / 2220C DATA DMACH(5) / Z'3FD34413509F79FF' / 2221C 2222C MACHINE CONSTANTS FOR THE CONVEX 2223C USING THE -p8 COMPILER OPTION 2224C 2225C DATA DMACH(1) / Z'00010000000000000000000000000000' / 2226C DATA DMACH(2) / Z'7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF' / 2227C DATA DMACH(3) / Z'3F900000000000000000000000000000' / 2228C DATA DMACH(4) / Z'3F910000000000000000000000000000' / 2229C DATA DMACH(5) / Z'3FFF34413509F79FEF311F12B35816F9' / 2230C 2231C MACHINE CONSTANTS FOR THE CRAY 2232C 2233C DATA SMALL(1) / 201354000000000000000B / 2234C DATA SMALL(2) / 000000000000000000000B / 2235C DATA LARGE(1) / 577767777777777777777B / 2236C DATA LARGE(2) / 000007777777777777774B / 2237C DATA RIGHT(1) / 376434000000000000000B / 2238C DATA RIGHT(2) / 000000000000000000000B / 2239C DATA DIVER(1) / 376444000000000000000B / 2240C DATA DIVER(2) / 000000000000000000000B / 2241C DATA LOG10(1) / 377774642023241175717B / 2242C DATA LOG10(2) / 000007571421742254654B / 2243C 2244C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 2245C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - 2246C STATIC DMACH(5) 2247C 2248C DATA SMALL / 20K, 3*0 / 2249C DATA LARGE / 77777K, 3*177777K / 2250C DATA RIGHT / 31420K, 3*0 / 2251C DATA DIVER / 32020K, 3*0 / 2252C DATA LOG10 / 40423K, 42023K, 50237K, 74776K / 2253C 2254C MACHINE CONSTANTS FOR THE DEC ALPHA 2255C USING G_FLOAT 2256C 2257C DATA DMACH(1) / '0000000000000010'X / 2258C DATA DMACH(2) / 'FFFFFFFFFFFF7FFF'X / 2259C DATA DMACH(3) / '0000000000003CC0'X / 2260C DATA DMACH(4) / '0000000000003CD0'X / 2261C DATA DMACH(5) / '79FF509F44133FF3'X / 2262C 2263C MACHINE CONSTANTS FOR THE DEC ALPHA 2264C USING IEEE_FORMAT 2265C 2266C DATA DMACH(1) / '0010000000000000'X / 2267C DATA DMACH(2) / '7FEFFFFFFFFFFFFF'X / 2268C DATA DMACH(3) / '3CA0000000000000'X / 2269C DATA DMACH(4) / '3CB0000000000000'X / 2270C DATA DMACH(5) / '3FD34413509F79FF'X / 2271C 2272C MACHINE CONSTANTS FOR THE DEC RISC 2273C 2274C DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000'/ 2275C DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF'/ 2276C DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000'/ 2277C DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000'/ 2278C DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413'/ 2279C 2280C MACHINE CONSTANTS FOR THE DEC VAX 2281C USING D_FLOATING 2282C (EXPRESSED IN INTEGER AND HEXADECIMAL) 2283C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS 2284C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS 2285C 2286C DATA SMALL(1), SMALL(2) / 128, 0 / 2287C DATA LARGE(1), LARGE(2) / -32769, -1 / 2288C DATA RIGHT(1), RIGHT(2) / 9344, 0 / 2289C DATA DIVER(1), DIVER(2) / 9472, 0 / 2290C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / 2291C 2292C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / 2293C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / 2294C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / 2295C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / 2296C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / 2297C 2298C MACHINE CONSTANTS FOR THE DEC VAX 2299C USING G_FLOATING 2300C (EXPRESSED IN INTEGER AND HEXADECIMAL) 2301C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS 2302C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS 2303C 2304C DATA SMALL(1), SMALL(2) / 16, 0 / 2305C DATA LARGE(1), LARGE(2) / -32769, -1 / 2306C DATA RIGHT(1), RIGHT(2) / 15552, 0 / 2307C DATA DIVER(1), DIVER(2) / 15568, 0 / 2308C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / 2309C 2310C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / 2311C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / 2312C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / 2313C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / 2314C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / 2315C 2316C MACHINE CONSTANTS FOR THE ELXSI 6400 2317C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) 2318C 2319C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / 2320C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / 2321C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / 2322C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / 2323C DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X / 2324C 2325C MACHINE CONSTANTS FOR THE HARRIS 220 2326C 2327C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / 2328C DATA LARGE(1), LARGE(2) / '37777777, '37777577 / 2329C DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 / 2330C DATA DIVER(1), DIVER(2) / '20000000, '00000334 / 2331C DATA LOG10(1), LOG10(2) / '23210115, '10237777 / 2332C 2333C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES 2334C 2335C DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 / 2336C DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 / 2337C DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 / 2338C DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 / 2339C DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 / 2340C 2341C MACHINE CONSTANTS FOR THE HP 730 2342C 2343C DATA DMACH(1) / Z'0010000000000000' / 2344C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / 2345C DATA DMACH(3) / Z'3CA0000000000000' / 2346C DATA DMACH(4) / Z'3CB0000000000000' / 2347C DATA DMACH(5) / Z'3FD34413509F79FF' / 2348C 2349C MACHINE CONSTANTS FOR THE HP 2100 2350C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 2351C 2352C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / 2353C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / 2354C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / 2355C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / 2356C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / 2357C 2358C MACHINE CONSTANTS FOR THE HP 2100 2359C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 2360C 2361C DATA SMALL(1), SMALL(2) / 40000B, 0 / 2362C DATA SMALL(3), SMALL(4) / 0, 1 / 2363C DATA LARGE(1), LARGE(2) / 77777B, 177777B / 2364C DATA LARGE(3), LARGE(4) / 177777B, 177776B / 2365C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / 2366C DATA RIGHT(3), RIGHT(4) / 0, 225B / 2367C DATA DIVER(1), DIVER(2) / 40000B, 0 / 2368C DATA DIVER(3), DIVER(4) / 0, 227B / 2369C DATA LOG10(1), LOG10(2) / 46420B, 46502B / 2370C DATA LOG10(3), LOG10(4) / 76747B, 176377B / 2371C 2372C MACHINE CONSTANTS FOR THE HP 9000 2373C 2374C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / 2375C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / 2376C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / 2377C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / 2378C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / 2379C 2380C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, 2381C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND 2382C THE PERKIN ELMER (INTERDATA) 7/32. 2383C 2384C DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 / 2385C DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / 2386C DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 / 2387C DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 / 2388C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / 2389C 2390C MACHINE CONSTANTS FOR THE IBM PC 2391C ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION 2392C ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087. 2393C 2394C DATA SMALL(1) / 2.23D-308 / 2395C DATA LARGE(1) / 1.79D+308 / 2396C DATA RIGHT(1) / 1.11D-16 / 2397C DATA DIVER(1) / 2.22D-16 / 2398C DATA LOG10(1) / 0.301029995663981195D0 / 2399C 2400C MACHINE CONSTANTS FOR THE IBM RS 6000 2401C 2402C DATA DMACH(1) / Z'0010000000000000' / 2403C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / 2404C DATA DMACH(3) / Z'3CA0000000000000' / 2405C DATA DMACH(4) / Z'3CB0000000000000' / 2406C DATA DMACH(5) / Z'3FD34413509F79FF' / 2407C 2408C MACHINE CONSTANTS FOR THE INTEL i860 2409C 2410C DATA DMACH(1) / Z'0010000000000000' / 2411C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / 2412C DATA DMACH(3) / Z'3CA0000000000000' / 2413C DATA DMACH(4) / Z'3CB0000000000000' / 2414C DATA DMACH(5) / Z'3FD34413509F79FF' / 2415C 2416C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR) 2417C 2418C DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 / 2419C DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 / 2420C DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 / 2421C DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 / 2422C DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 / 2423C 2424C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR) 2425C 2426C DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 / 2427C DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 / 2428C DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 / 2429C DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 / 2430C DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 / 2431C 2432C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING 2433C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). 2434C 2435C DATA SMALL(1), SMALL(2) / 8388608, 0 / 2436C DATA LARGE(1), LARGE(2) / 2147483647, -1 / 2437C DATA RIGHT(1), RIGHT(2) / 612368384, 0 / 2438C DATA DIVER(1), DIVER(2) / 620756992, 0 / 2439C DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 / 2440C 2441C DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 / 2442C DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 / 2443C DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 / 2444C DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 / 2445C DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 / 2446C 2447C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING 2448C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). 2449C 2450C DATA SMALL(1), SMALL(2) / 128, 0 / 2451C DATA SMALL(3), SMALL(4) / 0, 0 / 2452C DATA LARGE(1), LARGE(2) / 32767, -1 / 2453C DATA LARGE(3), LARGE(4) / -1, -1 / 2454C DATA RIGHT(1), RIGHT(2) / 9344, 0 / 2455C DATA RIGHT(3), RIGHT(4) / 0, 0 / 2456C DATA DIVER(1), DIVER(2) / 9472, 0 / 2457C DATA DIVER(3), DIVER(4) / 0, 0 / 2458C DATA LOG10(1), LOG10(2) / 16282, 8346 / 2459C DATA LOG10(3), LOG10(4) / -31493, -12296 / 2460C 2461C DATA SMALL(1), SMALL(2) / O000200, O000000 / 2462C DATA SMALL(3), SMALL(4) / O000000, O000000 / 2463C DATA LARGE(1), LARGE(2) / O077777, O177777 / 2464C DATA LARGE(3), LARGE(4) / O177777, O177777 / 2465C DATA RIGHT(1), RIGHT(2) / O022200, O000000 / 2466C DATA RIGHT(3), RIGHT(4) / O000000, O000000 / 2467C DATA DIVER(1), DIVER(2) / O022400, O000000 / 2468C DATA DIVER(3), DIVER(4) / O000000, O000000 / 2469C DATA LOG10(1), LOG10(2) / O037632, O020232 / 2470C DATA LOG10(3), LOG10(4) / O102373, O147770 / 2471C 2472C MACHINE CONSTANTS FOR THE SILICON GRAPHICS 2473C 2474C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / 2475C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / 2476C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / 2477C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / 2478C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / 2479C 2480C MACHINE CONSTANTS FOR THE SUN 2481C 2482C DATA DMACH(1) / Z'0010000000000000' / 2483C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / 2484C DATA DMACH(3) / Z'3CA0000000000000' / 2485C DATA DMACH(4) / Z'3CB0000000000000' / 2486C DATA DMACH(5) / Z'3FD34413509F79FF' / 2487C 2488C MACHINE CONSTANTS FOR THE SUN 2489C USING THE -r8 COMPILER OPTION 2490C 2491C DATA DMACH(1) / Z'00010000000000000000000000000000' / 2492C DATA DMACH(2) / Z'7FFEFFFFFFFFFFFFFFFFFFFFFFFFFFFF' / 2493C DATA DMACH(3) / Z'3F8E0000000000000000000000000000' / 2494C DATA DMACH(4) / Z'3F8F0000000000000000000000000000' / 2495C DATA DMACH(5) / Z'3FFD34413509F79FEF311F12B35816F9' / 2496C 2497C MACHINE CONSTANTS FOR THE SUN 386i 2498C 2499C DATA SMALL(1), SMALL(2) / Z'FFFFFFFD', Z'000FFFFF' / 2500C DATA LARGE(1), LARGE(2) / Z'FFFFFFB0', Z'7FEFFFFF' / 2501C DATA RIGHT(1), RIGHT(2) / Z'000000B0', Z'3CA00000' / 2502C DATA DIVER(1), DIVER(2) / Z'FFFFFFCB', Z'3CAFFFFF' 2503C DATA LOG10(1), LOG10(2) / Z'509F79E9', Z'3FD34413' / 2504C 2505C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER 2506C 2507C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / 2508C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / 2509C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / 2510C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / 2511C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 / 2512C 2513C***FIRST EXECUTABLE STATEMENT D1MACH 2514 IF (I .LT. 1 .OR. I .GT. 5) CALL XERMSG ('SLATEC', 'D1MACH', 2515 + 'I OUT OF BOUNDS', 1, 2) 2516C 2517 D1MACH = DMACH(I) 2518 RETURN 2519C 2520 END 2521*DECK XGETUA 2522 SUBROUTINE XGETUA (IUNITA, N) 2523C***BEGIN PROLOGUE XGETUA 2524C***PURPOSE Return unit number(s) to which error messages are being 2525C sent. 2526C***LIBRARY SLATEC (XERROR) 2527C***CATEGORY R3C 2528C***TYPE ALL (XGETUA-A) 2529C***KEYWORDS ERROR, XERROR 2530C***AUTHOR Jones, R. E., (SNLA) 2531C***DESCRIPTION 2532C 2533C Abstract 2534C XGETUA may be called to determine the unit number or numbers 2535C to which error messages are being sent. 2536C These unit numbers may have been set by a call to XSETUN, 2537C or a call to XSETUA, or may be a default value. 2538C 2539C Description of Parameters 2540C --Output-- 2541C IUNIT - an array of one to five unit numbers, depending 2542C on the value of N. A value of zero refers to the 2543C default unit, as defined by the I1MACH machine 2544C constant routine. Only IUNIT(1),...,IUNIT(N) are 2545C defined by XGETUA. The values of IUNIT(N+1),..., 2546C IUNIT(5) are not defined (for N .LT. 5) or altered 2547C in any way by XGETUA. 2548C N - the number of units to which copies of the 2549C error messages are being sent. N will be in the 2550C range from 1 to 5. 2551C 2552C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 2553C Error-handling Package, SAND82-0800, Sandia 2554C Laboratories, 1982. 2555C***ROUTINES CALLED J4SAVE 2556C***REVISION HISTORY (YYMMDD) 2557C 790801 DATE WRITTEN 2558C 861211 REVISION DATE from Version 3.2 2559C 891214 Prologue converted to Version 4.0 format. (BAB) 2560C 920501 Reformatted the REFERENCES section. (WRB) 2561C***END PROLOGUE XGETUA 2562 DIMENSION IUNITA(5) 2563C***FIRST EXECUTABLE STATEMENT XGETUA 2564 N = J4SAVE(5,0,.FALSE.) 2565 DO 30 I=1,N 2566 INDEX = I+4 2567 IF (I.EQ.1) INDEX = 3 2568 IUNITA(I) = J4SAVE(INDEX,0,.FALSE.) 2569 30 CONTINUE 2570 RETURN 2571 END 2572*DECK DSTEPS 2573 SUBROUTINE DSTEPS (DF, NEQN, Y, X, H, EPS, WT, START, HOLD, K, 2574 + KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, V, W, G, 2575 + PHASE1, NS, NORND, KSTEPS, TWOU, FOURU, XOLD, KPREV, IVC, IV, 2576 + KGI, GI, RPAR, IPAR) 2577C***BEGIN PROLOGUE DSTEPS 2578C***PURPOSE Integrate a system of first order ordinary differential 2579C equations one step. 2580C***LIBRARY SLATEC (DEPAC) 2581C***CATEGORY I1A1B 2582C***TYPE DOUBLE PRECISION (STEPS-S, DSTEPS-D) 2583C***KEYWORDS ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE, 2584C ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR 2585C***AUTHOR Shampine, L. F., (SNLA) 2586C Gordon, M. K., (SNLA) 2587C MODIFIED BY H.A. WATTS 2588C***DESCRIPTION 2589C 2590C Written by L. F. Shampine and M. K. Gordon 2591C 2592C Abstract 2593C 2594C Subroutine DSTEPS is normally used indirectly through subroutine 2595C DDEABM . Because DDEABM suffices for most problems and is much 2596C easier to use, using it should be considered before using DSTEPS 2597C alone. 2598C 2599C Subroutine DSTEPS integrates a system of NEQN first order ordinary 2600C differential equations one step, normally from X to X+H, using a 2601C modified divided difference form of the Adams Pece formulas. Local 2602C extrapolation is used to improve absolute stability and accuracy. 2603C The code adjusts its order and step size to control the local error 2604C per unit step in a generalized sense. Special devices are included 2605C to control roundoff error and to detect when the user is requesting 2606C too much accuracy. 2607C 2608C This code is completely explained and documented in the text, 2609C Computer Solution of Ordinary Differential Equations, The Initial 2610C Value Problem by L. F. Shampine and M. K. Gordon. 2611C Further details on use of this code are available in "Solving 2612C Ordinary Differential Equations with ODE, STEP, and INTRP", 2613C by L. F. Shampine and M. K. Gordon, SLA-73-1060. 2614C 2615C 2616C The parameters represent -- 2617C DF -- subroutine to evaluate derivatives 2618C NEQN -- number of equations to be integrated 2619C Y(*) -- solution vector at X 2620C X -- independent variable 2621C H -- appropriate step size for next step. Normally determined by 2622C code 2623C EPS -- local error tolerance 2624C WT(*) -- vector of weights for error criterion 2625C START -- logical variable set .TRUE. for first step, .FALSE. 2626C otherwise 2627C HOLD -- step size used for last successful step 2628C K -- appropriate order for next step (determined by code) 2629C KOLD -- order used for last successful step 2630C CRASH -- logical variable set .TRUE. when no step can be taken, 2631C .FALSE. otherwise. 2632C YP(*) -- derivative of solution vector at X after successful 2633C step 2634C KSTEPS -- counter on attempted steps 2635C TWOU -- 2.*U where U is machine unit roundoff quantity 2636C FOURU -- 4.*U where U is machine unit roundoff quantity 2637C RPAR,IPAR -- parameter arrays which you may choose to use 2638C for communication between your program and subroutine F. 2639C They are not altered or used by DSTEPS. 2640C The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G, 2641C W,P,IV and GI are required for the interpolation subroutine SINTRP. 2642C The remaining variables and arrays are included in the call list 2643C only to eliminate local retention of variables between calls. 2644C 2645C Input to DSTEPS 2646C 2647C First call -- 2648C 2649C The user must provide storage in his calling program for all arrays 2650C in the call list, namely 2651C 2652C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12), 2653C 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10), 2654C 2 RPAR(*),IPAR(*) 2655C 2656C **Note** 2657C 2658C The user must also declare START , CRASH , PHASE1 and NORND 2659C logical variables and DF an EXTERNAL subroutine, supply the 2660C subroutine DF(X,Y,YP) to evaluate 2661C DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN)) 2662C and initialize only the following parameters. 2663C NEQN -- number of equations to be integrated 2664C Y(*) -- vector of initial values of dependent variables 2665C X -- initial value of the independent variable 2666C H -- nominal step size indicating direction of integration 2667C and maximum size of step. Must be variable 2668C EPS -- local error tolerance per step. Must be variable 2669C WT(*) -- vector of non-zero weights for error criterion 2670C START -- .TRUE. 2671C YP(*) -- vector of initial derivative values 2672C KSTEPS -- set KSTEPS to zero 2673C TWOU -- 2.*U where U is machine unit roundoff quantity 2674C FOURU -- 4.*U where U is machine unit roundoff quantity 2675C Define U to be the machine unit roundoff quantity by calling 2676C the function routine D1MACH, U = D1MACH(4), or by 2677C computing U so that U is the smallest positive number such 2678C that 1.0+U .GT. 1.0. 2679C 2680C DSTEPS requires that the L2 norm of the vector with components 2681C LOCAL ERROR(L)/WT(L) be less than EPS for a successful step. The 2682C array WT allows the user to specify an error test appropriate 2683C for his problem. For example, 2684C WT(L) = 1.0 specifies absolute error, 2685C = ABS(Y(L)) error relative to the most recent value of the 2686C L-th component of the solution, 2687C = ABS(YP(L)) error relative to the most recent value of 2688C the L-th component of the derivative, 2689C = MAX(WT(L),ABS(Y(L))) error relative to the largest 2690C magnitude of L-th component obtained so far, 2691C = ABS(Y(L))*RELERR/EPS + ABSERR/EPS specifies a mixed 2692C relative-absolute test where RELERR is relative 2693C error, ABSERR is absolute error and EPS = 2694C MAX(RELERR,ABSERR) . 2695C 2696C Subsequent calls -- 2697C 2698C Subroutine DSTEPS is designed so that all information needed to 2699C continue the integration, including the step size H and the order 2700C K , is returned with each step. With the exception of the step 2701C size, the error tolerance, and the weights, none of the parameters 2702C should be altered. The array WT must be updated after each step 2703C to maintain relative error tests like those above. Normally the 2704C integration is continued just beyond the desired endpoint and the 2705C solution interpolated there with subroutine SINTRP . If it is 2706C impossible to integrate beyond the endpoint, the step size may be 2707C reduced to hit the endpoint since the code will not take a step 2708C larger than the H input. Changing the direction of integration, 2709C i.e., the sign of H , requires the user set START = .TRUE. before 2710C calling DSTEPS again. This is the only situation in which START 2711C should be altered. 2712C 2713C Output from DSTEPS 2714C 2715C Successful Step -- 2716C 2717C The subroutine returns after each successful step with START and 2718C CRASH set .FALSE. . X represents the independent variable 2719C advanced one step of length HOLD from its value on input and Y 2720C the solution vector at the new value of X . All other parameters 2721C represent information corresponding to the new X needed to 2722C continue the integration. 2723C 2724C Unsuccessful Step -- 2725C 2726C When the error tolerance is too small for the machine precision, 2727C the subroutine returns without taking a step and CRASH = .TRUE. . 2728C An appropriate step size and error tolerance for continuing are 2729C estimated and all other information is restored as upon input 2730C before returning. To continue with the larger tolerance, the user 2731C just calls the code again. A restart is neither required nor 2732C desirable. 2733C 2734C***REFERENCES L. F. Shampine and M. K. Gordon, Solving ordinary 2735C differential equations with ODE, STEP, and INTRP, 2736C Report SLA-73-1060, Sandia Laboratories, 1973. 2737C***ROUTINES CALLED D1MACH, DHSTRT 2738C***REVISION HISTORY (YYMMDD) 2739C 740101 DATE WRITTEN 2740C 890531 Changed all specific intrinsics to generic. (WRB) 2741C 890831 Modified array declarations. (WRB) 2742C 890831 REVISION DATE from Version 3.2 2743C 891214 Prologue converted to Version 4.0 format. (BAB) 2744C 920501 Reformatted the REFERENCES section. (WRB) 2745C***END PROLOGUE DSTEPS 2746C 2747 INTEGER I, IFAIL, IM1, IP1, IPAR, IQ, J, K, KM1, KM2, KNEW, 2748 1 KOLD, KP1, KP2, KSTEPS, L, LIMIT1, LIMIT2, NEQN, NS, NSM2, 2749 2 NSP1, NSP2 2750 DOUBLE PRECISION ABSH, ALPHA, BETA, BIG, D1MACH, 2751 1 EPS, ERK, ERKM1, ERKM2, ERKP1, ERR, 2752 2 FOURU, G, GI, GSTR, H, HNEW, HOLD, P, P5EPS, PHI, PSI, R, 2753 3 REALI, REALNS, RHO, ROUND, RPAR, SIG, TAU, TEMP1, 2754 4 TEMP2, TEMP3, TEMP4, TEMP5, TEMP6, TWO, TWOU, U, V, W, WT, 2755 5 X, XOLD, Y, YP 2756 LOGICAL START,CRASH,PHASE1,NORND 2757 DIMENSION Y(*),WT(*),PHI(NEQN,16),P(*),YP(*),PSI(12), 2758 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10), 2759 2 RPAR(*),IPAR(*) 2760 DIMENSION TWO(13),GSTR(13) 2761 EXTERNAL DF 2762 SAVE TWO, GSTR 2763C 2764 DATA TWO(1),TWO(2),TWO(3),TWO(4),TWO(5),TWO(6),TWO(7),TWO(8), 2765 1 TWO(9),TWO(10),TWO(11),TWO(12),TWO(13) 2766 2 /2.0D0,4.0D0,8.0D0,16.0D0,32.0D0,64.0D0,128.0D0,256.0D0, 2767 3 512.0D0,1024.0D0,2048.0D0,4096.0D0,8192.0D0/ 2768 DATA GSTR(1),GSTR(2),GSTR(3),GSTR(4),GSTR(5),GSTR(6),GSTR(7), 2769 1 GSTR(8),GSTR(9),GSTR(10),GSTR(11),GSTR(12),GSTR(13) 2770 2 /0.5D0,0.0833D0,0.0417D0,0.0264D0,0.0188D0,0.0143D0,0.0114D0, 2771 3 0.00936D0,0.00789D0,0.00679D0,0.00592D0,0.00524D0,0.00468D0/ 2772C 2773C *** BEGIN BLOCK 0 *** 2774C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE 2775C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A 2776C STARTING STEP SIZE. 2777C *** 2778C 2779C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE 2780C 2781C***FIRST EXECUTABLE STATEMENT DSTEPS 2782 CRASH = .TRUE. 2783 IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5 2784 H = SIGN(FOURU*ABS(X),H) 2785 RETURN 2786 5 P5EPS = 0.5D0*EPS 2787C 2788C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE 2789C 2790 ROUND = 0.0D0 2791 DO 10 L = 1,NEQN 2792 10 ROUND = ROUND + (Y(L)/WT(L))**2 2793 ROUND = TWOU*SQRT(ROUND) 2794 IF(P5EPS .GE. ROUND) GO TO 15 2795 EPS = 2.0D0*ROUND*(1.0D0 + FOURU) 2796 RETURN 2797 15 CRASH = .FALSE. 2798 G(1) = 1.0D0 2799 G(2) = 0.5D0 2800 SIG(1) = 1.0D0 2801 IF(.NOT.START) GO TO 99 2802C 2803C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP 2804C 2805C CALL DF(X,Y,YP,RPAR,IPAR) 2806C SUM = 0.0 2807 DO 20 L = 1,NEQN 2808 PHI(L,1) = YP(L) 2809 20 PHI(L,2) = 0.0D0 2810C20 SUM = SUM + (YP(L)/WT(L))**2 2811C SUM = SQRT(SUM) 2812C ABSH = ABS(H) 2813C IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM) 2814C H = SIGN(MAX(ABSH,FOURU*ABS(X)),H) 2815C 2816 U = D1MACH(4) 2817 BIG = SQRT(D1MACH(2)) 2818 CALL DHSTRT(DF,NEQN,X,X+H,Y,YP,WT,1,U,BIG, 2819 1 PHI(1,3),PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H) 2820C 2821 HOLD = 0.0D0 2822 K = 1 2823 KOLD = 0 2824 KPREV = 0 2825 START = .FALSE. 2826 PHASE1 = .TRUE. 2827 NORND = .TRUE. 2828 IF(P5EPS .GT. 100.0D0*ROUND) GO TO 99 2829 NORND = .FALSE. 2830 DO 25 L = 1,NEQN 2831 25 PHI(L,15) = 0.0D0 2832 99 IFAIL = 0 2833C *** END BLOCK 0 *** 2834C 2835C *** BEGIN BLOCK 1 *** 2836C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING 2837C THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED. 2838C *** 2839C 2840 100 KP1 = K+1 2841 KP2 = K+2 2842 KM1 = K-1 2843 KM2 = K-2 2844C 2845C NS IS THE NUMBER OF DSTEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT 2846C ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE 2847C 2848 IF(H .NE. HOLD) NS = 0 2849 IF (NS.LE.KOLD) NS = NS+1 2850 NSP1 = NS+1 2851 IF (K .LT. NS) GO TO 199 2852C 2853C COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH 2854C ARE CHANGED 2855C 2856 BETA(NS) = 1.0D0 2857 REALNS = NS 2858 ALPHA(NS) = 1.0D0/REALNS 2859 TEMP1 = H*REALNS 2860 SIG(NSP1) = 1.0D0 2861 IF(K .LT. NSP1) GO TO 110 2862 DO 105 I = NSP1,K 2863 IM1 = I-1 2864 TEMP2 = PSI(IM1) 2865 PSI(IM1) = TEMP1 2866 BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2 2867 TEMP1 = TEMP2 + H 2868 ALPHA(I) = H/TEMP1 2869 REALI = I 2870 105 SIG(I+1) = REALI*ALPHA(I)*SIG(I) 2871 110 PSI(K) = TEMP1 2872C 2873C COMPUTE COEFFICIENTS G(*) 2874C 2875C INITIALIZE V(*) AND SET W(*). 2876C 2877 IF(NS .GT. 1) GO TO 120 2878 DO 115 IQ = 1,K 2879 TEMP3 = IQ*(IQ+1) 2880 V(IQ) = 1.0D0/TEMP3 2881 115 W(IQ) = V(IQ) 2882 IVC = 0 2883 KGI = 0 2884 IF (K .EQ. 1) GO TO 140 2885 KGI = 1 2886 GI(1) = W(2) 2887 GO TO 140 2888C 2889C IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*) 2890C 2891 120 IF(K .LE. KPREV) GO TO 130 2892 IF (IVC .EQ. 0) GO TO 122 2893 JV = KP1 - IV(IVC) 2894 IVC = IVC - 1 2895 GO TO 123 2896 122 JV = 1 2897 TEMP4 = K*KP1 2898 V(K) = 1.0D0/TEMP4 2899 W(K) = V(K) 2900 IF (K .NE. 2) GO TO 123 2901 KGI = 1 2902 GI(1) = W(2) 2903 123 NSM2 = NS-2 2904 IF(NSM2 .LT. JV) GO TO 130 2905 DO 125 J = JV,NSM2 2906 I = K-J 2907 V(I) = V(I) - ALPHA(J+1)*V(I+1) 2908 125 W(I) = V(I) 2909 IF (I .NE. 2) GO TO 130 2910 KGI = NS - 1 2911 GI(KGI) = W(2) 2912C 2913C UPDATE V(*) AND SET W(*) 2914C 2915 130 LIMIT1 = KP1 - NS 2916 TEMP5 = ALPHA(NS) 2917 DO 135 IQ = 1,LIMIT1 2918 V(IQ) = V(IQ) - TEMP5*V(IQ+1) 2919 135 W(IQ) = V(IQ) 2920 G(NSP1) = W(1) 2921 IF (LIMIT1 .EQ. 1) GO TO 137 2922 KGI = NS 2923 GI(KGI) = W(2) 2924 137 W(LIMIT1+1) = V(LIMIT1+1) 2925 IF (K .GE. KOLD) GO TO 140 2926 IVC = IVC + 1 2927 IV(IVC) = LIMIT1 + 2 2928C 2929C COMPUTE THE G(*) IN THE WORK VECTOR W(*) 2930C 2931 140 NSP2 = NS + 2 2932 KPREV = K 2933 IF(KP1 .LT. NSP2) GO TO 199 2934 DO 150 I = NSP2,KP1 2935 LIMIT2 = KP2 - I 2936 TEMP6 = ALPHA(I-1) 2937 DO 145 IQ = 1,LIMIT2 2938 145 W(IQ) = W(IQ) - TEMP6*W(IQ+1) 2939 150 G(I) = W(1) 2940 199 CONTINUE 2941C *** END BLOCK 1 *** 2942C 2943C *** BEGIN BLOCK 2 *** 2944C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED 2945C SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K, 2946C K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED. 2947C *** 2948C 2949C INCREMENT COUNTER ON ATTEMPTED DSTEPS 2950C 2951 KSTEPS = KSTEPS + 1 2952C 2953C CHANGE PHI TO PHI STAR 2954C 2955 IF(K .LT. NSP1) GO TO 215 2956 DO 210 I = NSP1,K 2957 TEMP1 = BETA(I) 2958 DO 205 L = 1,NEQN 2959 205 PHI(L,I) = TEMP1*PHI(L,I) 2960 210 CONTINUE 2961C 2962C PREDICT SOLUTION AND DIFFERENCES 2963C 2964 215 DO 220 L = 1,NEQN 2965 PHI(L,KP2) = PHI(L,KP1) 2966 PHI(L,KP1) = 0.0D0 2967 220 P(L) = 0.0D0 2968 DO 230 J = 1,K 2969 I = KP1 - J 2970 IP1 = I+1 2971 TEMP2 = G(I) 2972 DO 225 L = 1,NEQN 2973 P(L) = P(L) + TEMP2*PHI(L,I) 2974 225 PHI(L,I) = PHI(L,I) + PHI(L,IP1) 2975 230 CONTINUE 2976 IF(NORND) GO TO 240 2977 DO 235 L = 1,NEQN 2978 TAU = H*P(L) - PHI(L,15) 2979 P(L) = Y(L) + TAU 2980 235 PHI(L,16) = (P(L) - Y(L)) - TAU 2981 GO TO 250 2982 240 DO 245 L = 1,NEQN 2983 245 P(L) = Y(L) + H*P(L) 2984 250 XOLD = X 2985 X = X + H 2986 ABSH = ABS(H) 2987 CALL DF(X,P,YP,RPAR,IPAR) 2988C 2989C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 2990C 2991 ERKM2 = 0.0D0 2992 ERKM1 = 0.0D0 2993 ERK = 0.0D0 2994 DO 265 L = 1,NEQN 2995 TEMP3 = 1.0D0/WT(L) 2996 TEMP4 = YP(L) - PHI(L,1) 2997 IF(KM2)265,260,255 2998 255 ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2 2999 260 ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2 3000 265 ERK = ERK + (TEMP4*TEMP3)**2 3001 IF(KM2)280,275,270 3002 270 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2) 3003 275 ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1) 3004 280 TEMP5 = ABSH*SQRT(ERK) 3005 ERR = TEMP5*(G(K)-G(KP1)) 3006 ERK = TEMP5*SIG(KP1)*GSTR(K) 3007 KNEW = K 3008C 3009C TEST IF ORDER SHOULD BE LOWERED 3010C 3011 IF(KM2)299,290,285 3012 285 IF(MAX(ERKM1,ERKM2) .LE. ERK) KNEW = KM1 3013 GO TO 299 3014 290 IF(ERKM1 .LE. 0.5D0*ERK) KNEW = KM1 3015C 3016C TEST IF STEP SUCCESSFUL 3017C 3018 299 IF(ERR .LE. EPS) GO TO 400 3019C *** END BLOCK 2 *** 3020C 3021C *** BEGIN BLOCK 3 *** 3022C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) . 3023C IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE 3024C THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR 3025C TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE 3026C PRECISION. 3027C *** 3028C 3029C RESTORE X, PHI(*,*) AND PSI(*) 3030C 3031 PHASE1 = .FALSE. 3032 X = XOLD 3033 DO 310 I = 1,K 3034 TEMP1 = 1.0D0/BETA(I) 3035 IP1 = I+1 3036 DO 305 L = 1,NEQN 3037 305 PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1)) 3038 310 CONTINUE 3039 IF(K .LT. 2) GO TO 320 3040 DO 315 I = 2,K 3041 315 PSI(I-1) = PSI(I) - H 3042C 3043C ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, USE OPTIMAL STEP 3044C SIZE 3045C 3046 320 IFAIL = IFAIL + 1 3047 TEMP2 = 0.5D0 3048 IF(IFAIL - 3) 335,330,325 3049 325 IF(P5EPS .LT. 0.25D0*ERK) TEMP2 = SQRT(P5EPS/ERK) 3050 330 KNEW = 1 3051 335 H = TEMP2*H 3052 K = KNEW 3053 NS = 0 3054 IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340 3055 CRASH = .TRUE. 3056 H = SIGN(FOURU*ABS(X),H) 3057 EPS = EPS + EPS 3058 RETURN 3059 340 GO TO 100 3060C *** END BLOCK 3 *** 3061C 3062C *** BEGIN BLOCK 4 *** 3063C THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION, EVALUATE 3064C THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE 3065C DIFFERENCES. DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP. 3066C *** 3067 400 KOLD = K 3068 HOLD = H 3069C 3070C CORRECT AND EVALUATE 3071C 3072 TEMP1 = H*G(KP1) 3073 IF(NORND) GO TO 410 3074 DO 405 L = 1,NEQN 3075 TEMP3 = Y(L) 3076 RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16) 3077 Y(L) = P(L) + RHO 3078 PHI(L,15) = (Y(L) - P(L)) - RHO 3079 405 P(L) = TEMP3 3080 GO TO 420 3081 410 DO 415 L = 1,NEQN 3082 TEMP3 = Y(L) 3083 Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1)) 3084 415 P(L) = TEMP3 3085 420 CALL DF(X,Y,YP,RPAR,IPAR) 3086C 3087C UPDATE DIFFERENCES FOR NEXT STEP 3088C 3089 DO 425 L = 1,NEQN 3090 PHI(L,KP1) = YP(L) - PHI(L,1) 3091 425 PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2) 3092 DO 435 I = 1,K 3093 DO 430 L = 1,NEQN 3094 430 PHI(L,I) = PHI(L,I) + PHI(L,KP1) 3095 435 CONTINUE 3096C 3097C ESTIMATE ERROR AT ORDER K+1 UNLESS: 3098C IN FIRST PHASE WHEN ALWAYS RAISE ORDER, 3099C ALREADY DECIDED TO LOWER ORDER, 3100C STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE 3101C 3102 ERKP1 = 0.0D0 3103 IF(KNEW .EQ. KM1 .OR. K .EQ. 12) PHASE1 = .FALSE. 3104 IF(PHASE1) GO TO 450 3105 IF(KNEW .EQ. KM1) GO TO 455 3106 IF(KP1 .GT. NS) GO TO 460 3107 DO 440 L = 1,NEQN 3108 440 ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2 3109 ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1) 3110C 3111C USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER 3112C FOR NEXT STEP 3113C 3114 IF(K .GT. 1) GO TO 445 3115 IF(ERKP1 .GE. 0.5D0*ERK) GO TO 460 3116 GO TO 450 3117 445 IF(ERKM1 .LE. MIN(ERK,ERKP1)) GO TO 455 3118 IF(ERKP1 .GE. ERK .OR. K .EQ. 12) GO TO 460 3119C 3120C HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE 3121C BEEN LOWERED IN BLOCK 2. THUS ORDER IS TO BE RAISED 3122C 3123C RAISE ORDER 3124C 3125 450 K = KP1 3126 ERK = ERKP1 3127 GO TO 460 3128C 3129C LOWER ORDER 3130C 3131 455 K = KM1 3132 ERK = ERKM1 3133C 3134C WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP 3135C 3136 460 HNEW = H + H 3137 IF(PHASE1) GO TO 465 3138 IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465 3139 HNEW = H 3140 IF(P5EPS .GE. ERK) GO TO 465 3141 TEMP2 = K+1 3142 R = (P5EPS/ERK)**(1.0D0/TEMP2) 3143 HNEW = ABSH*MAX(0.5D0,MIN(0.9D0,R)) 3144 HNEW = SIGN(MAX(HNEW,FOURU*ABS(X)),H) 3145 465 H = HNEW 3146 RETURN 3147C *** END BLOCK 4 *** 3148 END 3149*DECK FDUMP 3150 SUBROUTINE FDUMP 3151C***BEGIN PROLOGUE FDUMP 3152C***PURPOSE Symbolic dump (should be locally written). 3153C***LIBRARY SLATEC (XERROR) 3154C***CATEGORY R3 3155C***TYPE ALL (FDUMP-A) 3156C***KEYWORDS ERROR, XERMSG 3157C***AUTHOR Jones, R. E., (SNLA) 3158C***DESCRIPTION 3159C 3160C ***Note*** Machine Dependent Routine 3161C FDUMP is intended to be replaced by a locally written 3162C version which produces a symbolic dump. Failing this, 3163C it should be replaced by a version which prints the 3164C subprogram nesting list. Note that this dump must be 3165C printed on each of up to five files, as indicated by the 3166C XGETUA routine. See XSETUA and XGETUA for details. 3167C 3168C Written by Ron Jones, with SLATEC Common Math Library Subcommittee 3169C 3170C***REFERENCES (NONE) 3171C***ROUTINES CALLED (NONE) 3172C***REVISION HISTORY (YYMMDD) 3173C 790801 DATE WRITTEN 3174C 861211 REVISION DATE from Version 3.2 3175C 891214 Prologue converted to Version 4.0 format. (BAB) 3176C***END PROLOGUE FDUMP 3177C***FIRST EXECUTABLE STATEMENT FDUMP 3178 RETURN 3179 END 3180*DECK I1MACH 3181 INTEGER FUNCTION I1MACH (I) 3182C***BEGIN PROLOGUE I1MACH 3183C***PURPOSE Return integer machine dependent constants. 3184C***LIBRARY SLATEC 3185C***CATEGORY R1 3186C***TYPE INTEGER (I1MACH-I) 3187C***KEYWORDS MACHINE CONSTANTS 3188C***AUTHOR Fox, P. A., (Bell Labs) 3189C Hall, A. D., (Bell Labs) 3190C Schryer, N. L., (Bell Labs) 3191C***DESCRIPTION 3192C 3193C I1MACH can be used to obtain machine-dependent parameters for the 3194C local machine environment. It is a function subprogram with one 3195C (input) argument and can be referenced as follows: 3196C 3197C K = I1MACH(I) 3198C 3199C where I=1,...,16. The (output) value of K above is determined by 3200C the (input) value of I. The results for various values of I are 3201C discussed below. 3202C 3203C I/O unit numbers: 3204C I1MACH( 1) = the standard input unit. 3205C I1MACH( 2) = the standard output unit. 3206C I1MACH( 3) = the standard punch unit. 3207C I1MACH( 4) = the standard error message unit. 3208C 3209C Words: 3210C I1MACH( 5) = the number of bits per integer storage unit. 3211C I1MACH( 6) = the number of characters per integer storage unit. 3212C 3213C Integers: 3214C assume integers are represented in the S-digit, base-A form 3215C 3216C sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) 3217C 3218C where 0 .LE. X(I) .LT. A for I=0,...,S-1. 3219C I1MACH( 7) = A, the base. 3220C I1MACH( 8) = S, the number of base-A digits. 3221C I1MACH( 9) = A**S - 1, the largest magnitude. 3222C 3223C Floating-Point Numbers: 3224C Assume floating-point numbers are represented in the T-digit, 3225C base-B form 3226C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) 3227C 3228C where 0 .LE. X(I) .LT. B for I=1,...,T, 3229C 0 .LT. X(1), and EMIN .LE. E .LE. EMAX. 3230C I1MACH(10) = B, the base. 3231C 3232C Single-Precision: 3233C I1MACH(11) = T, the number of base-B digits. 3234C I1MACH(12) = EMIN, the smallest exponent E. 3235C I1MACH(13) = EMAX, the largest exponent E. 3236C 3237C Double-Precision: 3238C I1MACH(14) = T, the number of base-B digits. 3239C I1MACH(15) = EMIN, the smallest exponent E. 3240C I1MACH(16) = EMAX, the largest exponent E. 3241C 3242C To alter this function for a particular environment, the desired 3243C set of DATA statements should be activated by removing the C from 3244C column 1. Also, the values of I1MACH(1) - I1MACH(4) should be 3245C checked for consistency with the local operating system. 3246C 3247C***REFERENCES P. A. Fox, A. D. Hall and N. L. Schryer, Framework for 3248C a portable library, ACM Transactions on Mathematical 3249C Software 4, 2 (June 1978), pp. 177-188. 3250C***ROUTINES CALLED (NONE) 3251C***REVISION HISTORY (YYMMDD) 3252C 750101 DATE WRITTEN 3253C 891012 Added VAX G-floating constants. (WRB) 3254C 891012 REVISION DATE from Version 3.2 3255C 891214 Prologue converted to Version 4.0 format. (BAB) 3256C 900618 Added DEC RISC constants. (WRB) 3257C 900723 Added IBM RS 6000 constants. (WRB) 3258C 901009 Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16. 3259C (RWC) 3260C 910710 Added HP 730 constants. (SMR) 3261C 911114 Added Convex IEEE constants. (WRB) 3262C 920121 Added SUN -r8 compiler option constants. (WRB) 3263C 920229 Added Touchstone Delta i860 constants. (WRB) 3264C 920501 Reformatted the REFERENCES section. (WRB) 3265C 920625 Added Convex -p8 and -pd8 compiler option constants. 3266C (BKS, WRB) 3267C 930201 Added DEC Alpha and SGI constants. (RWC and WRB) 3268C 930618 Corrected I1MACH(5) for Convex -p8 and -pd8 compiler 3269C options. (DWL, RWC and WRB). 3270C 010817 Elevated IEEE to highest importance; see next set of 3271C comments below. (DWL) 3272C***END PROLOGUE I1MACH 3273C 3274C Initial data here correspond to the IEEE standard. If one of the 3275C sets of initial data below is preferred, do the necessary commenting 3276C and uncommenting. (DWL) 3277 INTEGER IMACH(16),OUTPUT 3278 DATA IMACH( 1) / 5 / 3279 DATA IMACH( 2) / 6 / 3280 DATA IMACH( 3) / 6 / 3281 DATA IMACH( 4) / 6 / 3282 DATA IMACH( 5) / 32 / 3283 DATA IMACH( 6) / 4 / 3284 DATA IMACH( 7) / 2 / 3285 DATA IMACH( 8) / 31 / 3286 DATA IMACH( 9) / 2147483647 / 3287 DATA IMACH(10) / 2 / 3288 DATA IMACH(11) / 24 / 3289 DATA IMACH(12) / -126 / 3290 DATA IMACH(13) / 127 / 3291 DATA IMACH(14) / 53 / 3292 DATA IMACH(15) / -1022 / 3293 DATA IMACH(16) / 1023 / 3294 SAVE IMACH 3295cc EQUIVALENCE (IMACH(4),OUTPUT) 3296C 3297C MACHINE CONSTANTS FOR THE AMIGA 3298C ABSOFT COMPILER 3299C 3300C DATA IMACH( 1) / 5 / 3301C DATA IMACH( 2) / 6 / 3302C DATA IMACH( 3) / 5 / 3303C DATA IMACH( 4) / 6 / 3304C DATA IMACH( 5) / 32 / 3305C DATA IMACH( 6) / 4 / 3306C DATA IMACH( 7) / 2 / 3307C DATA IMACH( 8) / 31 / 3308C DATA IMACH( 9) / 2147483647 / 3309C DATA IMACH(10) / 2 / 3310C DATA IMACH(11) / 24 / 3311C DATA IMACH(12) / -126 / 3312C DATA IMACH(13) / 127 / 3313C DATA IMACH(14) / 53 / 3314C DATA IMACH(15) / -1022 / 3315C DATA IMACH(16) / 1023 / 3316C 3317C MACHINE CONSTANTS FOR THE APOLLO 3318C 3319C DATA IMACH( 1) / 5 / 3320C DATA IMACH( 2) / 6 / 3321C DATA IMACH( 3) / 6 / 3322C DATA IMACH( 4) / 6 / 3323C DATA IMACH( 5) / 32 / 3324C DATA IMACH( 6) / 4 / 3325C DATA IMACH( 7) / 2 / 3326C DATA IMACH( 8) / 31 / 3327C DATA IMACH( 9) / 2147483647 / 3328C DATA IMACH(10) / 2 / 3329C DATA IMACH(11) / 24 / 3330C DATA IMACH(12) / -125 / 3331C DATA IMACH(13) / 129 / 3332C DATA IMACH(14) / 53 / 3333C DATA IMACH(15) / -1021 / 3334C DATA IMACH(16) / 1025 / 3335C 3336C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM 3337C 3338C DATA IMACH( 1) / 7 / 3339C DATA IMACH( 2) / 2 / 3340C DATA IMACH( 3) / 2 / 3341C DATA IMACH( 4) / 2 / 3342C DATA IMACH( 5) / 36 / 3343C DATA IMACH( 6) / 4 / 3344C DATA IMACH( 7) / 2 / 3345C DATA IMACH( 8) / 33 / 3346C DATA IMACH( 9) / Z1FFFFFFFF / 3347C DATA IMACH(10) / 2 / 3348C DATA IMACH(11) / 24 / 3349C DATA IMACH(12) / -256 / 3350C DATA IMACH(13) / 255 / 3351C DATA IMACH(14) / 60 / 3352C DATA IMACH(15) / -256 / 3353C DATA IMACH(16) / 255 / 3354C 3355C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM 3356C 3357C DATA IMACH( 1) / 5 / 3358C DATA IMACH( 2) / 6 / 3359C DATA IMACH( 3) / 7 / 3360C DATA IMACH( 4) / 6 / 3361C DATA IMACH( 5) / 48 / 3362C DATA IMACH( 6) / 6 / 3363C DATA IMACH( 7) / 2 / 3364C DATA IMACH( 8) / 39 / 3365C DATA IMACH( 9) / O0007777777777777 / 3366C DATA IMACH(10) / 8 / 3367C DATA IMACH(11) / 13 / 3368C DATA IMACH(12) / -50 / 3369C DATA IMACH(13) / 76 / 3370C DATA IMACH(14) / 26 / 3371C DATA IMACH(15) / -50 / 3372C DATA IMACH(16) / 76 / 3373C 3374C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS 3375C 3376C DATA IMACH( 1) / 5 / 3377C DATA IMACH( 2) / 6 / 3378C DATA IMACH( 3) / 7 / 3379C DATA IMACH( 4) / 6 / 3380C DATA IMACH( 5) / 48 / 3381C DATA IMACH( 6) / 6 / 3382C DATA IMACH( 7) / 2 / 3383C DATA IMACH( 8) / 39 / 3384C DATA IMACH( 9) / O0007777777777777 / 3385C DATA IMACH(10) / 8 / 3386C DATA IMACH(11) / 13 / 3387C DATA IMACH(12) / -50 / 3388C DATA IMACH(13) / 76 / 3389C DATA IMACH(14) / 26 / 3390C DATA IMACH(15) / -32754 / 3391C DATA IMACH(16) / 32780 / 3392C 3393C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE 3394C 3395C DATA IMACH( 1) / 5 / 3396C DATA IMACH( 2) / 6 / 3397C DATA IMACH( 3) / 7 / 3398C DATA IMACH( 4) / 6 / 3399C DATA IMACH( 5) / 64 / 3400C DATA IMACH( 6) / 8 / 3401C DATA IMACH( 7) / 2 / 3402C DATA IMACH( 8) / 63 / 3403C DATA IMACH( 9) / 9223372036854775807 / 3404C DATA IMACH(10) / 2 / 3405C DATA IMACH(11) / 47 / 3406C DATA IMACH(12) / -4095 / 3407C DATA IMACH(13) / 4094 / 3408C DATA IMACH(14) / 94 / 3409C DATA IMACH(15) / -4095 / 3410C DATA IMACH(16) / 4094 / 3411C 3412C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES 3413C 3414C DATA IMACH( 1) / 5 / 3415C DATA IMACH( 2) / 6 / 3416C DATA IMACH( 3) / 7 / 3417C DATA IMACH( 4) / 6LOUTPUT/ 3418C DATA IMACH( 5) / 60 / 3419C DATA IMACH( 6) / 10 / 3420C DATA IMACH( 7) / 2 / 3421C DATA IMACH( 8) / 48 / 3422C DATA IMACH( 9) / 00007777777777777777B / 3423C DATA IMACH(10) / 2 / 3424C DATA IMACH(11) / 47 / 3425C DATA IMACH(12) / -929 / 3426C DATA IMACH(13) / 1070 / 3427C DATA IMACH(14) / 94 / 3428C DATA IMACH(15) / -929 / 3429C DATA IMACH(16) / 1069 / 3430C 3431C MACHINE CONSTANTS FOR THE CELERITY C1260 3432C 3433C DATA IMACH( 1) / 5 / 3434C DATA IMACH( 2) / 6 / 3435C DATA IMACH( 3) / 6 / 3436C DATA IMACH( 4) / 0 / 3437C DATA IMACH( 5) / 32 / 3438C DATA IMACH( 6) / 4 / 3439C DATA IMACH( 7) / 2 / 3440C DATA IMACH( 8) / 31 / 3441C DATA IMACH( 9) / Z'7FFFFFFF' / 3442C DATA IMACH(10) / 2 / 3443C DATA IMACH(11) / 24 / 3444C DATA IMACH(12) / -126 / 3445C DATA IMACH(13) / 127 / 3446C DATA IMACH(14) / 53 / 3447C DATA IMACH(15) / -1022 / 3448C DATA IMACH(16) / 1023 / 3449C 3450C MACHINE CONSTANTS FOR THE CONVEX 3451C USING THE -fn COMPILER OPTION 3452C 3453C DATA IMACH( 1) / 5 / 3454C DATA IMACH( 2) / 6 / 3455C DATA IMACH( 3) / 7 / 3456C DATA IMACH( 4) / 6 / 3457C DATA IMACH( 5) / 32 / 3458C DATA IMACH( 6) / 4 / 3459C DATA IMACH( 7) / 2 / 3460C DATA IMACH( 8) / 31 / 3461C DATA IMACH( 9) / 2147483647 / 3462C DATA IMACH(10) / 2 / 3463C DATA IMACH(11) / 24 / 3464C DATA IMACH(12) / -127 / 3465C DATA IMACH(13) / 127 / 3466C DATA IMACH(14) / 53 / 3467C DATA IMACH(15) / -1023 / 3468C DATA IMACH(16) / 1023 / 3469C 3470C MACHINE CONSTANTS FOR THE CONVEX 3471C USING THE -fi COMPILER OPTION 3472C 3473C DATA IMACH( 1) / 5 / 3474C DATA IMACH( 2) / 6 / 3475C DATA IMACH( 3) / 7 / 3476C DATA IMACH( 4) / 6 / 3477C DATA IMACH( 5) / 32 / 3478C DATA IMACH( 6) / 4 / 3479C DATA IMACH( 7) / 2 / 3480C DATA IMACH( 8) / 31 / 3481C DATA IMACH( 9) / 2147483647 / 3482C DATA IMACH(10) / 2 / 3483C DATA IMACH(11) / 24 / 3484C DATA IMACH(12) / -125 / 3485C DATA IMACH(13) / 128 / 3486C DATA IMACH(14) / 53 / 3487C DATA IMACH(15) / -1021 / 3488C DATA IMACH(16) / 1024 / 3489C 3490C MACHINE CONSTANTS FOR THE CONVEX 3491C USING THE -p8 COMPILER OPTION 3492C 3493C DATA IMACH( 1) / 5 / 3494C DATA IMACH( 2) / 6 / 3495C DATA IMACH( 3) / 7 / 3496C DATA IMACH( 4) / 6 / 3497C DATA IMACH( 5) / 64 / 3498C DATA IMACH( 6) / 4 / 3499C DATA IMACH( 7) / 2 / 3500C DATA IMACH( 8) / 63 / 3501C DATA IMACH( 9) / 9223372036854775807 / 3502C DATA IMACH(10) / 2 / 3503C DATA IMACH(11) / 53 / 3504C DATA IMACH(12) / -1023 / 3505C DATA IMACH(13) / 1023 / 3506C DATA IMACH(14) / 113 / 3507C DATA IMACH(15) / -16383 / 3508C DATA IMACH(16) / 16383 / 3509C 3510C MACHINE CONSTANTS FOR THE CONVEX 3511C USING THE -pd8 COMPILER OPTION 3512C 3513C DATA IMACH( 1) / 5 / 3514C DATA IMACH( 2) / 6 / 3515C DATA IMACH( 3) / 7 / 3516C DATA IMACH( 4) / 6 / 3517C DATA IMACH( 5) / 64 / 3518C DATA IMACH( 6) / 4 / 3519C DATA IMACH( 7) / 2 / 3520C DATA IMACH( 8) / 63 / 3521C DATA IMACH( 9) / 9223372036854775807 / 3522C DATA IMACH(10) / 2 / 3523C DATA IMACH(11) / 53 / 3524C DATA IMACH(12) / -1023 / 3525C DATA IMACH(13) / 1023 / 3526C DATA IMACH(14) / 53 / 3527C DATA IMACH(15) / -1023 / 3528C DATA IMACH(16) / 1023 / 3529C 3530C MACHINE CONSTANTS FOR THE CRAY 3531C USING THE 46 BIT INTEGER COMPILER OPTION 3532C 3533C DATA IMACH( 1) / 100 / 3534C DATA IMACH( 2) / 101 / 3535C DATA IMACH( 3) / 102 / 3536C DATA IMACH( 4) / 101 / 3537C DATA IMACH( 5) / 64 / 3538C DATA IMACH( 6) / 8 / 3539C DATA IMACH( 7) / 2 / 3540C DATA IMACH( 8) / 46 / 3541C DATA IMACH( 9) / 1777777777777777B / 3542C DATA IMACH(10) / 2 / 3543C DATA IMACH(11) / 47 / 3544C DATA IMACH(12) / -8189 / 3545C DATA IMACH(13) / 8190 / 3546C DATA IMACH(14) / 94 / 3547C DATA IMACH(15) / -8099 / 3548C DATA IMACH(16) / 8190 / 3549C 3550C MACHINE CONSTANTS FOR THE CRAY 3551C USING THE 64 BIT INTEGER COMPILER OPTION 3552C 3553C DATA IMACH( 1) / 100 / 3554C DATA IMACH( 2) / 101 / 3555C DATA IMACH( 3) / 102 / 3556C DATA IMACH( 4) / 101 / 3557C DATA IMACH( 5) / 64 / 3558C DATA IMACH( 6) / 8 / 3559C DATA IMACH( 7) / 2 / 3560C DATA IMACH( 8) / 63 / 3561C DATA IMACH( 9) / 777777777777777777777B / 3562C DATA IMACH(10) / 2 / 3563C DATA IMACH(11) / 47 / 3564C DATA IMACH(12) / -8189 / 3565C DATA IMACH(13) / 8190 / 3566C DATA IMACH(14) / 94 / 3567C DATA IMACH(15) / -8099 / 3568C DATA IMACH(16) / 8190 / 3569C 3570C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 3571C 3572C DATA IMACH( 1) / 11 / 3573C DATA IMACH( 2) / 12 / 3574C DATA IMACH( 3) / 8 / 3575C DATA IMACH( 4) / 10 / 3576C DATA IMACH( 5) / 16 / 3577C DATA IMACH( 6) / 2 / 3578C DATA IMACH( 7) / 2 / 3579C DATA IMACH( 8) / 15 / 3580C DATA IMACH( 9) / 32767 / 3581C DATA IMACH(10) / 16 / 3582C DATA IMACH(11) / 6 / 3583C DATA IMACH(12) / -64 / 3584C DATA IMACH(13) / 63 / 3585C DATA IMACH(14) / 14 / 3586C DATA IMACH(15) / -64 / 3587C DATA IMACH(16) / 63 / 3588C 3589C MACHINE CONSTANTS FOR THE DEC ALPHA 3590C USING G_FLOAT 3591C 3592C DATA IMACH( 1) / 5 / 3593C DATA IMACH( 2) / 6 / 3594C DATA IMACH( 3) / 5 / 3595C DATA IMACH( 4) / 6 / 3596C DATA IMACH( 5) / 32 / 3597C DATA IMACH( 6) / 4 / 3598C DATA IMACH( 7) / 2 / 3599C DATA IMACH( 8) / 31 / 3600C DATA IMACH( 9) / 2147483647 / 3601C DATA IMACH(10) / 2 / 3602C DATA IMACH(11) / 24 / 3603C DATA IMACH(12) / -127 / 3604C DATA IMACH(13) / 127 / 3605C DATA IMACH(14) / 53 / 3606C DATA IMACH(15) / -1023 / 3607C DATA IMACH(16) / 1023 / 3608C 3609C MACHINE CONSTANTS FOR THE DEC ALPHA 3610C USING IEEE_FLOAT 3611C 3612C DATA IMACH( 1) / 5 / 3613C DATA IMACH( 2) / 6 / 3614C DATA IMACH( 3) / 6 / 3615C DATA IMACH( 4) / 6 / 3616C DATA IMACH( 5) / 32 / 3617C DATA IMACH( 6) / 4 / 3618C DATA IMACH( 7) / 2 / 3619C DATA IMACH( 8) / 31 / 3620C DATA IMACH( 9) / 2147483647 / 3621C DATA IMACH(10) / 2 / 3622C DATA IMACH(11) / 24 / 3623C DATA IMACH(12) / -125 / 3624C DATA IMACH(13) / 128 / 3625C DATA IMACH(14) / 53 / 3626C DATA IMACH(15) / -1021 / 3627C DATA IMACH(16) / 1024 / 3628C 3629C MACHINE CONSTANTS FOR THE DEC RISC 3630C 3631C DATA IMACH( 1) / 5 / 3632C DATA IMACH( 2) / 6 / 3633C DATA IMACH( 3) / 6 / 3634C DATA IMACH( 4) / 6 / 3635C DATA IMACH( 5) / 32 / 3636C DATA IMACH( 6) / 4 / 3637C DATA IMACH( 7) / 2 / 3638C DATA IMACH( 8) / 31 / 3639C DATA IMACH( 9) / 2147483647 / 3640C DATA IMACH(10) / 2 / 3641C DATA IMACH(11) / 24 / 3642C DATA IMACH(12) / -125 / 3643C DATA IMACH(13) / 128 / 3644C DATA IMACH(14) / 53 / 3645C DATA IMACH(15) / -1021 / 3646C DATA IMACH(16) / 1024 / 3647C 3648C MACHINE CONSTANTS FOR THE DEC VAX 3649C USING D_FLOATING 3650C 3651C DATA IMACH( 1) / 5 / 3652C DATA IMACH( 2) / 6 / 3653C DATA IMACH( 3) / 5 / 3654C DATA IMACH( 4) / 6 / 3655C DATA IMACH( 5) / 32 / 3656C DATA IMACH( 6) / 4 / 3657C DATA IMACH( 7) / 2 / 3658C DATA IMACH( 8) / 31 / 3659C DATA IMACH( 9) / 2147483647 / 3660C DATA IMACH(10) / 2 / 3661C DATA IMACH(11) / 24 / 3662C DATA IMACH(12) / -127 / 3663C DATA IMACH(13) / 127 / 3664C DATA IMACH(14) / 56 / 3665C DATA IMACH(15) / -127 / 3666C DATA IMACH(16) / 127 / 3667C 3668C MACHINE CONSTANTS FOR THE DEC VAX 3669C USING G_FLOATING 3670C 3671C DATA IMACH( 1) / 5 / 3672C DATA IMACH( 2) / 6 / 3673C DATA IMACH( 3) / 5 / 3674C DATA IMACH( 4) / 6 / 3675C DATA IMACH( 5) / 32 / 3676C DATA IMACH( 6) / 4 / 3677C DATA IMACH( 7) / 2 / 3678C DATA IMACH( 8) / 31 / 3679C DATA IMACH( 9) / 2147483647 / 3680C DATA IMACH(10) / 2 / 3681C DATA IMACH(11) / 24 / 3682C DATA IMACH(12) / -127 / 3683C DATA IMACH(13) / 127 / 3684C DATA IMACH(14) / 53 / 3685C DATA IMACH(15) / -1023 / 3686C DATA IMACH(16) / 1023 / 3687C 3688C MACHINE CONSTANTS FOR THE ELXSI 6400 3689C 3690C DATA IMACH( 1) / 5 / 3691C DATA IMACH( 2) / 6 / 3692C DATA IMACH( 3) / 6 / 3693C DATA IMACH( 4) / 6 / 3694C DATA IMACH( 5) / 32 / 3695C DATA IMACH( 6) / 4 / 3696C DATA IMACH( 7) / 2 / 3697C DATA IMACH( 8) / 32 / 3698C DATA IMACH( 9) / 2147483647 / 3699C DATA IMACH(10) / 2 / 3700C DATA IMACH(11) / 24 / 3701C DATA IMACH(12) / -126 / 3702C DATA IMACH(13) / 127 / 3703C DATA IMACH(14) / 53 / 3704C DATA IMACH(15) / -1022 / 3705C DATA IMACH(16) / 1023 / 3706C 3707C MACHINE CONSTANTS FOR THE HARRIS 220 3708C 3709C DATA IMACH( 1) / 5 / 3710C DATA IMACH( 2) / 6 / 3711C DATA IMACH( 3) / 0 / 3712C DATA IMACH( 4) / 6 / 3713C DATA IMACH( 5) / 24 / 3714C DATA IMACH( 6) / 3 / 3715C DATA IMACH( 7) / 2 / 3716C DATA IMACH( 8) / 23 / 3717C DATA IMACH( 9) / 8388607 / 3718C DATA IMACH(10) / 2 / 3719C DATA IMACH(11) / 23 / 3720C DATA IMACH(12) / -127 / 3721C DATA IMACH(13) / 127 / 3722C DATA IMACH(14) / 38 / 3723C DATA IMACH(15) / -127 / 3724C DATA IMACH(16) / 127 / 3725C 3726C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES 3727C 3728C DATA IMACH( 1) / 5 / 3729C DATA IMACH( 2) / 6 / 3730C DATA IMACH( 3) / 43 / 3731C DATA IMACH( 4) / 6 / 3732C DATA IMACH( 5) / 36 / 3733C DATA IMACH( 6) / 6 / 3734C DATA IMACH( 7) / 2 / 3735C DATA IMACH( 8) / 35 / 3736C DATA IMACH( 9) / O377777777777 / 3737C DATA IMACH(10) / 2 / 3738C DATA IMACH(11) / 27 / 3739C DATA IMACH(12) / -127 / 3740C DATA IMACH(13) / 127 / 3741C DATA IMACH(14) / 63 / 3742C DATA IMACH(15) / -127 / 3743C DATA IMACH(16) / 127 / 3744C 3745C MACHINE CONSTANTS FOR THE HP 730 3746C 3747C DATA IMACH( 1) / 5 / 3748C DATA IMACH( 2) / 6 / 3749C DATA IMACH( 3) / 6 / 3750C DATA IMACH( 4) / 6 / 3751C DATA IMACH( 5) / 32 / 3752C DATA IMACH( 6) / 4 / 3753C DATA IMACH( 7) / 2 / 3754C DATA IMACH( 8) / 31 / 3755C DATA IMACH( 9) / 2147483647 / 3756C DATA IMACH(10) / 2 / 3757C DATA IMACH(11) / 24 / 3758C DATA IMACH(12) / -125 / 3759C DATA IMACH(13) / 128 / 3760C DATA IMACH(14) / 53 / 3761C DATA IMACH(15) / -1021 / 3762C DATA IMACH(16) / 1024 / 3763C 3764C MACHINE CONSTANTS FOR THE HP 2100 3765C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 3766C 3767C DATA IMACH( 1) / 5 / 3768C DATA IMACH( 2) / 6 / 3769C DATA IMACH( 3) / 4 / 3770C DATA IMACH( 4) / 1 / 3771C DATA IMACH( 5) / 16 / 3772C DATA IMACH( 6) / 2 / 3773C DATA IMACH( 7) / 2 / 3774C DATA IMACH( 8) / 15 / 3775C DATA IMACH( 9) / 32767 / 3776C DATA IMACH(10) / 2 / 3777C DATA IMACH(11) / 23 / 3778C DATA IMACH(12) / -128 / 3779C DATA IMACH(13) / 127 / 3780C DATA IMACH(14) / 39 / 3781C DATA IMACH(15) / -128 / 3782C DATA IMACH(16) / 127 / 3783C 3784C MACHINE CONSTANTS FOR THE HP 2100 3785C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 3786C 3787C DATA IMACH( 1) / 5 / 3788C DATA IMACH( 2) / 6 / 3789C DATA IMACH( 3) / 4 / 3790C DATA IMACH( 4) / 1 / 3791C DATA IMACH( 5) / 16 / 3792C DATA IMACH( 6) / 2 / 3793C DATA IMACH( 7) / 2 / 3794C DATA IMACH( 8) / 15 / 3795C DATA IMACH( 9) / 32767 / 3796C DATA IMACH(10) / 2 / 3797C DATA IMACH(11) / 23 / 3798C DATA IMACH(12) / -128 / 3799C DATA IMACH(13) / 127 / 3800C DATA IMACH(14) / 55 / 3801C DATA IMACH(15) / -128 / 3802C DATA IMACH(16) / 127 / 3803C 3804C MACHINE CONSTANTS FOR THE HP 9000 3805C 3806C DATA IMACH( 1) / 5 / 3807C DATA IMACH( 2) / 6 / 3808C DATA IMACH( 3) / 6 / 3809C DATA IMACH( 4) / 7 / 3810C DATA IMACH( 5) / 32 / 3811C DATA IMACH( 6) / 4 / 3812C DATA IMACH( 7) / 2 / 3813C DATA IMACH( 8) / 32 / 3814C DATA IMACH( 9) / 2147483647 / 3815C DATA IMACH(10) / 2 / 3816C DATA IMACH(11) / 24 / 3817C DATA IMACH(12) / -126 / 3818C DATA IMACH(13) / 127 / 3819C DATA IMACH(14) / 53 / 3820C DATA IMACH(15) / -1015 / 3821C DATA IMACH(16) / 1017 / 3822C 3823C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, 3824C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND 3825C THE PERKIN ELMER (INTERDATA) 7/32. 3826C 3827C DATA IMACH( 1) / 5 / 3828C DATA IMACH( 2) / 6 / 3829C DATA IMACH( 3) / 7 / 3830C DATA IMACH( 4) / 6 / 3831C DATA IMACH( 5) / 32 / 3832C DATA IMACH( 6) / 4 / 3833C DATA IMACH( 7) / 2 / 3834C DATA IMACH( 8) / 31 / 3835C DATA IMACH( 9) / Z7FFFFFFF / 3836C DATA IMACH(10) / 16 / 3837C DATA IMACH(11) / 6 / 3838C DATA IMACH(12) / -64 / 3839C DATA IMACH(13) / 63 / 3840C DATA IMACH(14) / 14 / 3841C DATA IMACH(15) / -64 / 3842C DATA IMACH(16) / 63 / 3843C 3844C MACHINE CONSTANTS FOR THE IBM PC 3845C 3846C DATA IMACH( 1) / 5 / 3847C DATA IMACH( 2) / 6 / 3848C DATA IMACH( 3) / 0 / 3849C DATA IMACH( 4) / 0 / 3850C DATA IMACH( 5) / 32 / 3851C DATA IMACH( 6) / 4 / 3852C DATA IMACH( 7) / 2 / 3853C DATA IMACH( 8) / 31 / 3854C DATA IMACH( 9) / 2147483647 / 3855C DATA IMACH(10) / 2 / 3856C DATA IMACH(11) / 24 / 3857C DATA IMACH(12) / -125 / 3858C DATA IMACH(13) / 127 / 3859C DATA IMACH(14) / 53 / 3860C DATA IMACH(15) / -1021 / 3861C DATA IMACH(16) / 1023 / 3862C 3863C MACHINE CONSTANTS FOR THE IBM RS 6000 3864C 3865C DATA IMACH( 1) / 5 / 3866C DATA IMACH( 2) / 6 / 3867C DATA IMACH( 3) / 6 / 3868C DATA IMACH( 4) / 0 / 3869C DATA IMACH( 5) / 32 / 3870C DATA IMACH( 6) / 4 / 3871C DATA IMACH( 7) / 2 / 3872C DATA IMACH( 8) / 31 / 3873C DATA IMACH( 9) / 2147483647 / 3874C DATA IMACH(10) / 2 / 3875C DATA IMACH(11) / 24 / 3876C DATA IMACH(12) / -125 / 3877C DATA IMACH(13) / 128 / 3878C DATA IMACH(14) / 53 / 3879C DATA IMACH(15) / -1021 / 3880C DATA IMACH(16) / 1024 / 3881C 3882C MACHINE CONSTANTS FOR THE INTEL i860 3883C 3884C DATA IMACH( 1) / 5 / 3885C DATA IMACH( 2) / 6 / 3886C DATA IMACH( 3) / 6 / 3887C DATA IMACH( 4) / 6 / 3888C DATA IMACH( 5) / 32 / 3889C DATA IMACH( 6) / 4 / 3890C DATA IMACH( 7) / 2 / 3891C DATA IMACH( 8) / 31 / 3892C DATA IMACH( 9) / 2147483647 / 3893C DATA IMACH(10) / 2 / 3894C DATA IMACH(11) / 24 / 3895C DATA IMACH(12) / -125 / 3896C DATA IMACH(13) / 128 / 3897C DATA IMACH(14) / 53 / 3898C DATA IMACH(15) / -1021 / 3899C DATA IMACH(16) / 1024 / 3900C 3901C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR) 3902C 3903C DATA IMACH( 1) / 5 / 3904C DATA IMACH( 2) / 6 / 3905C DATA IMACH( 3) / 5 / 3906C DATA IMACH( 4) / 6 / 3907C DATA IMACH( 5) / 36 / 3908C DATA IMACH( 6) / 5 / 3909C DATA IMACH( 7) / 2 / 3910C DATA IMACH( 8) / 35 / 3911C DATA IMACH( 9) / "377777777777 / 3912C DATA IMACH(10) / 2 / 3913C DATA IMACH(11) / 27 / 3914C DATA IMACH(12) / -128 / 3915C DATA IMACH(13) / 127 / 3916C DATA IMACH(14) / 54 / 3917C DATA IMACH(15) / -101 / 3918C DATA IMACH(16) / 127 / 3919C 3920C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR) 3921C 3922C DATA IMACH( 1) / 5 / 3923C DATA IMACH( 2) / 6 / 3924C DATA IMACH( 3) / 5 / 3925C DATA IMACH( 4) / 6 / 3926C DATA IMACH( 5) / 36 / 3927C DATA IMACH( 6) / 5 / 3928C DATA IMACH( 7) / 2 / 3929C DATA IMACH( 8) / 35 / 3930C DATA IMACH( 9) / "377777777777 / 3931C DATA IMACH(10) / 2 / 3932C DATA IMACH(11) / 27 / 3933C DATA IMACH(12) / -128 / 3934C DATA IMACH(13) / 127 / 3935C DATA IMACH(14) / 62 / 3936C DATA IMACH(15) / -128 / 3937C DATA IMACH(16) / 127 / 3938C 3939C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING 3940C 32-BIT INTEGER ARITHMETIC. 3941C 3942C DATA IMACH( 1) / 5 / 3943C DATA IMACH( 2) / 6 / 3944C DATA IMACH( 3) / 5 / 3945C DATA IMACH( 4) / 6 / 3946C DATA IMACH( 5) / 32 / 3947C DATA IMACH( 6) / 4 / 3948C DATA IMACH( 7) / 2 / 3949C DATA IMACH( 8) / 31 / 3950C DATA IMACH( 9) / 2147483647 / 3951C DATA IMACH(10) / 2 / 3952C DATA IMACH(11) / 24 / 3953C DATA IMACH(12) / -127 / 3954C DATA IMACH(13) / 127 / 3955C DATA IMACH(14) / 56 / 3956C DATA IMACH(15) / -127 / 3957C DATA IMACH(16) / 127 / 3958C 3959C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING 3960C 16-BIT INTEGER ARITHMETIC. 3961C 3962C DATA IMACH( 1) / 5 / 3963C DATA IMACH( 2) / 6 / 3964C DATA IMACH( 3) / 5 / 3965C DATA IMACH( 4) / 6 / 3966C DATA IMACH( 5) / 16 / 3967C DATA IMACH( 6) / 2 / 3968C DATA IMACH( 7) / 2 / 3969C DATA IMACH( 8) / 15 / 3970C DATA IMACH( 9) / 32767 / 3971C DATA IMACH(10) / 2 / 3972C DATA IMACH(11) / 24 / 3973C DATA IMACH(12) / -127 / 3974C DATA IMACH(13) / 127 / 3975C DATA IMACH(14) / 56 / 3976C DATA IMACH(15) / -127 / 3977C DATA IMACH(16) / 127 / 3978C 3979C MACHINE CONSTANTS FOR THE SILICON GRAPHICS 3980C 3981C DATA IMACH( 1) / 5 / 3982C DATA IMACH( 2) / 6 / 3983C DATA IMACH( 3) / 6 / 3984C DATA IMACH( 4) / 6 / 3985C DATA IMACH( 5) / 32 / 3986C DATA IMACH( 6) / 4 / 3987C DATA IMACH( 7) / 2 / 3988C DATA IMACH( 8) / 31 / 3989C DATA IMACH( 9) / 2147483647 / 3990C DATA IMACH(10) / 2 / 3991C DATA IMACH(11) / 24 / 3992C DATA IMACH(12) / -125 / 3993C DATA IMACH(13) / 128 / 3994C DATA IMACH(14) / 53 / 3995C DATA IMACH(15) / -1021 / 3996C DATA IMACH(16) / 1024 / 3997C 3998C MACHINE CONSTANTS FOR THE SUN 3999C 4000C DATA IMACH( 1) / 5 / 4001C DATA IMACH( 2) / 6 / 4002C DATA IMACH( 3) / 6 / 4003C DATA IMACH( 4) / 6 / 4004C DATA IMACH( 5) / 32 / 4005C DATA IMACH( 6) / 4 / 4006C DATA IMACH( 7) / 2 / 4007C DATA IMACH( 8) / 31 / 4008C DATA IMACH( 9) / 2147483647 / 4009C DATA IMACH(10) / 2 / 4010C DATA IMACH(11) / 24 / 4011C DATA IMACH(12) / -125 / 4012C DATA IMACH(13) / 128 / 4013C DATA IMACH(14) / 53 / 4014C DATA IMACH(15) / -1021 / 4015C DATA IMACH(16) / 1024 / 4016C 4017C MACHINE CONSTANTS FOR THE SUN 4018C USING THE -r8 COMPILER OPTION 4019C 4020C DATA IMACH( 1) / 5 / 4021C DATA IMACH( 2) / 6 / 4022C DATA IMACH( 3) / 6 / 4023C DATA IMACH( 4) / 6 / 4024C DATA IMACH( 5) / 32 / 4025C DATA IMACH( 6) / 4 / 4026C DATA IMACH( 7) / 2 / 4027C DATA IMACH( 8) / 31 / 4028C DATA IMACH( 9) / 2147483647 / 4029C DATA IMACH(10) / 2 / 4030C DATA IMACH(11) / 53 / 4031C DATA IMACH(12) / -1021 / 4032C DATA IMACH(13) / 1024 / 4033C DATA IMACH(14) / 113 / 4034C DATA IMACH(15) / -16381 / 4035C DATA IMACH(16) / 16384 / 4036C 4037C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER 4038C 4039C DATA IMACH( 1) / 5 / 4040C DATA IMACH( 2) / 6 / 4041C DATA IMACH( 3) / 1 / 4042C DATA IMACH( 4) / 6 / 4043C DATA IMACH( 5) / 36 / 4044C DATA IMACH( 6) / 4 / 4045C DATA IMACH( 7) / 2 / 4046C DATA IMACH( 8) / 35 / 4047C DATA IMACH( 9) / O377777777777 / 4048C DATA IMACH(10) / 2 / 4049C DATA IMACH(11) / 27 / 4050C DATA IMACH(12) / -128 / 4051C DATA IMACH(13) / 127 / 4052C DATA IMACH(14) / 60 / 4053C DATA IMACH(15) / -1024 / 4054C DATA IMACH(16) / 1023 / 4055C 4056C MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR 4057C 4058C DATA IMACH( 1) / 1 / 4059C DATA IMACH( 2) / 1 / 4060C DATA IMACH( 3) / 0 / 4061C DATA IMACH( 4) / 1 / 4062C DATA IMACH( 5) / 16 / 4063C DATA IMACH( 6) / 2 / 4064C DATA IMACH( 7) / 2 / 4065C DATA IMACH( 8) / 15 / 4066C DATA IMACH( 9) / 32767 / 4067C DATA IMACH(10) / 2 / 4068C DATA IMACH(11) / 24 / 4069C DATA IMACH(12) / -127 / 4070C DATA IMACH(13) / 127 / 4071C DATA IMACH(14) / 56 / 4072C DATA IMACH(15) / -127 / 4073C DATA IMACH(16) / 127 / 4074C 4075C***FIRST EXECUTABLE STATEMENT I1MACH 4076 IF (I .LT. 1 .OR. I .GT. 16) GO TO 10 4077C 4078 I1MACH = IMACH(I) 4079 RETURN 4080C 4081 10 CONTINUE 4082 WRITE (UNIT = OUTPUT, FMT = 9000) 4083 9000 FORMAT ('1ERROR 1 IN I1MACH - I OUT OF BOUNDS') 4084C 4085C CALL FDUMP 4086C 4087 STOP 4088 END 4089*DECK DHSTRT 4090 SUBROUTINE DHSTRT (DF, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL, 4091 + BIG, SPY, PV, YP, SF, RPAR, IPAR, H) 4092C***BEGIN PROLOGUE DHSTRT 4093C***SUBSIDIARY 4094C***PURPOSE Subsidiary to DDEABM, DDEBDF and DDERKF 4095C***LIBRARY SLATEC 4096C***TYPE DOUBLE PRECISION (HSTART-S, DHSTRT-D) 4097C***AUTHOR Watts, H. A., (SNLA) 4098C***DESCRIPTION 4099C 4100C DHSTRT computes a starting step size to be used in solving initial 4101C value problems in ordinary differential equations. 4102C 4103C ********************************************************************** 4104C ABSTRACT 4105C 4106C Subroutine DHSTRT computes a starting step size to be used by an 4107C initial value method in solving ordinary differential equations. 4108C It is based on an estimate of the local Lipschitz constant for the 4109C differential equation (lower bound on a norm of the Jacobian) , 4110C a bound on the differential equation (first derivative) , and 4111C a bound on the partial derivative of the equation with respect to 4112C the independent variable. 4113C (all approximated near the initial point A) 4114C 4115C Subroutine DHSTRT uses a function subprogram DHVNRM for computing 4116C a vector norm. The maximum norm is presently utilized though it 4117C can easily be replaced by any other vector norm. It is presumed 4118C that any replacement norm routine would be carefully coded to 4119C prevent unnecessary underflows or overflows from occurring, and 4120C also, would not alter the vector or number of components. 4121C 4122C ********************************************************************** 4123C On input you must provide the following 4124C 4125C DF -- This is a subroutine of the form 4126C DF(X,U,UPRIME,RPAR,IPAR) 4127C which defines the system of first order differential 4128C equations to be solved. For the given values of X and the 4129C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must 4130C evaluate the NEQ components of the system of differential 4131C equations DU/DX=DF(X,U) and store the derivatives in the 4132C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for 4133C equations I=1,...,NEQ. 4134C 4135C Subroutine DF must not alter X or U(*). You must declare 4136C the name DF in an external statement in your program that 4137C calls DHSTRT. You must dimension U and UPRIME in DF. 4138C 4139C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter 4140C arrays which you can use for communication between your 4141C program and subroutine DF. They are not used or altered by 4142C DHSTRT. If you do not need RPAR or IPAR, ignore these 4143C parameters by treating them as dummy arguments. If you do 4144C choose to use them, dimension them in your program and in 4145C DF as arrays of appropriate length. 4146C 4147C NEQ -- This is the number of (first order) differential equations 4148C to be integrated. 4149C 4150C A -- This is the initial point of integration. 4151C 4152C B -- This is a value of the independent variable used to define 4153C the direction of integration. A reasonable choice is to 4154C set B to the first point at which a solution is desired. 4155C You can also use B, if necessary, to restrict the length 4156C of the first integration step because the algorithm will 4157C not compute a starting step length which is bigger than 4158C ABS(B-A), unless B has been chosen too close to A. 4159C (it is presumed that DHSTRT has been called with B 4160C different from A on the machine being used. Also see the 4161C discussion about the parameter SMALL.) 4162C 4163C Y(*) -- This is the vector of initial values of the NEQ solution 4164C components at the initial point A. 4165C 4166C YPRIME(*) -- This is the vector of derivatives of the NEQ 4167C solution components at the initial point A. 4168C (defined by the differential equations in subroutine DF) 4169C 4170C ETOL -- This is the vector of error tolerances corresponding to 4171C the NEQ solution components. It is assumed that all 4172C elements are positive. Following the first integration 4173C step, the tolerances are expected to be used by the 4174C integrator in an error test which roughly requires that 4175C ABS(LOCAL ERROR) .LE. ETOL 4176C for each vector component. 4177C 4178C MORDER -- This is the order of the formula which will be used by 4179C the initial value method for taking the first integration 4180C step. 4181C 4182C SMALL -- This is a small positive machine dependent constant 4183C which is used for protecting against computations with 4184C numbers which are too small relative to the precision of 4185C floating point arithmetic. SMALL should be set to 4186C (approximately) the smallest positive DOUBLE PRECISION 4187C number such that (1.+SMALL) .GT. 1. on the machine being 4188C used. The quantity SMALL**(3/8) is used in computing 4189C increments of variables for approximating derivatives by 4190C differences. Also the algorithm will not compute a 4191C starting step length which is smaller than 4192C 100*SMALL*ABS(A). 4193C 4194C BIG -- This is a large positive machine dependent constant which 4195C is used for preventing machine overflows. A reasonable 4196C choice is to set big to (approximately) the square root of 4197C the largest DOUBLE PRECISION number which can be held in 4198C the machine. 4199C 4200C SPY(*),PV(*),YP(*),SF(*) -- These are DOUBLE PRECISION work 4201C arrays of length NEQ which provide the routine with needed 4202C storage space. 4203C 4204C RPAR,IPAR -- These are parameter arrays, of DOUBLE PRECISION and 4205C INTEGER type, respectively, which can be used for 4206C communication between your program and the DF subroutine. 4207C They are not used or altered by DHSTRT. 4208C 4209C ********************************************************************** 4210C On Output (after the return from DHSTRT), 4211C 4212C H -- is an appropriate starting step size to be attempted by the 4213C differential equation method. 4214C 4215C All parameters in the call list remain unchanged except for 4216C the working arrays SPY(*),PV(*),YP(*), and SF(*). 4217C 4218C ********************************************************************** 4219C 4220C***SEE ALSO DDEABM, DDEBDF, DDERKF 4221C***ROUTINES CALLED DHVNRM 4222C***REVISION HISTORY (YYMMDD) 4223C 820301 DATE WRITTEN 4224C 890531 Changed all specific intrinsics to generic. (WRB) 4225C 890831 Modified array declarations. (WRB) 4226C 890911 Removed unnecessary intrinsics. (WRB) 4227C 891024 Changed references from DVNORM to DHVNRM. (WRB) 4228C 891214 Prologue converted to Version 4.0 format. (BAB) 4229C 900328 Added TYPE section. (WRB) 4230C 910722 Updated AUTHOR section. (ALS) 4231C***END PROLOGUE DHSTRT 4232C 4233 INTEGER IPAR, J, K, LK, MORDER, NEQ 4234 DOUBLE PRECISION A, ABSDX, B, BIG, DA, DELF, DELY, 4235 1 DFDUB, DFDXB, DHVNRM, 4236 2 DX, DY, ETOL, FBND, H, PV, RELPER, RPAR, SF, SMALL, SPY, 4237 3 SRYDPB, TOLEXP, TOLMIN, TOLP, TOLSUM, Y, YDPB, YP, YPRIME 4238 DIMENSION Y(*),YPRIME(*),ETOL(*),SPY(*),PV(*),YP(*), 4239 1 SF(*),RPAR(*),IPAR(*) 4240 EXTERNAL DF 4241C 4242C .................................................................. 4243C 4244C BEGIN BLOCK PERMITTING ...EXITS TO 160 4245C***FIRST EXECUTABLE STATEMENT DHSTRT 4246 DX = B - A 4247 ABSDX = ABS(DX) 4248 RELPER = SMALL**0.375D0 4249C 4250C ............................................................... 4251C 4252C COMPUTE AN APPROXIMATE BOUND (DFDXB) ON THE PARTIAL 4253C DERIVATIVE OF THE EQUATION WITH RESPECT TO THE 4254C INDEPENDENT VARIABLE. PROTECT AGAINST AN OVERFLOW. 4255C ALSO COMPUTE A BOUND (FBND) ON THE FIRST DERIVATIVE 4256C LOCALLY. 4257C 4258 DA = SIGN(MAX(MIN(RELPER*ABS(A),ABSDX), 4259 1 100.0D0*SMALL*ABS(A)),DX) 4260 IF (DA .EQ. 0.0D0) DA = RELPER*DX 4261 CALL DF(A+DA,Y,SF,RPAR,IPAR) 4262 DO 10 J = 1, NEQ 4263 YP(J) = SF(J) - YPRIME(J) 4264 10 CONTINUE 4265 DELF = DHVNRM(YP,NEQ) 4266 DFDXB = BIG 4267 IF (DELF .LT. BIG*ABS(DA)) DFDXB = DELF/ABS(DA) 4268 FBND = DHVNRM(SF,NEQ) 4269C 4270C ............................................................... 4271C 4272C COMPUTE AN ESTIMATE (DFDUB) OF THE LOCAL LIPSCHITZ 4273C CONSTANT FOR THE SYSTEM OF DIFFERENTIAL EQUATIONS. THIS 4274C ALSO REPRESENTS AN ESTIMATE OF THE NORM OF THE JACOBIAN 4275C LOCALLY. THREE ITERATIONS (TWO WHEN NEQ=1) ARE USED TO 4276C ESTIMATE THE LIPSCHITZ CONSTANT BY NUMERICAL DIFFERENCES. 4277C THE FIRST PERTURBATION VECTOR IS BASED ON THE INITIAL 4278C DERIVATIVES AND DIRECTION OF INTEGRATION. THE SECOND 4279C PERTURBATION VECTOR IS FORMED USING ANOTHER EVALUATION OF 4280C THE DIFFERENTIAL EQUATION. THE THIRD PERTURBATION VECTOR 4281C IS FORMED USING PERTURBATIONS BASED ONLY ON THE INITIAL 4282C VALUES. COMPONENTS THAT ARE ZERO ARE ALWAYS CHANGED TO 4283C NON-ZERO VALUES (EXCEPT ON THE FIRST ITERATION). WHEN 4284C INFORMATION IS AVAILABLE, CARE IS TAKEN TO ENSURE THAT 4285C COMPONENTS OF THE PERTURBATION VECTOR HAVE SIGNS WHICH ARE 4286C CONSISTENT WITH THE SLOPES OF LOCAL SOLUTION CURVES. 4287C ALSO CHOOSE THE LARGEST BOUND (FBND) FOR THE FIRST 4288C DERIVATIVE. 4289C 4290C PERTURBATION VECTOR SIZE IS HELD 4291C CONSTANT FOR ALL ITERATIONS. COMPUTE 4292C THIS CHANGE FROM THE 4293C SIZE OF THE VECTOR OF INITIAL 4294C VALUES. 4295 DELY = RELPER*DHVNRM(Y,NEQ) 4296 IF (DELY .EQ. 0.0D0) DELY = RELPER 4297 DELY = SIGN(DELY,DX) 4298 DELF = DHVNRM(YPRIME,NEQ) 4299 FBND = MAX(FBND,DELF) 4300 IF (DELF .EQ. 0.0D0) GO TO 30 4301C USE INITIAL DERIVATIVES FOR FIRST PERTURBATION 4302 DO 20 J = 1, NEQ 4303 SPY(J) = YPRIME(J) 4304 YP(J) = YPRIME(J) 4305 20 CONTINUE 4306 GO TO 50 4307 30 CONTINUE 4308C CANNOT HAVE A NULL PERTURBATION VECTOR 4309 DO 40 J = 1, NEQ 4310 SPY(J) = 0.0D0 4311 YP(J) = 1.0D0 4312 40 CONTINUE 4313 DELF = DHVNRM(YP,NEQ) 4314 50 CONTINUE 4315C 4316 DFDUB = 0.0D0 4317 LK = MIN(NEQ+1,3) 4318 DO 140 K = 1, LK 4319C DEFINE PERTURBED VECTOR OF INITIAL VALUES 4320 DO 60 J = 1, NEQ 4321 PV(J) = Y(J) + DELY*(YP(J)/DELF) 4322 60 CONTINUE 4323 IF (K .EQ. 2) GO TO 80 4324C EVALUATE DERIVATIVES ASSOCIATED WITH PERTURBED 4325C VECTOR AND COMPUTE CORRESPONDING DIFFERENCES 4326 CALL DF(A,PV,YP,RPAR,IPAR) 4327 DO 70 J = 1, NEQ 4328 PV(J) = YP(J) - YPRIME(J) 4329 70 CONTINUE 4330 GO TO 100 4331 80 CONTINUE 4332C USE A SHIFTED VALUE OF THE INDEPENDENT VARIABLE 4333C IN COMPUTING ONE ESTIMATE 4334 CALL DF(A+DA,PV,YP,RPAR,IPAR) 4335 DO 90 J = 1, NEQ 4336 PV(J) = YP(J) - SF(J) 4337 90 CONTINUE 4338 100 CONTINUE 4339C CHOOSE LARGEST BOUNDS ON THE FIRST DERIVATIVE 4340C AND A LOCAL LIPSCHITZ CONSTANT 4341 FBND = MAX(FBND,DHVNRM(YP,NEQ)) 4342 DELF = DHVNRM(PV,NEQ) 4343C ...EXIT 4344 IF (DELF .GE. BIG*ABS(DELY)) GO TO 150 4345 DFDUB = MAX(DFDUB,DELF/ABS(DELY)) 4346C ......EXIT 4347 IF (K .EQ. LK) GO TO 160 4348C CHOOSE NEXT PERTURBATION VECTOR 4349 IF (DELF .EQ. 0.0D0) DELF = 1.0D0 4350 DO 130 J = 1, NEQ 4351 IF (K .EQ. 2) GO TO 110 4352 DY = ABS(PV(J)) 4353 IF (DY .EQ. 0.0D0) DY = DELF 4354 GO TO 120 4355 110 CONTINUE 4356 DY = Y(J) 4357 IF (DY .EQ. 0.0D0) DY = DELY/RELPER 4358 120 CONTINUE 4359 IF (SPY(J) .EQ. 0.0D0) SPY(J) = YP(J) 4360 IF (SPY(J) .NE. 0.0D0) DY = SIGN(DY,SPY(J)) 4361 YP(J) = DY 4362 130 CONTINUE 4363 DELF = DHVNRM(YP,NEQ) 4364 140 CONTINUE 4365 150 CONTINUE 4366C 4367C PROTECT AGAINST AN OVERFLOW 4368 DFDUB = BIG 4369 160 CONTINUE 4370C 4371C .................................................................. 4372C 4373C COMPUTE A BOUND (YDPB) ON THE NORM OF THE SECOND DERIVATIVE 4374C 4375 YDPB = DFDXB + DFDUB*FBND 4376C 4377C .................................................................. 4378C 4379C DEFINE THE TOLERANCE PARAMETER UPON WHICH THE STARTING STEP 4380C SIZE IS TO BE BASED. A VALUE IN THE MIDDLE OF THE ERROR 4381C TOLERANCE RANGE IS SELECTED. 4382C 4383 TOLMIN = BIG 4384 TOLSUM = 0.0D0 4385 DO 170 K = 1, NEQ 4386 TOLEXP = LOG10(ETOL(K)) 4387 TOLMIN = MIN(TOLMIN,TOLEXP) 4388 TOLSUM = TOLSUM + TOLEXP 4389 170 CONTINUE 4390 TOLP = 10.0D0**(0.5D0*(TOLSUM/NEQ + TOLMIN)/(MORDER+1)) 4391C 4392C .................................................................. 4393C 4394C COMPUTE A STARTING STEP SIZE BASED ON THE ABOVE FIRST AND 4395C SECOND DERIVATIVE INFORMATION 4396C 4397C RESTRICT THE STEP LENGTH TO BE NOT BIGGER 4398C THAN ABS(B-A). (UNLESS B IS TOO CLOSE 4399C TO A) 4400 H = ABSDX 4401C 4402 IF (YDPB .NE. 0.0D0 .OR. FBND .NE. 0.0D0) GO TO 180 4403C 4404C BOTH FIRST DERIVATIVE TERM (FBND) AND SECOND 4405C DERIVATIVE TERM (YDPB) ARE ZERO 4406 IF (TOLP .LT. 1.0D0) H = ABSDX*TOLP 4407 GO TO 200 4408 180 CONTINUE 4409C 4410 IF (YDPB .NE. 0.0D0) GO TO 190 4411C 4412C ONLY SECOND DERIVATIVE TERM (YDPB) IS ZERO 4413 IF (TOLP .LT. FBND*ABSDX) H = TOLP/FBND 4414 GO TO 200 4415 190 CONTINUE 4416C 4417C SECOND DERIVATIVE TERM (YDPB) IS NON-ZERO 4418 SRYDPB = SQRT(0.5D0*YDPB) 4419 IF (TOLP .LT. SRYDPB*ABSDX) H = TOLP/SRYDPB 4420 200 CONTINUE 4421C 4422C FURTHER RESTRICT THE STEP LENGTH TO BE NOT 4423C BIGGER THAN 1/DFDUB 4424 IF (H*DFDUB .GT. 1.0D0) H = 1.0D0/DFDUB 4425C 4426C FINALLY, RESTRICT THE STEP LENGTH TO BE NOT 4427C SMALLER THAN 100*SMALL*ABS(A). HOWEVER, IF 4428C A=0. AND THE COMPUTED H UNDERFLOWED TO ZERO, 4429C THE ALGORITHM RETURNS SMALL*ABS(B) FOR THE 4430C STEP LENGTH. 4431 H = MAX(H,100.0D0*SMALL*ABS(A)) 4432 IF (H .EQ. 0.0D0) H = SMALL*ABS(B) 4433C 4434C NOW SET DIRECTION OF INTEGRATION 4435 H = SIGN(H,DX) 4436C 4437 RETURN 4438 END 4439*DECK DHVNRM 4440 DOUBLE PRECISION FUNCTION DHVNRM (V, NCOMP) 4441C***BEGIN PROLOGUE DHVNRM 4442C***SUBSIDIARY 4443C***PURPOSE Subsidiary to DDEABM, DDEBDF and DDERKF 4444C***LIBRARY SLATEC 4445C***TYPE DOUBLE PRECISION (HVNRM-S, DHVNRM-D) 4446C***AUTHOR Watts, H. A., (SNLA) 4447C***DESCRIPTION 4448C 4449C Compute the maximum norm of the vector V(*) of length NCOMP and 4450C return the result as DHVNRM 4451C 4452C***SEE ALSO DDEABM, DDEBDF, DDERKF 4453C***ROUTINES CALLED (NONE) 4454C***REVISION HISTORY (YYMMDD) 4455C 820301 DATE WRITTEN 4456C 890531 Changed all specific intrinsics to generic. (WRB) 4457C 890831 Modified array declarations. (WRB) 4458C 891024 Changed references from DVNORM to DHVNRM. (WRB) 4459C 891024 Changed routine name from DVNORM to DHVNRM. (WRB) 4460C 891214 Prologue converted to Version 4.0 format. (BAB) 4461C 900328 Added TYPE section. (WRB) 4462C 910722 Updated AUTHOR section. (ALS) 4463C***END PROLOGUE DHVNRM 4464C 4465 INTEGER K, NCOMP 4466 DOUBLE PRECISION V 4467 DIMENSION V(*) 4468C***FIRST EXECUTABLE STATEMENT DHVNRM 4469 DHVNRM = 0.0D0 4470 DO 10 K = 1, NCOMP 4471 DHVNRM = MAX(DHVNRM,ABS(V(K))) 4472 10 CONTINUE 4473 RETURN 4474 END 4475*DECK J4SAVE 4476 FUNCTION J4SAVE (IWHICH, IVALUE, ISET) 4477C***BEGIN PROLOGUE J4SAVE 4478C***SUBSIDIARY 4479C***PURPOSE Save or recall global variables needed by error 4480C handling routines. 4481C***LIBRARY SLATEC (XERROR) 4482C***TYPE INTEGER (J4SAVE-I) 4483C***KEYWORDS ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR 4484C***AUTHOR Jones, R. E., (SNLA) 4485C***DESCRIPTION 4486C 4487C Abstract 4488C J4SAVE saves and recalls several global variables needed 4489C by the library error handling routines. 4490C 4491C Description of Parameters 4492C --Input-- 4493C IWHICH - Index of item desired. 4494C = 1 Refers to current error number. 4495C = 2 Refers to current error control flag. 4496C = 3 Refers to current unit number to which error 4497C messages are to be sent. (0 means use standard.) 4498C = 4 Refers to the maximum number of times any 4499C message is to be printed (as set by XERMAX). 4500C = 5 Refers to the total number of units to which 4501C each error message is to be written. 4502C = 6 Refers to the 2nd unit for error messages 4503C = 7 Refers to the 3rd unit for error messages 4504C = 8 Refers to the 4th unit for error messages 4505C = 9 Refers to the 5th unit for error messages 4506C IVALUE - The value to be set for the IWHICH-th parameter, 4507C if ISET is .TRUE. . 4508C ISET - If ISET=.TRUE., the IWHICH-th parameter will BE 4509C given the value, IVALUE. If ISET=.FALSE., the 4510C IWHICH-th parameter will be unchanged, and IVALUE 4511C is a dummy parameter. 4512C --Output-- 4513C The (old) value of the IWHICH-th parameter will be returned 4514C in the function value, J4SAVE. 4515C 4516C***SEE ALSO XERMSG 4517C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 4518C Error-handling Package, SAND82-0800, Sandia 4519C Laboratories, 1982. 4520C***ROUTINES CALLED (NONE) 4521C***REVISION HISTORY (YYMMDD) 4522C 790801 DATE WRITTEN 4523C 891214 Prologue converted to Version 4.0 format. (BAB) 4524C 900205 Minor modifications to prologue. (WRB) 4525C 900402 Added TYPE section. (WRB) 4526C 910411 Added KEYWORDS section. (WRB) 4527C 920501 Reformatted the REFERENCES section. (WRB) 4528C***END PROLOGUE J4SAVE 4529 LOGICAL ISET 4530 INTEGER IPARAM(9) 4531 SAVE IPARAM 4532 DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/ 4533 DATA IPARAM(5)/1/ 4534 DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/ 4535C***FIRST EXECUTABLE STATEMENT J4SAVE 4536 J4SAVE = IPARAM(IWHICH) 4537 IF (ISET) IPARAM(IWHICH) = IVALUE 4538 RETURN 4539 END 4540*DECK XERCNT 4541 SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL) 4542C***BEGIN PROLOGUE XERCNT 4543C***SUBSIDIARY 4544C***PURPOSE Allow user control over handling of errors. 4545C***LIBRARY SLATEC (XERROR) 4546C***CATEGORY R3C 4547C***TYPE ALL (XERCNT-A) 4548C***KEYWORDS ERROR, XERROR 4549C***AUTHOR Jones, R. E., (SNLA) 4550C***DESCRIPTION 4551C 4552C Abstract 4553C Allows user control over handling of individual errors. 4554C Just after each message is recorded, but before it is 4555C processed any further (i.e., before it is printed or 4556C a decision to abort is made), a call is made to XERCNT. 4557C If the user has provided his own version of XERCNT, he 4558C can then override the value of KONTROL used in processing 4559C this message by redefining its value. 4560C KONTRL may be set to any value from -2 to 2. 4561C The meanings for KONTRL are the same as in XSETF, except 4562C that the value of KONTRL changes only for this message. 4563C If KONTRL is set to a value outside the range from -2 to 2, 4564C it will be moved back into that range. 4565C 4566C Description of Parameters 4567C 4568C --Input-- 4569C LIBRAR - the library that the routine is in. 4570C SUBROU - the subroutine that XERMSG is being called from 4571C MESSG - the first 20 characters of the error message. 4572C NERR - same as in the call to XERMSG. 4573C LEVEL - same as in the call to XERMSG. 4574C KONTRL - the current value of the control flag as set 4575C by a call to XSETF. 4576C 4577C --Output-- 4578C KONTRL - the new value of KONTRL. If KONTRL is not 4579C defined, it will remain at its original value. 4580C This changed value of control affects only 4581C the current occurrence of the current message. 4582C 4583C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 4584C Error-handling Package, SAND82-0800, Sandia 4585C Laboratories, 1982. 4586C***ROUTINES CALLED (NONE) 4587C***REVISION HISTORY (YYMMDD) 4588C 790801 DATE WRITTEN 4589C 861211 REVISION DATE from Version 3.2 4590C 891214 Prologue converted to Version 4.0 format. (BAB) 4591C 900206 Routine changed from user-callable to subsidiary. (WRB) 4592C 900510 Changed calling sequence to include LIBRARY and SUBROUTINE 4593C names, changed routine name from XERCTL to XERCNT. (RWC) 4594C 920501 Reformatted the REFERENCES section. (WRB) 4595C***END PROLOGUE XERCNT 4596 CHARACTER*(*) LIBRAR, SUBROU, MESSG 4597C***FIRST EXECUTABLE STATEMENT XERCNT 4598 RETURN 4599 END 4600*DECK XERHLT 4601 SUBROUTINE XERHLT (MESSG) 4602C***BEGIN PROLOGUE XERHLT 4603C***SUBSIDIARY 4604C***PURPOSE Abort program execution and print error message. 4605C***LIBRARY SLATEC (XERROR) 4606C***CATEGORY R3C 4607C***TYPE ALL (XERHLT-A) 4608C***KEYWORDS ABORT PROGRAM EXECUTION, ERROR, XERROR 4609C***AUTHOR Jones, R. E., (SNLA) 4610C***DESCRIPTION 4611C 4612C Abstract 4613C ***Note*** machine dependent routine 4614C XERHLT aborts the execution of the program. 4615C The error message causing the abort is given in the calling 4616C sequence, in case one needs it for printing on a dayfile, 4617C for example. 4618C 4619C Description of Parameters 4620C MESSG is as in XERMSG. 4621C 4622C***REFERENCES R. E. Jones and D. K. Kahaner, XERROR, the SLATEC 4623C Error-handling Package, SAND82-0800, Sandia 4624C Laboratories, 1982. 4625C***ROUTINES CALLED (NONE) 4626C***REVISION HISTORY (YYMMDD) 4627C 790801 DATE WRITTEN 4628C 861211 REVISION DATE from Version 3.2 4629C 891214 Prologue converted to Version 4.0 format. (BAB) 4630C 900206 Routine changed from user-callable to subsidiary. (WRB) 4631C 900510 Changed calling sequence to delete length of character 4632C and changed routine name from XERABT to XERHLT. (RWC) 4633C 920501 Reformatted the REFERENCES section. (WRB) 4634C***END PROLOGUE XERHLT 4635 CHARACTER*(*) MESSG 4636C***FIRST EXECUTABLE STATEMENT XERHLT 4637 STOP 4638 END 4639