c 01/04/2000 c----------------------------------------------------------------------------- c----------------------------------------------------------------------------- c c MEBDF VERSION THAT USES MA28 (HARWELL SPARSE MATRIX ROUTINES) c c----------------------------------------------------------------------------- c----------------------------------------------------------------------------- SUBROUTINE MEBDFSD(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,LWORK, + WORK,LIWORK,IWORK,MAXDER,ITOL,RTOL,ATOL,F,PDERV,MAS, + U,NZ,LICN,LIRN,NZMAX,NZMAS,IPAR,RPAR,IERR) c *********************************************************************** c *********************************************************************** c Written by T.J. Abdulla, J.R. Cash and M.T. Diamantakis, c Department of Mathematics, c Imperial College, c London SW7 2AZ c England c t.abdulla@ic.ac.uk or j.cash@ic.ac.uk c The authors would be pleased to receive any comments c good or bad!! c c *********************************************************************** c *********************************************************************** c IMPLICIT DOUBLE PRECISION(A-H,O-Z) C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C THIS SUBROUTINE IS FOR THE PURPOSE * C OF SPLITTING UP THE WORK ARRAYS WORK AND IWORK * C FOR USE INSIDE THE INTEGRATOR STIFF * C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C .. SCALAR ARGUMENTS .. INTEGER IDID,LIWORK,LOUT,LWORK,MF,N,MAXDER,ITOL C .. C .. ARRAY ARGUMENTS .. DIMENSION WORK(LWORK),Y0(N),RTOL(1),ATOL(1),IWORK(LIWORK), & IPAR(*),RPAR(*) C .. C .. LOCAL SCALARS .. INTEGER I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14, + I15,I16,I17,I18,JK1,JK2,JK3,JK4,JK5,JK6,JK7,JK8,JK9, + JK10,JK11,JK12 C LOGICAL MA28B C COMMON BLOCKS COMMON /MA28FL/MA28B C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL OVDRIV,F,PDERV,MAS,SPARSDET C .. C .. SAVE STATEMENT .. SAVE I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14, + I15,I16,I17,I18,JK1,JK2,JK3,JK4,JK5,JK6,J7,JK8, + JK9,JK10,JK11,JK12,NZPW C IF (IDID.EQ.1) THEN IF (N.LE.0) THEN WRITE (LOUT,9020) N IDID = -4 ELSE I1 = N*12+ 3 I2 = I1 + N*12 I3 = I2 + N*2 I4 = I3 + N I5 = I4 + N I6 = I5 + N I7 = I6 + N I8 = I7 + N I9 = I8 + N I10 = I9 + N I11 = I10 + LICN I12 = I11 + LICN I13 = i12 + NZMAS I14 = I13 + N I15 = I14 + N I16 = I15 + 2*N I17 = I16 + N I18 = I17 + N C C.....INTEGER WORKSPACE C K = NZMAX + NZMAS JK1 = LICN + 15 JK2 = JK1 + LIRN JK3 = JK2 + 5*N JK4 = JK3 + 8*N JK5 = JK4 + NZMAS JK6 = JK5 + NZMAS JK7 = JK6 + N+1 JK8 = JK7 + N+1 JK9 = JK8 + K JK10 = JK9 + K JK11 = JK10 + 2*N+1 JK12 = JK11 + NZMAX UROUND=DLAMCH('Epsilon') WORK(1)=UROUND EPSJAC=SQRT(WORK(1)) C IF (LWORK.LT.(I18-1)) THEN IDID = -11 WRITE (LOUT,9000) I18 - 1 ENDIF IF (LIWORK.LT.JK12-1) THEN IDID = -12 WRITE (LOUT,9010) JK12-1 ENDIF ENDIF IF (IDID.LT.0) RETURN MA28B=.FALSE. ENDIF C C THE DIMENSION OF THE REAL WORKSPACE, WORK, HAS TO BE AT LEAST C (39*N+2*LICN+NZMAS+3). WHILE THE DIMENSION OF THE INTEGER c WORKSPACE HAS TO BE AT LEAST (17*N+LICN+LIRN+4*NZMAS+3*NZMAX+18), C CALL OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,WORK(3),WORK(I1), + WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7), + WORK(I8),WORK(I9),WORK(I10),WORK(I11),IWORK(15),IWORK(JK1), + IWORK(JK2),IWORK(JK3),MAXDER,ITOL,RTOL,ATOL,F,PDERV,MAS, + IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5),IWORK(6), + IWORK(7),IWORK(8),IWORK(9),IWORK(10),IWORK(11),IWORK(12), + IWORK(13),IWORK(14),UROUND,WORK(2),EPSJAC,U,NZ,NZPW,LICN, + LIRN,WORK(I12),IWORK(JK4),IWORK(JK5),NZMAS,IPAR,RPAR,IERR, + IWORK(JK6),WORK(I13),WORK(I14),WORK(I15),WORK(I16), + IWORK(JK7),NZMAX,IWORK(JK8),IWORK(JK9),IWORK(JK10), + WORK(I17),IWORK(JK11)) C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C WORK() HOUSES THE FOLLOWING ARRAYS C Y(N,12) , YHOLD(N,12) , YNHOLD(N,2) , YMAX(N) C ERRORS(N) , SAVE1(N) ,SAVE2(N),SCALE(N),ARH(N),PW(NZ+NZMAS) C PWCOPY(NZ+NZMAS) C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< RETURN 9000 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN DRIVER',/,/, + ' >>> REAL WORKSPACE IS INSUFFICIENT <<< ',/, + ' WORKSPACE MUST BE AT LEAST ',I8,' ELEMENTS LONG') 9010 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN DRIVER',/,/, + ' >>> INTEGER WORKSPACE IS INSUFFICIENT <<< ',/, + ' WORKSPACE MUST BE AT LEAST ',I8,' ELEMENTS LONG') 9020 FORMAT (/,/,' ***** ERROR ***** INTEGRATION HALTED IN DRIVER',/,/, + ' >>> ILLEGAL VALUE FOR NUMBER OF EQUATIONS <<< ',/, + ' WITH N = ',I8) END c-------------------------------------------------------------------- SUBROUTINE OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,Y,YHOLD, + YNHOLD,YMAX,ERRORS,SAVE1,SAVE2,SCALE,ARH,W,PW,PWCOPY, + ICN,IRN,IKEEP,IW,MAXDER,ITOL,RTOL,ATOL,F,PDERV,MAS, + NIND1,NIND2,NIND3,NQUSED,NSTEP,NFAIL,NFE,NJE,NDEC, + NBSOL,NPSET,NCOSET,MAXORD,MAXSTP,UROUND,HUSED,EPSJAC, + U,NZ,NZPW,LICN,LIRN,AM,JAM,IAM,NZMAS,IPAR,RPAR,IERR, + IAMPT,HDUMMY,YDOT,W1,Z,IP,NZMAX,JCOL,IROW,IG,HMAX1,IRNT) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C THIS IS THE NOVEMBER 6th 1998 VERSION OF OVDRIV, A PACKAGE FOR C THE SOLUTION OF THE INITIAL VALUE PROBLEM FOR SYSTEMS OF C ORDINARY DIFFERENTIAL EQUATIONS WITH SPARSE JACOBIANS C DY/DT = F(Y,T), Y=(Y(1),Y(2),Y(3), . . . ,Y(N)) C SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE C C REFERENCES C C 1. J. R. CASH, THE INTEGRATION OF STIFF INITIAL VALUE PROBLEMS C IN O.D.E.S USING MODIFIED EXTENDED BACKWARD DIFFERENTIATION C FORMULAE, COMP. AND MATHS. WITH APPLICS., 9, 645-657, (1983). C 2. J.R. CASH AND S. CONSIDINE, AN MEBDF CODE FOR STIFF C INITIAL VALUE PROBLEMS, ACM TRANS MATH SOFTWARE, 142-158, C (1992). C 3. J.R. CASH, STABLE RECURSIONS WITH APPLICATIONS TO THE C NUMERICAL SOLUTION OF STIFF SYSTEMS, ACADEMIC PRESS,(1979). C 4. A.C. HINDMARSH, ODEPACK, A SYSTEMISED COLLECTION OF ODE C SOLVERS, in SCIENTIFIC COMPUTING, R.S. STEPLEMAN et. al. C (eds) North-Holland, AMSTERDAM, pp55-64 , (1983). C 5. E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL C EQUATIONS II, STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS, C SPRINGER 1996, page 267. C ---------------------------------------------------------------- C OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND C IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR STIFF. C C THE INPUT PARAMETERS ARE .. C N = THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE C (USED ONLY ON THE FIRST CALL) C HO = THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE C FIRST CALL) C Y0 = A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF Y C (USED FOR INPUT ONLY ON FIRST CALL) C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT C TEND = END OF THE RANGE OF INTEGRATION. C MF = THE METHOD FLAG. AT PRESENT MF=25 IS C ALLOWED. THESE ARE EXTENDED BACKWARD DIFFERENTIATION C FORMULAE USING THE CHORD METHOD WITH ANALYTIC JACOBIAN. C IDID = THE INTEGER USED ON INPUT TO INDICATE THE TYPE OF CALL. C 1 THIS IS THE FIRST CALL FOR THE PROBLEM. C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM C AND INTEGRATION IS TO CONTINUE. C -1 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, C AND THE USER HAS RESET N, RTOL, ATOL, AND/OR MF. C 2 SAME AS 0 EXCEPT THAT TOUT HAS TO BE HIT C EXACTLY (NO INTERPOLATION IS DONE). C ASSUMES TOUT .GE. THE CURRENT T. C 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING C PROGRAM AFTER ONE STEP. TOUT IS IGNORED, UNTIL THE C INTEGRATION REACHES TOUT OR BEYOND. IF IT PASSES TOUT C THE PROGRAM INTERPOLATES THE SOLUTION VALUES AND C RETURNS THE SOLUTION VALUE AT TOUT. C SINCE THE NORMAL OUTPUT VALUE OF IDID IS 0, C IT NEED NOT BE RESET FOR NORMAL CONTINUATION. C LOUT = THE LOGICAL OUTPUT CHANNEL FOR MESSAGE PASSING. C MAXDER= THE MAXIMUM ORDER IS MAXDER + 1. C THE VALUE OF MAXDER CANNOT EXCEED 7. THIS IS THE C VALUE RECOMMENDED UNLESS IT IS BELIEVED THAT THERE C ARE SEVERE STABILITY PROBLEMS IN WHICH CASE MAXDER=3 C OR 4 SHOULD BE TRIED. C ITOL = AN INDICATOR OF THE TYPE OF ERROR CONTROL. SEE C DESCRIPTION BELOW UNDER ATOL. C RTOL = A RELATIVE ERROR TOLERANCE PARAMETER. CAN BE EITHER A C SCALAR OR AN ARRAY OF LENGTH N. SEE DESCRIPTION C BELOW UNDER ATOL. C ATOL = THE ABSOLUTE ERROR BOUND. C THE INPUT PARAMETERS ITOL, RTOL AND ATOL DETERMINE C THE ERROR CONTROL PERFORMED BY THE SOLVER. THE C SOLVER WILL CONTROL THE VECTOR e = (e(i)) OF ESTIMATED C LOCAL ERRORS IN y ACCORDING TO AN INEQUALITY OF THE FORM C RMS-NORM OF (e(i)/ewt(i)) .LE. 1 C THE ROOT MEAN SQUARE NORM IS C RMS-NORM(V) = SQRT((SUM v(i)**2)/N). HERE C ewt = (ewt(i)) IS A VECTOR OF WEIGHTS WHICH MUST C ALWAYS BE POSITIVE, AND THE VALUES OF RTOL AND ATOL C SHOULD BE NON-NEGATIVE. IF ITOL = 1 THEN SINGLE STEP ERROR C ESTIMATES DIVIDED BY YMAX(I) WILL BE KEPT LESS THAN 1 C IN ROOT-MEAN-SQUARE NORM. THE VECTOR YMAX OF WEIGHTS IS C COMPUTED IN OVDRIV. INITIALLY YMAX(I) IS SET AS C THE MAXIMUM OF 1 AND ABS(Y(I)). THEREAFTER YMAX(I) IS C THE LARGEST VALUE OF ABS(Y(I)) SEEN SO FAR, OR THE C INITIAL VALUE YMAX(I) IF THAT IS LARGER. C IF ITOL = 1 THE USER NEEDS TO SET ATOL = RTOL = C THE PRECISION REQUIRED. THEN C ewt(i) = RTOL(1)*YMAX(i) C IF ITOL IS GREATER THAN 1 THEN C ewt(i) = rtol(i)*abs(y(i)) + atol(i) C THE FOLLOWING TABLE GIVES THE TYPES (SCALAR/ARRAY) C OF RTOL AND ATOL, AND THE CORRESPONDING FORM OF ewt(i) C ITOL RTOL ATOL ewt(i) C 2 SCALAR SCALAR rtol*abs(y(i)) + atol C 3 SCALAR ARRAY rtol*abs(y(i)) + atol(i) C 4 ARRAY SCALAR rtol(i)*abs(y(i))+ atol C 5 ARRAY ARRAY rtol(i)*abs(y(i))+ atol(i) C IF EITHER OF THESE PARAMETERS IS A SCALAR, IT NEED C NOT BE DIMENSIONED IN THE USER'S CALLING PROGRAM. C C NIND1 = THE NUMBER OF VARIABLES OF INDEX 1,2,3 RESPECTIVELY. C NIND2,NIND3 THESE ARE SET IN IWORK(1),(2),(3). C THE EQUATIONS MUST BE DEFINED SO THAT THE INDEX 1 C VARIABLES PRECEDE THE INDEX 2 VARIABLES WHICH IN C TURN PRECEDE THE INDEX 3 VARIABLES. C RPAR, IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) C WHICH CAN BE USED FOR COMMUNICATION BETWEEN THE USER'S C CALLING PROGRAM AND THE F, PDERV AND MAS SUBROUTINES. C IERR IERR IS AN INTEGER FLAG WHICH IS ALWAYS EQUAL TO ZERO C ON INPUT. SUBROUTINES F, PDERV AND MAS SHOULD ALTER C IERR ONLY IF ONE OF THEM ENCOUNTERS AN ILLEGAL OPERATION SUCH C AS THE SQUARE ROOT OF A NEGATIVE NUMBER OR EXPONENT C OVERFLOW. THE USER CAN THEN ALTER H AND CALL THE C SUBROUTINE AGAIN WITH IDID=-1 IF HE WISHES. C C C NZ THE NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX. C C LICN INTEGER VARIABLE. IT IS THE LENGHT OF ARRAY PW, ICN AND IRN C IT SHOULD BE BIG ENOGH TO HOLD ALL THE NONZEROS OF L AND U C AND LEAVE SOME ROOM. RECOMENDED VALUE BE BETWEEN 2*NZ AND C 3*NZ. C C U IS A REAL VARIABLE BETWEEN ZERO AND ONE TO DETERMINE THE C BALANCE BETWEEN SPARSITY AND STABILITY PIVOTING IN THE SPARSE C SOLVERS. A VALUES NEAR 0.0 EMPHASIZING SPARSITY AND VALUES NEAR C 1.0 EMPHASIZING STABILITY. C C IROW IS AN INTEGER ARRAY OF LENGTH NZ+NZMAS. THE FIRST NZ C ENTRIES HOLD THE ROW NUMBERS OF THE NONZERO JACOBIAN ELEMENTS. C THE USER SETS THEM. THE LAST N ENTRIES ARE PROVIDED FOR C STORING THE ELEMENTS OF THE MAS MATRIX M WHICH MAY HAVE C COORDINATES DIFFERENT THAN IROW(J),JCOL(J) J<=NZ. C C JCOL IS AN INTEGER ARRAY OF LENGTH NZ+NZMAS. THE FIRST NZ ENTRIES C HOLD THE COLUMN NUMBERS OF THE NONZERO JACOBIAN ENTRIES. C IROW,JCOL DEFINE THE SPARSE JACOBIAN STRUCTURE AND REMAIN C CONSTANT DURING THE INTEGRATION. C C NZ THE NUMBER OF NONZERO ELEMENTS IN THE MASS MATRIX. C C IAM IS AN INTEGER ARRAY OF LENGTH NZMAS HOLD THE ROW NUMBERS OF C THE NONZERO MASS ELEMENTS. C C JAM IS AN INTEGER ARRAY OF LENGTH NZMAS HOLD THE COLUMN NUMBERS OF C THE NONZERO MASS ELEMENTS. C C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURED AND A NORMAL C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. C A CHANGE OF PARAMETERS WITH IDID = -1 CAN BE MADE AFTER C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN. C C THE OUTPUT PARAMETERS ARE.. C T0 = THE VALUE OF T WHICH RELATES TO THE CURRENT SOLUTION C POINT Y0() C HO = THE STEPSIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT. C Y0 = THE COMPUTED VALUES OF Y AT T = TOUT C TOUT = UNCHANGED FROM ITS INPUT VALUE. C IDID = INTEGER USED ON OUTPUT TO INDICATE RESULTS, WITH C THE FOLLOWING VALUES AND MEANINGS.. C C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. C C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE C ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF C 1.E10 FROM ITS INITIAL VALUE. C C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY C A TEST ON RTOL/ATOL. TOO MUCH ACCURACY HAS BEEN REQUESTED. C C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE C CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A C FACTOR OF 1.E10 FROM ITS INITIAL VALUE. C C -4 IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT C PARAMETERS. SEE PRINTED MESSAGE. C C -5 IDID WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT C WAS NOT BEYOND T. INTERPOLATION AT T = TOUT WAS C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, C SIMPLY CALL AGAIN WITH IDID = -1 AND A NEW TOUT. C C C -6 MAXIMUM ALLOWABLE NUMBER OF INTEGRATION STEPS EXCEEDED. C TO CONTINUE THE USER SHOULD RESET IWORK(14). C C C -7 STEPSIZE IS TOO SMALL (LESS THAN SQRT(UROUND)/100) C C C -11 INSUFFICIENT REAL WORKSPACE FOR THE INTEGRATION C C -12 INSUFFICIENT INTEGER WORKSPACE FOR THE INTEGRATION C C IN ADDITION TO OVDRIVE, THE FOLLOWING ROUTINES ARE PROVIDED C IN THE PACKAGE.. C C INTERP( - ) INTERPOLATES TO GET THE OUTPUT VALUES C AT T=TOUT FROM THE DATA IN THE Y ARRAY. C STIFF( - ) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A C SINGLE STEP AND ASSOCIATED ERROR CONTROL. C COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION C SCHEMES FOR USE IN THE CORE INTEGRATOR. C PSET( - ) COMPUTES AND PROCESSES THE JACOBIAN C MATRIX J = DF/DY C HAS BEEN CALLED FOR THE MATRIX A C MA28AD( - ) PREORDERS THE MATRIX INTO BLOCK TRIANGULAR FROM C AND THEN PEROFORMS A SPARSITY LU FACTORIZATION. C MA28BD( - ) PEROFORMS A SPARSITY LU FACTORIZATION TO A MATRIX C PREVIOUSLY FACTORED BY MA28AD. C MA28CD( - ) SOLVES A SPARSE LINEAR SYSTEM A*x=b USING THE C FACTORIZATION OBTAINED BY MA28AD OR MA28BD. C AND HARWELL ROUTINES C MC13D, MC23E, MC20AD, MC21A,MC21B, MA22AD, C MC23DA, MC24AD, MCA30AD, MA30DD, MA30CD, MA30D C C THE FOLLOWING ROUTINES ARE TO BE SUPPLIED BY THE USER AND C SHOULD BE DECLARED AS EXTERNAL. C C F(N,T,Y,YDOT,IPAR,RPAR,IERR) C COMPUTES THE FUNCTION YDOT = F(Y,T), THE C RIGHT HAND SIDE OF THE O.D.E. C HERE Y AND YDOT ARE VECTORS OF LENGTH N C C PDERV(N,T,Y,NZ,PD,IROW,JCOL,IPAR,RPAR,IERR) C COMPUTES THE NONZERO ENTRIES OF THE JACOBIAN MATRIX C AND STORES IT IN PD. C IF THE JACOBIAN IS SPARSE THE ARRAYS C PD(I) (I<=NZ) CONTAINS THE JACOBIAN ELEMENT C AT LOCATION (IROW(I),JCOL(I)). C PDERV IS CALLED ONLY IF MITER = 5. C OTHERWISE A DUMMY ROUTINE CAN BE SUBSTITUTED. C C MAS(N,NZMAS,AM,IAM,JAM,IPAR,RPAR) C COMPUTES THE NONZER ENTRIES OF THE MASS MATRIX M C AND STORES THEM IN THE ARRAYS AM(I) AT LOCATION C (IAM(I),JAM(I)). C C C THE DIMENSIONS OF YMAX,ERROR,SAVE1,SAVE2,IPIV AND THE FIRST DIMENSION C OF Y SHOULD ALL BE AT LEAST N. C C C C UROUND THIS IS THE UNIT ROUNDOFF AND HAS TO BE SET AS C UROUND = DLAMCH('Epsilon') C EPSJAC = sqrt(UROUND). C C HUSED (=WORK(2)) LAST STEPSIZE SUCCESSFULLY USED BY THE INTEGRATOR C NQUSED (=IWORK(4)) LAST ORDER SUCCESSFULLY USED C NSTEP (=IWORK(5)) NUMBER OF SUCCESSFUL STEPS TAKEN SO FAR C NFAIL (=IWORK(6)) NUMBER OF FAILED STEPS C NFE (=IWORK(7)) NUMBER OF FUNCTION EVALUATIONS SO FAR C NJE (=IWORK(8)) NUMBER OF JACOBIAN EVALUATIONS SO FAR C NDEC (=IWORK(9)) NUMBER OF LU DECOMPOSITIONS SO FAR C NBSOL (=IWORK(10)) NUMBER OF 'BACKSOLVES' SO FAR C NPSET (=IWORK(11)) NUMBER OF TIMES A NEW COEFFICIENT MATRIX HAS BEEN C FORMED SO FAR C NCOSET (=IWORK(12)) NUMBER OF TIMES THE ORDER OF THE METHOD USED HAS C BEEN CHANGED SO FAR C MAXORD (=IWORK(13)) THE MAXIMUM ORDER USED SO FAR IN THE INTEGRATION c MAXSTP (=IWORK(14)) THE MAXIMUM ALLOWED NUMBER OF STEPS SET BY THE USER C IF IT IS ONLY REQUIRED TO CONTROL THE ACCURACY IN THE C DIFFERENTIAL VARIABLES THEN THE USER SHOULD FIND THE C STRING 'AMMEND' AND MAKE CHANGES THERE C C START OF PROGRAM PROPER C C .. SCALAR ARGUMENTS .. INTEGER IDID,LOUT,MF,N,MAXDER,ITOL,NZ C .. C .. ARRAY ARGUMENTS .. DIMENSION ARH(N),ERRORS(N),PW(*),PWCOPY(*),SAVE1(N), + SAVE2(N),Y(N,12),Y0(N),YHOLD(N,12),YMAX(N),YNHOLD(N,2), + RTOL(1),ATOL(1),SCALE(N),ICN(LICN),IRN(LIRN),IKEEP(N,5), + IW(N,8),W(N),AM(*),JAM(*),IAM(*),JCOL(*),IROW(*), + IPAR(*),RPAR(*),IAMPT(N+1),HDUMMY(N),YDOT(N), + IP(N+1),W1(2*N),Z(N),IG(2*N+1),HMAX1(N),IRNT(NZMAX) C .. C .. LOCAL SCALARS .. INTEGER I,KGO,NHCUT,LICN,LIRN,NZMAX C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INTERP,STIFF,F,PDERV,MAS C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1 C .. INTEGER JSTART,KFLAG,MAXORD, + NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP SAVE T,H,HMIN,HMAX,KFLAG,JSTART C .. IF (IDID.EQ.0) THEN C I.E. NORMAL CONTINUATION OF INTEGRATION T0=T HMAX = DABS(TEND-T0)*10.0D+0 IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT THE OUTPUT POINT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG T0 = TOUT HO = H RETURN ENDIF ELSE IF (IDID.EQ.2) THEN C I.E. CONTINUING INTEGRATION BUT WISH TO HIT TOUT T0 = T HMAX = DABS(TEND-T0)*10.0D+0 IF (((T+H)-TOUT)*H.GT.0.0D+0) THEN C WE HAVE ALREADY OVERSHOT THE OUTPUT POINT OR WE WILL C DO SO ON THE NEXT STEP IF (((T-TOUT)*H.GE.0.0D+0) .OR. (DABS(T-TOUT).LE. + 100.0D+0*UROUND*HMAX)) THEN C HAVE OVERSHOT THE OUTPUT POINT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT HO = H IDID = KFLAG RETURN ELSE C WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE C SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY' H = (TOUT-T)* (1.0D+0-4.0D+0*UROUND) JSTART = -1 ENDIF ENDIF ELSE IF (IDID.EQ.-1) THEN C NOT FIRST CALL BUT PARAMETERS RESET H=HO IF(H.LT.EPSJAC/100.0D+0) THEN WRITE(6,9160) IDID = -7 RETURN ENDIF T0 = T IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT TOUT WRITE (LOUT,9080) T,TOUT,H CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) HO = H T0 = TOUT IDID = -5 RETURN ELSE JSTART = -1 ENDIF ELSE IF (IDID.EQ.3) THEN T0 = T IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT,SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG T0 = TOUT HO = H RETURN ENDIF ELSE C IDID SHOULD BE 1 AND THIS IS THE FIRST CALL FOR THIS PROBLEM C CHECK THE ARGUMENTS THAT WERE PASSED FOR CORRECTNESS IF (IDID.NE.1) THEN C VALUE OF IDID NOT ALLOWED WRITE (LOUT,9070) IDID IDID = -4 ENDIF NN=N IF(ITOL.LE.3) NN = 1 DO 1 I=1,NN IF (RTOL(I).LT.0.0D+0) THEN C ILLEGAL VALUE FOR RELATIVE ERROR TOLERENCE WRITE (LOUT,9040) IDID = -4 ENDIF 1 CONTINUE NN=N IF(ITOL.EQ.1.OR.ITOL.EQ.2.OR.ITOL.EQ.4) NN=1 DO 2 I=1,NN IF (ATOL(I).LT.0.0D+0) THEN C ILLEGAL ABSOLUTE ERROR TOLERANCE WRITE(LOUT,9045) IDID=-4 ENDIF 2 CONTINUE IF(ITOL.EQ.1.AND.RTOL(1).EQ.0) THEN C ILLEGAL ERROR TOLERANCE WRITE(LOUT,9040) IDID = -4 ENDIF IF(ITOL.NE.1) THEN VHOLD = 0.0D+0 DO 3 I=1,N IF(ITOL.EQ.2) THEN VHOLD = DMAX1(RTOL(1),ATOL(1)) ELSE IF (ITOL.EQ.3) THEN VHOLD = DMAX1(RTOL(1),ATOL(I)) ELSE IF (ITOL.EQ.4) THEN VHOLD = DMAX1(RTOL(I),ATOL(1)) ELSE IF (ITOL.EQ.5) THEN VHOLD = DMAX1(RTOL(I),ATOL(I)) ENDIF IF(VHOLD.LE.0.0D+0) THEN WRITE(LOUT,9040) IDID = -4 ENDIF 3 CONTINUE ENDIF IF (N.LE.0) THEN C ILLEGAL VALUE FOR THE NUMBER OF EQUATIONS WRITE (LOUT,9050) IDID = -4 ENDIF IF ((T0-TOUT)*HO.GE.0.0D+0) THEN C PARAMETERS FOR INTEGRATION ARE ILLEGAL WRITE (LOUT,9060) IDID = -4 ENDIF IF ((MF.NE.25).AND.(MF.NE.26))THEN C ILLEGAL VALUE FOR METHOD FLAG WRITE (LOUT,9090) MF IDID = -4 ENDIF IF(ITOL.LT.1.OR.ITOL.GT.5) THEN C ILLEGAL VALUE FOR ERROR CONTROL PARAMETER WRITE (LOUT,9110) IDID=-4 ENDIF IF(MAXDER.LT.1.OR.MAXDER.GT.7) THEN C ILLEGAL VALUE FOR MAXIMUM ORDER WRITE(LOUT,9120) IDID = -4 ENDIF C IF(NIND1.EQ.0) NIND1=N IF(NIND1 + NIND2 + NIND3 .NE. N) THEN C SUM OF VARIABLES OF DIFFERENT INDEX SHOULD BE N. WRITE(LOUT,9140) IDID = -4 ENDIF IF (IDID.NE.1) THEN RETURN ELSE C THE INITIAL PARAMETERS ARE O.K. SO INITIALISE EVERYTHING C ELSE NECESSARY FOR THE INTEGRATION. C IF VALUES OF YMAX OTHER THAN THOSE SET BELOW ARE DESIRED, C THEY SHOULD BE SET HERE. ALL YMAX(I) MUST BE POSITIVE. IF C VALUES FOR HMIN OR HMAX, THE BOUNDS ON DABS(H), OTHER THAN C THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW. IF(ITOL.EQ.1) THEN DO 10 I = 1,N YMAX(I) = DABS(Y0(I)) YMAX(I)=DMAX1(YMAX(I),1.0D0) 10 CONTINUE ENDIF DO 15 I=1,N Y(I,1)=Y0(I) 15 CONTINUE T = T0 H = HO HMIN = DABS(HO) HMAX = DABS(T0-TEND)*10.0D+0 JSTART = 0 NHCUT = 0 C.... IF (MF .EQ. 26) then DO 5 I=1,N AYI = DABS(Y(I,1)) IF(ITOL.EQ.1) THEN SCALE(I) = YMAX(I) ELSE IF(ITOL.EQ.2) THEN SCALE(I) = RTOL(1)*AYI + ATOL(1) ELSE IF(ITOL.EQ.3) THEN SCALE(I) = RTOL(1)*AYI + ATOL(I) ELSE IF(ITOL.EQ.4) THEN SCALE(I) = RTOL(I)*AYI + ATOL(1) ELSE IF(ITOL.EQ.5) THEN SCALE(I) = RTOL(I)*AYI + ATOL(I) ENDIF 5 CONTINUE C... C... COMPUTE THE SPARSITY STRUCTURE OF THE JACOBIAN MATRIX C... CALL SPARSDET(N,T0,Y0,HO,SCALE,HDUMMY,YDOT,IP,NZMAX, & IRNT,F,NZ,JCOL,IROW,W1,IPAR,RPAR,IERR) ELSE NZ = NZMAX ENDIF C... C... CHANGING THE FORMAT OF THE MASS MATRIX TO THE UNCOMPRESS FROMAT C... CALL MAS(N,NZMAS,AM,IAM,JAM) DO J =1, N IAMPT(J) = 0 ENDDO C... C... COUNT TH ENUMBER OF THE NONZER IN EACH ROW C... DO K=1,NZMAS J = JAM(K) IAMPT(J) = IAMPT(J)+1 ENDDO C... C... SET THE IAMPT ARRAY ( WHICH POINTS TO THE FIRST NON-ZERO ELEMENT C... IN EACH ROW OF THE MASS MATRIX) C... KPT = 1 DO J=1,N KR = KPT + IAMPT(J) IAMPT(J) = KPT KPT = KR ENDDO IAMPT(N+1) = KR C... ENDIF ENDIF C... C <<<<<<<<<<<<<<<<< C < TAKE A STEP > C <<<<<<<<<<<<<<<<< 20 IF ((T+H).EQ.T) THEN WRITE (LOUT,9000) ENDIF c CALL STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,T,TOUT,TEND,Y,N, + YMAX,ERRORS,SAVE1,SAVE2,SCALE,PW,PWCOPY,YHOLD,YNHOLD, + ARH,LOUT,MAXDER,ITOL,RTOL,ATOL,F,PDERV,MAS,NQUSED, + NSTEP,NFAIL,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD, + UROUND,EPSJAC,HUSED,ICN,IRN,IROW,JCOL,IKEEP,IW,W,LICN, + LIRN,U,NZ,NZPW,IP,NIND1,NIND2,NIND3,AM,JAM,IAM,NZMAS, + IAMPT,IPAR,RPAR,IERR,HDUMMY,YDOT,Z,IG,W1,HMAX1,IRNT) C c IF(IERR .NE. 0) THEN c WRITE(LOUT,9150) c RETURN c ENDIF c KGO = 1 - KFLAG IF (KGO.EQ.1) THEN C NORMAL RETURN FROM STIFF GO TO 30 ELSE IF (KGO.EQ.2) THEN C COULD NOT ACHIEVE REQUIRED PRECISION WITH HMIN C SO CHOP HMIN IF WE HAVEN'T DONE SO 10 TIMES GO TO 60 ELSE IF (KGO.EQ.3) THEN C ERROR REQUIREMENT SMALLER THAN CAN BE HANDLED FOR THIS PROBLEM WRITE (LOUT,9010) T,H GO TO 70 ELSE IF (KGO.EQ.4) THEN C COULD NOT ACHIEVE CONVERGENCE WITH HMIN WRITE (LOUT,9030) T GO TO 60 ENDIF 30 CONTINUE C --------------------------------------------------------------------- C NORMAL RETURN FROM THE INTEGRATOR. C C THE WEIGHTS YMAX(I) ARE UPDATED IF ITOL=1. C IF DIFFERENT VALUES ARE DESIRED, THEY SHOULD BE SET HERE. C C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY C STEP SHOULD BE INSERTED HERE. C C IF IDID = 3, Y0 IS SET TO THE CURRENT Y VALUES ON RETURN. C IF IDID = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF C ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON RETURN. C FOR ANY OTHER VALUE OF IDID, CONTROL RETURNS TO THE INTEGRATOR C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF Y ARE C COMPUTED AND STORED IN Y0 ON RETURN. C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520. C -------------------------------------------------------------------- IF(NSTEP.GT.MAXSTP) THEN KGO=5 KLAG=4 c TOO MUCH WORK WRITE(LOUT,9130) IDID=-6 GOTO 70 ENDIF IF(ITOL.EQ.1) THEN D = 0.0D+0 DO 40 I = 1,N AYI = DABS(Y(I,1)) YMAX(I)=DMAX1(YMAX(I),AYI) 40 CONTINUE ENDIF IF (IDID.EQ.3.OR.IDID.EQ.1) GO TO 70 IF (DABS(T-TOUT).LE.DABS(15.0D+0*UROUND*TOUT)) THEN C EFFECTIVELY WE HAVE HIT TOUT IDID = KFLAG T0 = TOUT DO 50 I = 1,N Y0(I) = Y(I,1) 50 CONTINUE HO = H RETURN ENDIF IF (IDID.EQ.2) THEN C CONTINUING INTEGRATION BUT MUST HIT TOUT EXACTLY IF (((T+H)-TOUT)*H.GT.0.0D+0) THEN C WE HAVE ALREADY OVERSHOT THE OUTPUT POINT OR WE WILL DO C SO ON THE NEXT STEP IF (((T-TOUT)*H.GE.0.0D+0) .OR. (DABS(T-TOUT).LE. + 100.0D+0*UROUND*HMAX)) THEN C HAVE OVERSHOT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT HO = H IDID = KFLAG RETURN ELSE C WILL PASS TOUT ON NEXT STEP WITH CURRENT STEPSIZE C SO REDUCE STEPSIZE TO HIT TOUT 'EXACTLY' H = (TOUT-T)* (1.0D+0-4.0D+0*UROUND) JSTART = -1 ENDIF ENDIF ELSE IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT, SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG HO = H T0 = TOUT RETURN ENDIF GO TO 20 C ------------------------------------------------------------------- C ON AN ERROR RETURN FROM THE INTEGRATOR, AN IMMEDIATE RETURN OCCURS C IF KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. C H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10 TIMES C BEFORE GIVING UP. C -------------------------------------------------------------------- 60 CONTINUE IF (NHCUT.EQ.10) THEN C HAVE REDUCED H TEN TIMES WRITE (LOUT,9100) GO TO 70 ENDIF NHCUT = NHCUT + 1 HMIN = 0.1D+0*HMIN H = 0.1D+0*H JSTART = -1 GO TO 20 70 IF(DABS(T-TOUT).GT.1000.0D+0*UROUND) THEN DO 80 I = 1,N Y0(I) = Y(I,1) 80 CONTINUE T0 = T ELSE C HAVE PASSED TOUT SO INTERPOLATE CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT IDID = KFLAG ENDIF HO = H IF(KFLAG.NE.0) IDID = KFLAG RETURN C -------------------------- END OF SUBROUTINE OVDRIV ----------------- 9000 FORMAT (' WARNING.. T + H = T ON NEXT STEP.') 9010 FORMAT (/,/,' KFLAG = -2 FROM INTEGRATOR AT T = ',E16.8,' H =', + E16.8,/, + ' THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED',/,/) 9020 FORMAT (/,/,' INTEGRATION HALTED BY DRIVER AT T = ',E16.8,/, + ' ERROR TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION', + /) 9030 FORMAT (/,/,' KFLAG = -3 FROM INTEGRATOR AT T = ',E16.8,/, + ' CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED',/) 9035 FORMAT (/,/,' KFLAG = -4 FROM INTEGRATOR AT T = ',E16.8,/, + ' MA28 FATAL ERROR ',/) 9040 FORMAT (/,/,' ILLEGAL INPUT.. RTOL .LE. 0.',/,/) 9045 FORMAT (/,/,' ILLEGAL INPUT.. ATOL .LE. 0.',/,/) 9050 FORMAT (/,/,' ILLEGAL INPUT.. N .LE. 0',/,/) 9060 FORMAT (/,/,' ILLEGAL INPUT.. (T0-TOUT)*H .GE. 0.',/,/) 9070 FORMAT (/,/,' ILLEGAL INPUT.. IDID =',I5,/,/) 9080 FORMAT (/,/,' IDID = -1 ON INPUT WITH (T-TOUT)*H .GE. 0.',/, + ' T =',E16.8,' TOUT =',E16.8,' H =',E16.8,/, + ' INTERPOLATION WAS DONE AS ON NORMAL RETURN.',/, + ' DESIRED PARAMETER CHANGES WERE NOT MADE.') 9090 FORMAT (/,/,' ILLEGAL INPUT.. METHOD FLAG, MF, = ',I6,/, + ' ALLOWED VALUES ARE 25 OR 26',/) 9100 FORMAT (/,/,' PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT',/, + ' HMIN REDUCED BY A FACTOR OF 1.0E10',/,/) 9110 FORMAT (/,/,' ILLEGAL VALUE FOR ITOL',/,/) 9120 FORMAT (/,/,' ILLEGAL VALUE FOR MAXDER',/,/) 9130 FORMAT (/,/,' NUMBER OF STEPS EXCEEDS MAXIMUM',/,/) 9140 FORMAT (/,/,'BAD INPUT FOR NUMBER OF VARIABLES OF INDEX 1,2,3' + ,/,/) 9150 FORMAT (' IERR IS NOT ZERO BECAUSE OF AN ILLEGAL FUNCTION CALL') 9160 FORMAT (/,/,'STEPSIZE IS TOO SMALL') END C-------------------------------------------------------------------- SUBROUTINE INTERP(N,JSTART,H,T,Y,TOUT,Y0) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER JSTART,N C .. C .. ARRAY ARGUMENTS .. DIMENSION Y(N,12),Y0(N) C .. C .. LOCAL SCALARS .. INTEGER I,J,L C .. C .. INTRINSIC FUNCTIONS .. C .. DO 10 I = 1,N Y0(I) = Y(I,1) 10 CONTINUE L = JSTART + 2 S = (TOUT-T)/H S1 = 1.0D+0 DO 30 J = 2,L S1 = S1* (S+DBLE(FLOAT(J-2)))/DBLE(FLOAT(J-1)) DO 20 I = 1,N Y0(I) = Y0(I) + S1*Y(I,J) 20 CONTINUE 30 CONTINUE RETURN C -------------- END OF SUBROUTINE INTERP --------------------------- END C SUBROUTINE COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C -------------------------------------------------------------------- C COSET IS CALLED BY THE INTEGRATOR AND SETS THE COEFFICIENTS USED C BY THE CONVENTIONAL BACKWARD DIFFERENTIATION SCHEME AND THE C MODIFIED EXTENDED BACKWARD DIFFERENTIATION SCHEME. THE VECTOR C EL OF LENGTH NQ+1 DETERMINES THE BASIC BDF METHOD WHILE THE VECTOR C ELST OF LENGTH NQ+2 DETERMINES THE MEBDF. THE VECTOR TQ OF C LENGTH 4 IS INVOLVED IN ADJUSTING THE STEPSIZE IN RELATION TO THE C TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE PERTST ARRAY. THE C VECTORS EL AND TQ BOTH DEPEND ON METH AND NQ. THE C COEFFICIENTS IN PERTST NEED TO BE GIVEN TO ONLY ABOUT ONE PERCENT C ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS: C COEFFICIENTS FOR ORDER NQ-1, COEFFICIENTS FOR ORDER NQ, C COEFFICIENTS FOR ORDER NQ+1. C ------------------------------------------------------------------- C .. SCALAR ARGUMENTS .. INTEGER NQ C .. C .. ARRAY ARGUMENTS .. DIMENSION EL(10),ELST(10),TQ(5) C .. C .. LOCAL SCALARS .. INTEGER K C .. C .. LOCAL ARRAYS .. DIMENSION PERTST(8,3) C .. C .. INTRINSIC FUNCTIONS .. C .. C .. COMMON BLOCKS .. INTEGER MAXORD,NCOSET C .. C .. DATA STATEMENTS .. DATA PERTST(1,1)/1./,PERTST(2,1)/2./,PERTST(3,1)/4.5/, + PERTST(4,1)/7.333/,PERTST(5,1)/10.42/, + PERTST(6,1)/13.7/,PERTST(7,1)/17.15/, + PERTST(8,1)/20.74/ DATA PERTST(1,2)/2./,PERTST(2,2)/4.5/, + PERTST(3,2)/7.333/,PERTST(4,2)/10.42/, + PERTST(5,2)/13.7/,PERTST(6,2)/17.15/, + PERTST(7,2)/20.74/,PERTST(8,2)/24.46/ DATA PERTST(1,3)/4.5/,PERTST(2,3)/7.333/, + PERTST(3,3)/10.42/,PERTST(4,3)/13.7/, + PERTST(5,3)/17.15/,PERTST(6,3)/20.74/, + PERTST(7,3)/24.46/,PERTST(8,3)/1./ C .. C ------------------------------------------------------------------- C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY. C THEIR DERIVATION IS GIVEN IN REFERENCE 1. C ------------------------------------------------------------------- IF (NQ.GT.MAXORD) MAXORD = NQ NCOSET = NCOSET + 1 GO TO (10,20,30,40,50,60,70) NQ 10 EL(1) = 1.0D+0 ELST(1) = 1.5D+0 ELST(3) = -0.5D+0 GO TO 80 20 EL(1) = 6.6666666666667D-01 EL(3) = 3.3333333333333D-01 ELST(1) = 9.5652173913043D-01 ELST(3) = 2.1739130434782D-01 ELST(4) = -1.7391304347826D-01 GO TO 80 30 EL(1) = 5.4545454545455D-01 EL(3) = 4.5454545454545D-01 EL(4) = 1.8181818181818D-01 ELST(1) = 7.6142131979695D-01 ELST(3) = 3.2994923857868D-01 ELST(4) = 8.6294416243654D-02 ELST(5) = -9.1370558375634D-02 GO TO 80 40 EL(1) = 0.48D+0 EL(3) = 0.52D+0 EL(4) = 0.28D+0 EL(5) = 0.12D+0 ELST(1) = 6.5733706517393D-01 ELST(3) = 4.0023990403838D-01 ELST(4) = 1.5793682526989D-01 ELST(5) = 4.4382247101159D-02 ELST(6) = -5.7576969212315D-02 GO TO 80 50 EL(1) = 4.3795620437956D-01 EL(3) = 5.62043795620436D-01 EL(4) = 3.43065693430656D-01 EL(5) = 1.97080291970802D-01 EL(6) = 8.75912408759123D-02 ELST(1) = 5.9119243917152D-01 ELST(3) = 4.4902473356122D-01 ELST(4) = 2.1375427307460D-01 ELST(5) = 9.0421610027481503D-02 ELST(6) = 2.6409276761177D-02 ELST(7) = -4.0217172732757D-02 GO TO 80 60 EL(1) = 4.08163265306120D-01 EL(3) = 5.91836734693874D-01 EL(4) = 3.87755102040813D-01 EL(5) = 2.51700680272107D-01 EL(6) = 1.49659863945577D-01 EL(7) = 6.80272108843534D-02 ELST(1) = 5.4475876041119D-01 ELST(3) = 4.8525549636077D-01 ELST(4) = 2.5789750131312D-01 ELST(5) = 1.3133738525800D-01 ELST(6) = 5.7677396763462D-02 ELST(7) = 1.7258197643881D-02 ELST(8) = -3.0014256771967D-02 GO TO 80 70 EL(1) = 3.85674931129476D-01 EL(3) = 6.14325068870521D-01 EL(4) = 4.21487603305783D-01 EL(5) = 2.9292929292929D-01 EL(6) = 1.96510560146923D-01 EL(7) = 1.19375573921028D-01 EL(8) = 5.50964187327820D-02 ELST(1) = 5.0999746293734D-01 ELST(3) = 5.1345839935281D-01 ELST(4) = 2.9364346131937D-01 ELST(5) = 1.6664672120553D-01 ELST(6) = 8.8013735242353D-02 ELST(7) = 3.9571794884069D-02 ELST(8) = 1.2039080338722D-02 ELST(9) = -2.3455862290154D-02 80 DO 90 K = 1,3 TQ(K) = PERTST(NQ,K) 90 CONTINUE TQ(4) = 0.5D+0*TQ(2)/DBLE(FLOAT(NQ)) IF(NQ.NE.1) TQ(5)=PERTST(NQ-1,1) RETURN C --------------------- END OF SUBROUTINE COSET --------------------- END C SUBROUTINE PSET(Y,N,H,T,UROUND,EPSJAC,CON,MITER,IER,F,PDERV, + MAS,NRENEW,YMAX,SAVE1,SAVE2,PW,PWCOPY,WRKSPC,ITOL,RTOL, + ATOL,NPSET,NJE,NFE,NDEC,ICN,IRN,IROW,JCOL,IKEEP,IW,W, + LICN,LIRN,NZ,NZPW,U,IP,NIND1,NIND2,NIND3,AM,JAM,IAM, + NZMAS,IPAR,RPAR,IERR,HDUMMY,YDOT,Z,IG,W1,HMAX1,IRNT) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------- C PSET IS CALLED BY STIFF TO COMPUTE AND PROCESS THE MATRIX C I/(H*EL(1)) - J WHERE J IS AN APPROXIMATION TO THE RELEVANT JACOBIAN C AND I IS THE IDENTITY MATRIX. THIS MATRIX IS THEN SUBJECTED TO LU C DECOMPOSITION IN PREPARATION FOR LATER SOLUTION OF LINEAR SYSTEMS C OF ALGEBRAIC EQUATIONS WITH LU AS THE COEFFICIENT MATRIX. C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH C PSET USES THE FOLLOWING .. C EPSJAC = DSQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS. C ******************************************************************* C THE ARGUMENT NRENEW IS USED TO SIGNAL WHETHER OR NOT C WE REQUIRE A NEW JACOBIAN TO BE CALCULATED. C C IF NRENEW > 0 THEN WE REQUIRE A NEW J TO BE COMPUTED C = 0 THEN USE A COPY OF THE LAST J COMPUTED C ******************************************************************* C THE FLAG IER IS SET TO 1 WHEN A STRUCTURAL SINGULARITY OCCURS C AND TO 2 WHEN A NUMERICAL SINGULARITY OCCURS C ******************************************************************* C C .. SCALAR ARGUMENTS .. INTEGER IER,MITER,N,NRENEW C .. C .. ARRAY ARGUMENTS .. DIMENSION PW(*),PWCOPY(*),SAVE2(N),WRKSPC(N),Y(N,12),HMAX1(n), & YMAX(N),SAVE1(N),RTOL(1),ATOL(1),W(1),RPAR(*),IPAR(*) INTEGER IROW(NZPW),JCOL(NZPW),ICN(LICN),IRN(LIRN), & IKEEP(N,5),IW(N,8),IRN1(1),IP(N+1),IG(2*N+1), & IAM(NZMAS),JAM(NZMAS),IRNT(*) DOUBLE PRECISION YDOT(N),W1(2*N),HDUMMY(N),Z(N),AM(*) C LOGICAL MA28B C C .. EXTERNAL SUBROUTINES .. EXTERNAL F,PDERV,MA28AD,MA28BD,TD02AD,TD02BD,MAS C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSQRT C .. C .. COMMON BLOCKS .. COMMON /MA28FL/MA28B INTEGER NDEC,NFE,NJE,NPSET C... NPSET = NPSET + 1 IF (NRENEW.EQ.0) THEN IF ((MITER.EQ.5) .OR. (MITER .EQ. 6)) THEN DO I = 1,NZ PW(I) = -PWCOPY(I) ENDDO ENDIF GO TO 70 ENDIF C NJE = NJE + 1 IF (MITER .EQ. 5 ) THEN CALL PDERV(N,T,Y,NZ,PWCOPY,IROW,JCOL,IPAR,RPAR,IERR) DO 22 I = 1,NZ PW(I) = - PWCOPY(I) 22 CONTINUE C... c... NUMERICAL APRROXIMATION TO THE SPARSE JACOBIAN c... ELSE D = 0.D0 DO 41 I = 1,N D = D + (SAVE2(I))**2 41 CONTINUE IF (ITOL .EQ. 1) D = D*RTOL(1)**2 R0 = DABS(H)*DSQRT(D)*1.0D+03*UROUND IF (R0.EQ.0.D0) R0 = 1.D0 DO 25 I = 1,N R = EPSJAC*DABS(YMAX(I)) HDUMMY(I) = DMAX1(R,R0,1.D-11) 25 HMAX1(I) = 0.1D0 IF (NJE .EQ. 1) THEN IG(1) = 0 ENDIF CALL TD02AD(N,N,IRNT,IP,F,HDUMMY,Y,T,SAVE2, & HMAX1,PWCOPY,IG,W1,Z,IPAR,RPAR,IERR) NFE = NFE + W1(1) DO 30 I = 1,NZ PW(I) = - PWCOPY(I) 30 CONTINUE ENDIF C... 70 CONTINUE CALL MAS(N,NZMAS,AM,IAM,JAM) C... NZPW = NZ + NZMAS DO I= NZ+1,NZPW PW(I) = AM(I-NZ) /CON IROW(I) = IAM(I-NZ) JCOL(I) = JAM(I-NZ) ENDDO C... IF ((MITER .EQ. 5).OR.(MITER.EQ.6)) THEN IF ( MA28B ) THEN CALL MA28BD(N,NZPW,PW,LICN,IROW,JCOL,ICN,IKEEP,IW, & W,IFLAG) C C.....CHECK IF MATRIX IS STRUCTURALLY SINGULAR IF (IFLAG.EQ.-1) THEN WRITE(6, *) 'STRUCTURALLY SINGULAR MATRIX' IER=1 MA28B=.FALSE. RETURN ELSE IF (IFLAG.LT.0.AND.IFLAG.NE.-14. & AND.IFLAG.NE.-2) THEN WRITE(6,*) 'MA28 - FATAL ERROR, IFLAG:',IFLAG STOP ENDIF IF ( IFLAG.EQ.-2 ) MA28B = .FALSE. ENDIF IF ( .NOT.MA28B ) THEN DO I=1,NZPW IRN(I)=IROW(I) ICN(I)=JCOL(I) ENDDO CALL MA28AD(N,NZPW,PW,LICN,IRN,LIRN,ICN,U,IKEEP, & IW,W,IFLAG) IF ( IFLAG.EQ.0.OR.IFLAG.EQ.-14 ) THEN MA28B=.TRUE. ELSE IF ( IFLAG.EQ.-1 ) THEN WRITE(6,*) 'STRUCTURAL SINGULARITY, IFLAG:',IFLAG IER=1 ELSE IF (IFLAG.EQ.-2 ) THEN WRITE(6,*) 'NUMERICAL SINGULARITY, AT TIME:',T IER=2 ELSE IF ( IFLAG.LT.-2 ) THEN WRITE(6,*) 'MA28 - FATAL ERROR, IFLAG:',IFLAG STOP ENDIF ENDIF ENDIF NDEC = NDEC + 1 ENDIF RETURN END C ---------------------- END OF SUBROUTINE PSET --------------------- c SUBROUTINE ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C *************************************************** C C THIS ROUTINE CALCULATES ERRORS USED IN TESTS C IN STIFF . C C *************************************************** C .. SCALAR ARGUMENTS .. INTEGER N C .. C .. ARRAY ARGUMENTS .. DIMENSION TQ(5) C .. C .. LOCAL SCALARS .. C .. C .. INTRINSIC FUNCTIONS .. C .. SQHOL = DBLE(FLOAT(N)) EDN = TQ(1)*TQ(1)*SQHOL C C ** ERROR ASSOCIATED WITH METHOD OF ORDER ONE LOWER. C E = TQ(2)*TQ(2)*SQHOL C C ** ERROR ASSOCIATED WITH PRESENT ORDER C EUP = TQ(3)*TQ(3)*SQHOL C C ** ERROR ASSOCIATED WITH HIGHER ORDER METHOD C BND = TQ(4)*TQ(4)*SQHOL*0.5D+0 EDDN=TQ(5)*TQ(5)*SQHOL C C ** ERROR ASSOCIATED WITH METHOD OF ORDER TWO LOWER. RETURN END C---------------------------------------------------------------------- SUBROUTINE PRDICT(T,H,Y,L,N,YPRIME,NFE,F,IPAR,RPAR,IERR) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ********************************************************************** C PREDICTS A VALUE FOR Y AT (T+H) GIVEN THE HISTORY ARRAY AT T C THEN EVALUATES THE DERIVATIVE AT THIS POINT, THE RESULT OF THIS C EVALUATION BEING STORED IN YPRIME() C ********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER L,N C .. C .. ARRAY ARGUMENTS .. DIMENSION Y(N,12),YPRIME(N),RPAR(*),IPAR(*) C .. C .. LOCAL SCALARS .. INTEGER I,J2 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL F C .. C .. COMMON BLOCKS .. INTEGER NFE C .. DO 20 J2 = 2,L DO 10 I = 1,N Y(I,1) = Y(I,1) + Y(I,J2) 10 CONTINUE 20 CONTINUE T = T + H CALL F(N,T,Y,YPRIME,IPAR,RPAR,IERR) NFE = NFE + 1 RETURN END C C---------------------------------------------------------------------- C SUBROUTINE ITRAT2(QQQ,Y,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE, + M,WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,LMB,ITOL,RTOL, + ATOL,HUSED,NBSOL,NFE,NQUSED,F,LICN,U,ICN,IKEEP,W, + NIND1,NIND2,NIND3,AM,JAM,IAM,IAMPT,NZMAS,IPAR,RPAR,IERR) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER M,N,ITOL LOGICAL WORKED C .. C .. ARRAY ARGUMENTS .. DIMENSION ARH(N),ERROR(N),PW(*),SAVE1(N),SAVE2(N),Y(N,12), & YMAX(N),RTOL(1),ATOL(1),SCALE(N),W(N),ICN(LICN), & IKEEP(N,5),AM(*),JAM(*),IAM(*),IPAR(*),RPAR(*),IAMPT(N+1) C C .. EXTERNAL SUBROUTINES .. EXTERNAL F, MA28CD C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DMAX1,DMIN1 C .. C .. COMMON BLOCKS .. INTEGER NBSOL,NFE C .. C .. DATA STATEMENTS .. DATA ZERO/0.0D+0/ C .. DO 5 I=1,N AYI = DABS(Y(I,1)) IF(ITOL.EQ.1) THEN SCALE(I) = YMAX(I) ELSE IF(ITOL.EQ.2) THEN SCALE(I) = RTOL(1)*AYI + ATOL(1) ELSE IF(ITOL.EQ.3) THEN SCALE(I) = RTOL(1)*AYI + ATOL(I) ELSE IF(ITOL.EQ.4) THEN SCALE(I) = RTOL(I)*AYI + ATOL(1) ELSE IF(ITOL.EQ.5) THEN SCALE(I) = RTOL(I)*AYI + ATOL(I) ENDIF 5 CONTINUE IF(NIND2.NE.0) THEN DO 11 I = NIND1+1,NIND2+NIND1 SCALE(I)=SCALE(I)/HUSED 11 CONTINUE ENDIF IF(NIND3.NE.0) THEN DO 12 I = NIND1 +NIND2 + 1,NIND3+NIND2+NIND1 SCALE(I)=SCALE(I)/(HUSED**2) 12 CONTINUE ENDIF IF(LMB.EQ.1) GOTO 25 DO 14 I=1,N SAVE1(I) = (-SAVE1(I) + ARH(I)) 14 CONTINUE DO 9911 I=1,N SAVE2(I) = SAVE2(I)*HBETA DO K = IAMPT(I), IAMPT(I+1) -1 SAVE2(I) = SAVE2(I) + AM(K)*SAVE1(JAM(K)) ENDDO 9911 CONTINUE DO 8812 I=1,N SAVE1(I) = SAVE2(I)/QQQ 8812 CONTINUE C IF (( MF .EQ. 25 ) .OR. (MF .EQ. 26))THEN CALL MA28CD(N,PW,LICN,ICN,IKEEP,SAVE1,W,1) NBSOL = NBSOL + 1 ENDIF D = ZERO DO 20 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/(SCALE(I)))**2 SAVE1(I) = Y(I,1) + ERROR(I) 20 CONTINUE IF(ITOL.EQ.1) D=D/(RTOL(1)**2) TCRATE = TCRATE + CRATE D1 = D M = 1 CALL F(N,T,SAVE1,SAVE2,IPAR,RPAR,IERR) NFE = NFE + 1 25 CONTINUE WORKED = .TRUE. 30 CONTINUE C DO 44 I=1,N SAVE1(I) = (-SAVE1(I) + ARH(I)) 44 CONTINUE DO 41 I=1,N SAVE2(I) = SAVE2(I)*HBETA DO K= IAMPT(I), IAMPT(I+1)-1 SAVE2(I) = SAVE2(I) + AM(K)*SAVE1(JAM(K)) ENDDO 41 CONTINUE DO 42 I=1,N SAVE1(I) = SAVE2(I)/QQQ 42 CONTINUE C... C... IF WE ARE HERE THEN PARTIALS ARE O.K. C... IF(( MF.EQ. 25 ) .OR. (MF .EQ. 26)) THEN CALL MA28CD(N,PW,LICN,ICN,IKEEP,SAVE1,W,1) NBSOL=NBSOL + 1 ENDIF C C WE NOW CALCULATE A WEIGHTED RMS TYPE NORM C D = ZERO DO 50 I = 1,N ERROR(I) = ERROR(I) + SAVE1(I) D = D + (SAVE1(I)/(SCALE(I)))**2 SAVE1(I) = Y(I,1) + ERROR(I) 50 CONTINUE IF(ITOL.EQ.1) D=D/(RTOL(1)**2) C ------------------------------------------------------------------- C TEST FOR CONVERGENCE. IF M.GT.0 , AN ESTIMATE OF THE CONVERGENCE C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. C ------------------------------------------------------------------- IF (M.NE.0) THEN IF (D1.NE.ZERO) CRATE = DMAX1(0.9D+0*CRATE,D/D1) ENDIF TCRATE = TCRATE + CRATE IF ((D*DMIN1(1.0D+0,2.0D+0*CRATE)).LT.ERRBND/DFLOAT(NQUSED)) + RETURN IF (M.NE.0) THEN IF (D.GT.D1) THEN WORKED = .FALSE. RETURN ENDIF ENDIF D1 = D IF (M.EQ.4) THEN WORKED = .FALSE. RETURN ENDIF M = M + 1 CALL F(N,T,SAVE1,SAVE2,IPAR,RPAR,IERR) NFE = NFE + 1 GO TO 30 END C... C---------------------------------------------------------------------- C SUBROUTINE STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,T,TOUT,TEND,Y,N, + YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,PWCOPY,YHOLD,YNHOLD,ARH, + LOUT,MAXDER,ITOL,RTOL,ATOL,F,PDERV,MAS,NQUSED,NSTEP,NFAIL, + NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,UROUND,EPSJAC, + HUSED,ICN,IRN,IROW,JCOL,IKEEP,IW,W,LICN,LIRN,U,NZ,NZPW, + IP,NIND1,NIND2,NIND3,AM,JAM,IAM,NZMAS,IAMPT,IPAR,RPAR,IERR, + HDUMMY,YDOT,Z,IG,W1,HMAX1,IRNT) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------ C THE SUBROUTINE STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN C INITIAL VALUE PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL C EQUATIONS OR LINEARLY IMPLICIT DIFFERENTIAL ALGEBRAIC EQUATIONS. C COMMUNICATION WITH STIFF IS DONE WITH THE FOLLOWING VARIABLES.. C Y AN N BY LMAX+3 ARRAY CONTAINING THE DEPENDENT VARIABLES C AND THEIR BACKWARD DIFFERENCES. MAXDER (=LMAX-1) IS THE C MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. C Y(I,J+1) CONTAINS THE JTH BACKWARD DIFFERENCE OF Y(I) C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. C H THE STEPSIZE TO BE ATTEMPTED ON THE NEXT STEP. C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING C THE PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE BUT C ITS SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. C HMIN THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEPSIZE C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY C TIME BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. C RTOL,ATOL THE ERROR BOUNDS. SEE DESCRIPTION IN OVDRIV. C UROUND THE UNIT ROUNDOFF OF THE MACHINE. C N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. C MF THE METHOD FLAG. MUST BE SET TO 21,22,23 OR 24 AT PRESENT C KFLAG A COMPLETION FLAG WITH THE FOLLOWING MEANINGS.. C 0 THE STEP WAS SUCCESSFUL C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED C WITH ABS(H) = HMIN. C -2 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR ABS(H)=HMIN. C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND C THE Y ARRAY ARE AS AT THE BEGINNING OF THE LAST C STEP ATTEMPTED, AND H IS THE LAST STEP SIZE ATTEMPTED. C JSTART AN INTEGER USED ON INPUT AND OUTPUT. C ON INPUT IT HAS THE FOLLOWING VALUES AND MEANINGS.. C 0 PERFORM THE FIRST STEP. C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF H OR N. C ON EXIT, JSTART IS NQUSED, THE ORDER OF THE METHOD LAST USED. C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL C ERRORS IN Y ARE COMPARED C ERROR AN ARRAY OF N ELEMENTS. C SAVE1,2 TWO ARRAYS OF WORKING SPACE BOTH OF LENGTH N. C PW A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES C IPIV AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT INFORMATION. C JNEWIM IS TO INDICATE IF PRESENT ITERATION MATRIX C WAS FORMED USING A NEW J OR OLD J. C JSNOLD KEEPS TRACK OF NO. OF STEPS TAKEN WITH C PRESENT ITERATION MATRIX (BE IT FORMED BY C A NEW J OR NOT). C AVNEWJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION C MATRIX WAS FORMED BY A NEW J. C AVOLDJ STORES VALUE FOR AVERAGE CRATE WHEN ITERATION C MATRIX WAS FORMED BY AN OLD J. C NRENEW FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE PSET. C IF NRENEW > 0 THEN FORM A NEW JACOBIAN BEFORE C COMPUTING THE COEFFICIENT MATRIX FOR C THE NEWTON-RAPHSON ITERATION C = 0 FORM THE COEFFICIENT MATRIX USING A C COPY OF AN OLD JACOBIAN C NEWPAR FLAG USED IN THIS SUBROUTINE TO INDICATE IF A JACOBIAN C HAS BEEN EVALUATED FOR THE CURRENT STEP C ********************************************************************** C .. SCALAR ARGUMENTS .. INTEGER JSTART,KFLAG,LOUT,MF,N,ITOL,NZ C .. C .. ARRAY ARGUMENTS .. DIMENSION HSTPSZ(2,14) DIMENSION ARH(N),ERROR(N),PW(*),PWCOPY(*),SAVE1(N),SAVE2(N), + Y(N,12), YHOLD(N,12),YMAX(N),YNHOLD(N,2),RTOL(1),ATOL(1), + SCALE(N),ICN(LICN),IRN(LIRN),IKEEP(N,5),W(N),IROW(NZPW), + JCOL(NZPW),IW(N,8),IP(*),AM(*),IAM(*),JAM(*), + IPAR(*),RPAR(*),IAMPT(N+1),IRNT(*) C .. C .. LOCAL SCALARS .. INTEGER I,IER,IITER,IITER2,IREDO,ITST,J,J1,J2,KFAIL,LL,LMAX, + M3STEP,MAXDER,METH,MFOLD,MITER,NEWQ,LICN,LIRN,NZPW LOGICAL FINISH,OVRIDE,WORKED C .. C .. LOCAL ARRAYS .. DIMENSION EL(10),ELST(10),TQ(5) C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL COSET,CPYARY,ERRORS,F,PDERV,HCHOSE,ITRAT2, + PRDICT,PSET,RSCALE,MA28CD,MAS C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DMIN1 C .. C .. COMMON BLOCKS .. COMMON/STPSZE/HSTPSZ INTEGER IDOUB,ISAMP,IWEVAL,JCHANG,JSINUP,JSNOLD,L,M1,M2, + MAXORD, MEQC1,MEQC2,MQ1TMP,MQ2TMP,NBSOL,NCOSET, + NDEC,NEWPAR,NFE,NJE,NPSET,NQ,NQUSED,NRENEW,NSTEP C LOGICAL CFAIL,JNEWIM,SAMPLE C .. C .. SAVE STATEMENT .. C .. SAVE MFOLD,HOLD,LMAX,EDN,EUP,BND,EDDN,EL,TQ,ELST,E,IWEVAL, + NEWPAR,NRENEW,JSINUP,JSNOLD,IDOUB,JCHANG,L,NQ,MEQC1,QI,QQQ, + MEQC2,MQ1TMP,MQ2TMP,ISAMP,RH,RMAX,TOLD,CRATE1,CRATE2,IER, + TCRAT1,TCRAT2,AVNEW2,AVOLD2,AVNEWJ,AVOLDJ,MITER,IBND,UPBND, + CFAIL,KFAIL,RC,JNEWIM,VTOL,SAMPLE,IEMB,OLDLO,IREDO,OVRIDE C C .. DATA STATEMENTS .. DATA EL(2),ELST(2),OLDLO/3*1.0D+0/ DATA ZERO,ONE/0.0D+0,1.0D+0/ C .. TOLD = T KFLAG = 0 IF (JSTART.GT.0) GO TO 60 IF (JSTART.NE.0) GO TO 30 C ------------------------------------------------------------------ C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT C IS CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE C INCREASED IN A SINGLE STEP. RMAX IS SET EQUAL TO 1.D4 INITIALLY C TO COMPENSATE FOR THE SMALL INITIAL H, BUT THEN IS NORMALLY = 10. C IF A FAILURE OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), C RMAX IS SET AT 2. FOR THE NEXT INCREASE. C ------------------------------------------------------------------ CALL F(N,T,Y,SAVE1,IPAR,RPAR,IERR) IF (IERR .NE. 0) GOTO 8000 DO 10 I = 1,N Y(I,2) = H*SAVE1(I) 10 CONTINUE METH = 2 MITER = MF - 10*METH IBND=5 UPBND=0.1D+0 NQ = 1 NQUSED=NQ IER=0 L = 2 IDOUB = 3 KFAIL = 0 RMAX = 10000.0D+0 RC = ZERO CRATE1 = 0.1D+0 CRATE2 = 0.1D+0 JSNOLD = 0 JNEWIM = .TRUE. TCRAT1 = ZERO TCRAT2 = ZERO VTOL=DMAX1(RTOL(1),ATOL(1))/10.0D+0 DO 15 I=1,12 HSTPSZ(1,I)=1.0D+0 HSTPSZ(2,I)=VTOL 15 CONTINUE HOLD = H MFOLD = MF NSTEP = 0 NFE = 1 NJE = 0 NDEC = 0 NPSET = 0 NCOSET = 0 MAXORD = 1 NFAIL = 0 CFAIL = .TRUE. AVNEWJ = ZERO AVOLDJ = ZERO AVNEW2 = ZERO AVOLD2 = ZERO SAMPLE = .FALSE. ISAMP = 0 IEMB = 0 C ************************************************** C CFAIL=.TRUE. ENSURES THAT WE CALCULATE A NEW C J ON THE FIRST CALL. C ************************************************** MEQC1 = 0 MEQC2 = 0 MQ1TMP= 0 MQ2TMP= 0 NBSOL = 0 HUSED = H C ----------------------------------------------------------------- C IF THE CALLER HAS CHANGED N , THE CONSTANTS E, EDN, EUP C AND BND MUST BE RESET. E IS A COMPARISON FOR ERRORS AT THE C CURRENT ORDER NQ. EUP IS TO TEST FOR INCREASING THE ORDER, C EDN FOR DECREASING THE ORDER. BND IS USED TO TEST FOR CONVERGENCE C OF THE CORRECTOR ITERATES. IF THE CALLER HAS CHANGED H, Y MUST C BE RE-SCALED. IF H IS CHANGED, IDOUB IS SET TO L+1 TO PREVENT C FURTHER CHANGES IN H FOR THAT MANY STEPS. C ----------------------------------------------------------------- CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) IWEVAL = MITER NRENEW = 1 NEWPAR = 0 C ***************************************************** C NRENEW AND NEWPAR ARE TO INSTRUCT ROUTINE THAT C WE WISH A NEW J TO BE CALCULATED FOR THIS STEP. C ***************************************************** CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) DO 20 I = 1,N ARH(I) = EL(2)*Y(I,1) 20 CONTINUE CALL CPYARY(N*L,Y,YHOLD) QI = H*EL(1) QQ = ONE/QI CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,F,IPAR,RPAR,IERR) IF (IERR .NE. 0) GOTO 8000 GO TO 110 C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DIFFERENT PARAMETERS ON THIS CALL < C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 30 CALL CPYARY(N*L,YHOLD,Y) IF (MF.NE.MFOLD) THEN METH = MF/10 MITER = MF - 10*METH MFOLD = MF IWEVAL = MITER ENDIF IF(NSTEP.GT.0) GOTO 35 NJE = 0 NFE = 1 CFAIL = .TRUE. NEWPAR = 0 MQ1TMP = 0 MQ2TMP = 0 MEQC1 = 0 MEQC2 = 0 TCRAT1 = 0.0D+0 TCRAT2 = 0.0D+0 CRATE1 = 1.0D-1 CRATE2 = 1.0D-1 NSTEP = 0 NBSOL = 0 NPSET = 0 NCOSET = 0 NDEC = 0 35 CONTINUE IF (H.NE.HOLD) THEN RH = H/HOLD H = HOLD IREDO = 3 GO TO 50 ELSE GO TO 60 ENDIF C ********************************************* C RE-SCALE Y AFTER A CHANGE OF STEPSIZE * C ********************************************* 40 RH = DMAX1(RH,HMIN/DABS(H)) 50 RH = DMIN1(RH,HMAX/DABS(H),RMAX) CALL RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH RC = RC*RH C IF (JSNOLD.GT.IBND) THEN CFAIL = .TRUE. NEWPAR = 0 RC = ZERO C ********************************************************************** C CFAIL=TRUE AND NEWPAR=0 SHOULD FORCE A NEW J TO BE EVALUATED C AFTER 7 STEPS WITH AN OLD J, IF WE HAVE HAD A FAILURE OF ANY C KIND ON THE FIRST, SECOND OR THIRD STAGE OF THE CURRENT STEP C ********************************************************************** ENDIF IDOUB = L + 1 CALL CPYARY(N*L,Y,YHOLD) 60 IF (DABS(RC-ONE).GT.UPBND) IWEVAL = MITER HUSED = H C ------------------------------------------------------------------ C THIS SECTION COMPUTES THE PREDICTED VALUES OF Y C AND THE RHS, ARH, FOR USE IN THE NEWTON ITERATION SCHEME. C RC IS THE RATIO OF THE NEW TO OLD VALUES OF THE COEFFICIENT C H*EL(1). WHEN RC DIFFERS FROM 1 BY MORE THAN 20 PERCENT, IWEVAL IS C SET TO MITER TO FORCE THE PARTIALS TO BE UPDATED. C ------------------------------------------------------------------ QI = H*EL(1) QQ = ONE/QI DO 70 I = 1,N ARH(I) = EL(2)*Y(I,1) 70 CONTINUE DO 90 J1 = 2,NQ JP1 =J1+1 DO 80 I = 1,N ARH(I) = ARH(I) + EL(JP1)*Y(I,J1) 80 CONTINUE 90 CONTINUE IF (JCHANG.EQ.1) THEN C IF WE HAVE CHANGED STEPSIZE THEN PREDICT A VALUE FOR Y(T+H) C AND EVALUATE THE DERIVATIVE THERE (STORED IN SAVE2()) CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,F,IPAR,RPAR,IERR) IF (IERR .NE. 0) GOTO 8000 ELSE C ELSE USE THE VALUES COMPUTED FOR THE SECOND BDF FROM THE LAST C STEP. Y( ,LMAX+2) HOLDS THE VALUE FOR THE DERIVATIVE AT (T+H) C AND Y( ,LMAX+3) HOLDS THE APPROXIMATION TO Y AT THIS POINT. LMP2=LMAX+2 LMP3=LMAX+3 DO 100 I = 1,N SAVE2(I) = Y(I,LMP2) Y(I,1) = Y(I,LMP3) 100 CONTINUE T = T + H ENDIF 110 IF (IWEVAL.LE.0) GO TO 120 C ------------------------------------------------------------------- C IF INDICATED, THE MATRIX P = I/(H*EL(2)) - J IS RE-EVALUATED BEFORE C STARTING THE CORRECTOR ITERATION. IWEVAL IS SET = 0 TO INDICATE C THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN PSET. C THE PROCESSED MATRIX IS STORED IN PW C ------------------------------------------------------------------- IWEVAL = 0 RC = ONE IITER = MEQC1 - MQ1TMP IITER2 = MEQC2 - MQ2TMP IF (JNEWIM) THEN IF (JSNOLD.GE.3) THEN AVNEWJ = TCRAT1/DBLE(FLOAT(IITER)) AVNEW2 = TCRAT2/DBLE(FLOAT(IITER2)) ELSE AVNEWJ = ONE AVNEW2 = ONE END IF ELSE C C MATRIX P WAS FORMED WITH A COPY OF J C IF (JSNOLD.GE.3) THEN AVOLDJ = TCRAT1/DBLE(FLOAT(IITER)) AVOLD2 = TCRAT2/DBLE(FLOAT(IITER2)) IF (AVOLDJ.LT.AVNEWJ) THEN AVNEWJ = AVOLDJ ELSE IF (((DABS(AVOLDJ-AVNEWJ)).GT.0.3D+0) .OR. + ((AVOLDJ.GT.0.85D+0).AND. (AVOLDJ.NE.ONE))) THEN C C SINCE IN CERTAIN INSTANCES AVOLDJ WILL C BE 1.0 AND THERE WILL BE NO NEED TO C UPDATE J. C CFAIL = .TRUE. CRATE1 = 0.1D+0 CRATE2 = 0.1D+0 END IF ELSE CFAIL = .TRUE. CRATE1 = 0.1D+0 CRATE2 = 0.1D+0 C C ********************************************************* C IF WE HAVE REACHED HERE THINGS MUST HAVE GONE WRONG C ********************************************************* C END IF END IF TCRAT1 = ZERO TCRAT2 = ZERO IF (CFAIL) then NRENEW = 1 NEWPAR = 1 JSINUP = -1 JNEWIM = .TRUE. endif CFAIL = .FALSE. JSNOLD = 0 MQ1TMP = MEQC1 MQ2TMP = MEQC2 C... CALL PSET(Y,N,H,T,UROUND,EPSJAC,QI,MITER,IER,F,PDERV,MAS, + NRENEW,YMAX,SAVE1,SAVE2,PW,PWCOPY,ERROR,ITOL,RTOL,ATOL, + NPSET,NJE,NFE,NDEC,ICN,IRN,IROW,JCOL,IKEEP,IW,W,LICN, + LIRN,NZ,NZPW,U,IP,NIND1,NIND2,NIND3,AM,JAM,IAM,NZMAS, + IPAR,RPAR,IERR,HDUMMY,YDOT,Z,IG,W1,HMAX1,IRNT) IF (IERR .NE. 0 ) GOTO 8000 QQQ=QI C C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET IF (IER.NE.0) THEN C IF IER>0 THEN WE HAVE HAD A SINGULARITY IN THE ITERATION MATRIX IJUS=1 RED=0.5D+0 NFAIL = NFAIL + 1 GO TO 450 ENDIF 120 DO 130 I = 1,N SAVE1(I) = Y(I,1) ERROR(I) = ZERO 130 CONTINUE M1 = 0 C ********************************************************************** C UP TO 4 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS MADE C ON THE R.M.S. NORM OF EACH CORRECTION ,USING BND, WHICH DEPENDS C ON ATOL AND RTOL. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE C VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C ********************************************************************** IF (.NOT.SAMPLE) THEN CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1, + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,1, + ITOL,RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,LICN, + U,ICN,IKEEP,W,NIND1,NIND2,NIND3,AM,JAM,IAM,IAMPT, + NZMAS,IPAR,RPAR,IERR) C IF (IERR .NE. 0) GOTO 8000 ITST = 2 ELSE CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1, + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,0, + ITOL,RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,LICN, + U,ICN,IKEEP,W,NIND1,NIND2,NIND3,AM,JAM,IAM,IAMPT, + NZMAS,IPAR,RPAR,IERR) C IF (IERR .NE. 0) GOTO 8000 ITST = 3 ENDIF MEQC1 = MEQC1 + M1 + 1 C C NOW TEST TO SEE IF IT WAS SUCCESSFUL OR NOT C IF (.NOT.WORKED) THEN NFAIL = NFAIL + 1 C ********************************************************************** C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 4 TRIES. IF C PARTIALS ARE NOT UP TO DATE, THEY ARE RE-EVALUATED FOR THE C NEXT TRY. OTHERWISE THE Y ARRAY IS REPLACED BY ITS VALUES C BEFORE PREDICTION AND H IS REDUCED IF POSSIBLE. IF NOT A C NON-CONVERGENCE EXIT IS TAKEN C ********************************************************************** IF (IWEVAL.EQ.-1) THEN C HAVE BEEN USING OLD PARTIALS, UPDATE THEM AND TRY AGAIN IWEVAL = MITER CFAIL = .TRUE. CALL F(N,T,Y,SAVE2,IPAR,RPAR,IERR) IF (IERR .NE. 0) GOTO 8000 NFE = NFE + 1 GO TO 110 ENDIF IJUS=0 RED=0.5D+0 C *** failed at step 1 because of Newton GO TO 450 ENDIF IWEVAL = -1 HUSED = H NQUSED = NQ DO 140 I = 1,N Y(I,1) = (SAVE1(I)-ARH(I)) 140 CONTINUE C DO 141 I=1,N SAVE2(I) = 0.0d0 DO K= IAMPT(I),IAMPT(I+1)-1 SAVE2(I) = SAVE2(I) + AM(K)*Y(JAM(K),1) ENDDO 141 CONTINUE DO 142 I=1,N Y(I,1) = SAVE1(I) SAVE2(I) = SAVE2(I)*QQ 142 CONTINUE C C UPDATE THE DIFFERENCES AT N+1 C DO 160 J = 2,L JM1 = J-1 DO 150 I = 1,N Y(I,J) = Y(I,JM1) - YHOLD(I,JM1) 150 CONTINUE 160 CONTINUE C C COMPUTE ERROR IN THE SOLUTION C DO 161 I=1,N AYI = DABS(Y(I,1)) IF(ITOL.EQ.1) THEN SCALE(I) = YMAX(I) ELSE IF(ITOL.EQ.2) THEN SCALE(I) = RTOL(1)*AYI + ATOL(1) ELSE IF(ITOL.EQ.3) THEN SCALE(I) = RTOL(1)*AYI + ATOL(I) ELSE IF(ITOL.EQ.4) THEN SCALE(I) = RTOL(I)*AYI + ATOL(1) ELSE IF(ITOL.EQ.5) THEN SCALE(I) = RTOL(I)*AYI + ATOL(I) ENDIF 161 CONTINUE IF(NIND2.NE.0) THEN DO 162 I = NIND1+1,NIND2+NIND1 SCALE(I)=SCALE(I)/HUSED 162 CONTINUE ENDIF IF(NIND3.NE.0) THEN DO 163 I = NIND1 +NIND2 + 1,NIND1+NIND2+NIND3 SCALE(I)=SCALE(I)/(HUSED**2) 163 CONTINUE ENDIF C **** C **** AMMEND C **** CHANGE 1,N BELOW TO 1,NVARS C **** D = ZERO DO 170 I = 1,N D = D + ((Y(I,L)-YHOLD(I,L))/SCALE(I))**2 170 CONTINUE c c STORING Y FROM FIRST STEP FOR USE IN THIRD STEP. C IF(ITOL .EQ. 1) D = D/(RTOL(1)**2) IF(D.GT.E) GOTO 330 DO 180 I = 1,N YNHOLD(I,1) = Y(I,1) YNHOLD(I,2) = SAVE2(I) 180 CONTINUE KFAIL = 0 IREDO = 0 DO 190 I = 1,N ARH(I) = EL(2)*Y(I,1) 190 CONTINUE DO 210 J1 = 2,NQ JP1 = J1+1 DO 200 I = 1,N ARH(I) = ARH(I) + EL(JP1)*Y(I,J1) 200 CONTINUE 210 CONTINUE CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,F,IPAR,RPAR,IERR) IF (IERR .NE. 0 ) GOTO 8000 DO 220 I = 1,N SAVE1(I) = Y(I,1) ERROR(I) = ZERO 220 CONTINUE M2 = 0 C C FOR NOW WILL ASSUME THAT WE DO NOT WISH TO SAMPLE C AT THE N+2 STEP POINT C CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2, + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,1, + ITOL,RTOL,ATOL,HUSED,NBSOL,NFE,NQUSED,F,LICN, + U,ICN,IKEEP,W,NIND1,NIND2,NIND3,AM,JAM,IAM,IAMPT, + NZMAS,IPAR,RPAR,IERR) C IF (IERR .NE. 0) GOTO 8000 MEQC2 = MEQC2 + M2 + 1 C C NOW CHECK TO SEE IF IT WAS SUCCESSFUL OR NOT C IF (.NOT.WORKED) THEN NFAIL = NFAIL + 1 IJUS=0 RED=0.5D+0 C ***have failed on step 2 GOTO 450 ENDIF C C IF WE ARE DOWN TO HERE THEN THINGS MUST HAVE CONVERGED C LMP2=LMAX+2 LMP3=LMAX+3 DO 230 I = 1,N Y(I,LMP3) = (SAVE1(I)-ARH(I)) 230 CONTINUE C DO 231 I=1,N Y(I,LMP2) = 0.0D+0 DO K=IAMPT(I), IAMPT(I+1)-1 Y(I,LMP2) = Y(I,LMP2) + AM(K)*Y(JAM(K),LMP3) ENDDO 231 CONTINUE DO 232 I=1,N Y(I,LMP2) = Y(I,LMP2)*QQ Y(I,LMP3) = SAVE1(I) 232 CONTINUE C C WE ARE NOW COMPUTING THE THIRD STAGE C LL = L + 1 T = TOLD + H DELST = ELST(1)-EL(1) NQP2 = NQ+2 LMP2=LMAX+2 DO 280 I=1,N SCALE(I) = H*(ELST(NQP2)*Y(I,LMP2)+DELST*YNHOLD(I,2)) ARH(I) = SCALE(I) DO 270 J1 = 1,NQ ARH(I) = ARH(I) + ELST(J1+1)*YHOLD(I,J1) 270 CONTINUE 280 CONTINUE DO 290 I = 1,N SAVE2(I) = YNHOLD(I,2) Y(I,1) = YNHOLD(I,1) 290 CONTINUE M3STEP = 0 300 CONTINUE C DO 315 I=1,N SAVE1(I) = -Y(I,1) + ARH(I) - H*(ELST(NQP2)*Y(I,LMP2)+ + DELST*YNHOLD(I,2)) 315 CONTINUE DO 311 I=1,N SAVE2(I) = (SAVE2(I)*QI + H*(ELST(NQP2)*Y(I,LMP2)+ + DELST*YNHOLD(I,2))) DO K=IAMPT(I), IAMPT(I+1)-1 SAVE2(I) = SAVE2(I) + AM(K)*SAVE1(JAM(K)) ENDDO 311 CONTINUE DO 3111 I=1,N SAVE1(I) = SAVE2(I)/QQQ 3111 CONTINUE C IF ((MF .EQ. 25) .OR. (MF .EQ. 26)) THEN CALL MA28CD(N,PW,LICN,ICN,IKEEP,SAVE1,W,1) NBSOL = NBSOL + 1 ENDIF C DO 321 I=1,N AYI = DABS(Y(I,1)) IF(ITOL.EQ.1) THEN SCALE(I) = YMAX(I) ELSE IF(ITOL.EQ.2) THEN SCALE(I) = RTOL(1)*AYI + ATOL(1) ELSE IF(ITOL.EQ.3) THEN SCALE(I) = RTOL(1)*AYI + ATOL(I) ELSE IF(ITOL.EQ.4) THEN SCALE(I) = RTOL(I)*AYI + ATOL(1) ELSE IF(ITOL.EQ.5) THEN SCALE(I) = RTOL(I)*AYI + ATOL(I) ENDIF 321 CONTINUE IF(NIND2.NE.0) THEN DO 322 I = NIND1+1,NIND2+NIND1 SCALE(I)=SCALE(I)/HUSED 322 CONTINUE ENDIF IF(NIND3.NE.0) THEN DO 133 I = NIND1+NIND2 + 1,NIND1+NIND2+NIND3 SCALE(I)=SCALE(I)/(HUSED**2) 133 CONTINUE ENDIF D = ZERO DO 320 I = 1,N D = D + (SAVE1(I)/SCALE(I))**2 Y(I,1) = Y(I,1) + SAVE1(I) 320 CONTINUE IF(ITOL .EQ. 1) D = D/(RTOL(1)**2) IF ((D*DMIN1(ONE,2.0D+0*CRATE1)).LE.BND) GO TO 360 IF (M3STEP.EQ.5) THEN IJUS=1 RED=0.5D+0 C **** step 3 fails NFAIL = NFAIL + 1 GO TO 450 ENDIF M3STEP = M3STEP + 1 CALL F(N,T,Y,SAVE2,IPAR,RPAR,IERR) IF (IERR .NE. 0) GOTO 8000 NFE = NFE + 1 GO TO 300 330 KFAIL = KFAIL - 1 C ********************************************************************** C THE ERROR TEST FAILED. KFAIL KEEPS TRACK OF MULTIPLE FAILURES. C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES AND PREPARE TO C TRY THE STEP AGAIN. COMPUTE THE OPTIMAL STEP SIZE FOR THIS ORDER C AND ONE ORDER LOWER. C ********************************************************************** C *** failed on step 1 because of accuracy C COMPUTE ERROR IN THE SOLUTION C NFAIL = NFAIL + 1 DO 561 I=1,N AYI = DABS(Y(I,1)) IF(ITOL.EQ.1) THEN SCALE(I) = YMAX(I) ELSE IF(ITOL.EQ.2) THEN SCALE(I) = RTOL(1)*AYI + ATOL(1) ELSE IF(ITOL.EQ.3) THEN SCALE(I) = RTOL(1)*AYI + ATOL(I) ELSE IF(ITOL.EQ.4) THEN SCALE(I) = RTOL(I)*AYI + ATOL(1) ELSE IF(ITOL.EQ.5) THEN SCALE(I) = RTOL(I)*AYI + ATOL(I) ENDIF 561 CONTINUE IF(NIND2.NE.0) THEN DO 562 I = NIND1+1,NIND2+NIND1 SCALE(I)=SCALE(I)/HUSED 562 CONTINUE ENDIF IF(NIND3.NE.0) THEN DO 563 I = NIND1+NIND2 + 1,NIND1+NIND2+NIND3 SCALE(I)=SCALE(I)/(HUSED**2) 563 CONTINUE ENDIF DDOWN = ZERO TWODOWN = ZERO DO 1700 I=1,N DDOWN = DDOWN + ((Y(I,L))/SCALE(I))**2 TWODOWN = TWODOWN + ((Y(I,l-1))/SCALE(I))**2 1700 CONTINUE IF(ITOL .EQ. 1) D = D/(RTOL(1)**2) T = TOLD HOLD = H IF(NQ.GT.1) FFAIL = 0.5D+0/DBLE(FLOAT(NQ)) IF(NQ.GT.2) FRFAIL = 0.5D+0/DBLE(FLOAT(NQ-1)) EFAIL = 0.5D+0/DBLE(FLOAT(L)) CALL CPYARY(N*L,YHOLD,Y) RMAX = 2.0D+0 IF (DABS(H).LE.HMIN*1.00001D+0) THEN C C REQUESTED ERROR NOT POSSIBLE WITH GIVEN HMIN C KFLAG = -1 HOLD = H RETURN ENDIF IF (KFAIL.LE.-3) GO TO 340 IREDO = 2 C C PREDICTING A NEW H AFTER INSUFFICIENT ACCURACY C PRFAIL = ((D/(0.2D+0*E))**EFAIL)*1.5D+0 + 1.6D-6 PLFAIL = ((DDOWN/(0.2D+0*EDN))**FFAIL)*1.5D+0+1.7D-6 IF(NQ.GT.2) PLLFAIL =((TWODOWN/(0.2D+0*EDDN))**FRFAIL)* + 1.5D+0+1.7d-6 IF(PLLFAIL.GT.PLFAIL) PLFAIL=PLLFAIL IF(PLFAIL.LT.PRFAIL.AND.NQ.NE.1) THEN NEWQ=NQ-1 NQ=NEWQ RH=ONE/(PLFAIL*DBLE(FLOAT(-KFAIL))) L=NQ+1 CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) RC=RC*EL(1)/OLDLO OLDLO=EL(1) CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) ELSE NEWQ = NQ RH = ONE/ (PRFAIL*DBLE(FLOAT(-KFAIL))) ENDIF GO TO 40 C ********************************************************************** C CONTROL REACHES THIS STAGE IF 3 OR MORE FAILURES HAVE OCCURED. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE Y C ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST DERIVATIVE C IS RE-COMPUTED, AND THE ORDER IS SET TO 1. THEN H IS REDUCED BY A C FACTOR OF 10, AND THE STEP IS RETRIED. AFTER A TOTAL OF 7 C FAILURES AN EXIT IS TAKEN WITH KFLAG=-2. C ********************************************************************** 340 IF (KFAIL.EQ.-7) THEN C ERROR SMALLER THAN CAN BE HANDLED FOR PROBLEM KFLAG = -2 HOLD = H RETURN ENDIF C ********************************* C START FROM ORDER 1 AGAIN * C ********************************* JCHANG = 1 RH = DMAX1(HMIN/DABS(H),0.1D+0) CALL HCHOSE(RH,H,OVRIDE) H = H*RH CALL F(N,T,YHOLD,SAVE1,IPAR,RPAR,IERR) IF (IERR .NE. 0) GOTO 8000 NFE = NFE + 1 DO 350 I = 1,N Y(I,1) = YHOLD(I,1) Y(I,2) = H*SAVE1(I) YHOLD(I,2) = Y(I,2) 350 CONTINUE IWEVAL = MITER CFAIL = .TRUE. C SINCE WE HAVE HAD PROBLEMS PROCEED WITH THIS ORDER C FOR 10 STEPS (IF WE CAN) IDOUB = 10 IF (NQ.EQ.1) GO TO 60 NQ = 1 L = 2 C RESET ORDER, RECALCULATE ERROR BOUNDS CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) C NOW JUMP TO NORMAL CONTINUATION POINT GO TO 60 C ********************************************************************** C THE ITERATION FOR THE CORRECTED SOLUTION HAS CONVERGED. C UPDATE THE Y ARRAY. C ********************************************************************** 360 CONTINUE C **** C **** AMMEND **** C **** CHANGE 1,N BELOW TO 1,NVARS C **** DEMB=0.0D+0 DO 361 I=1,N DEMB=DEMB+((Y(I,1)-YNHOLD(I,1))/SCALE(I))**2 361 CONTINUE IF(DEMB.GT.4.0D+0*DBLE(FLOAT(N))) THEN IEMB=1 IJUS=1 RED=0.5D+0 C *** failed because of embedded error estimate NFAIL = NFAIL + 1 GOTO 450 ENDIF DO 380 J2 = 2,LL J2M1=J2-1 DO 370 I = 1,N Y(I,J2) = Y(I,J2M1) - YHOLD(I,J2M1) 370 CONTINUE 380 CONTINUE C --------------------------------------------------------------------- C IF THE COUNTER IDOUB EQUALS 2 AND WE ARE NOT ALREADY USING THE C MAXIMUM ALLOWABLE ORDER , STORE Y(I,LMAX+4) WHICH IS USED IN C ASSESSING THE POSSIBILITY OF INCREASING THE ORDER. IF IDOUB = 0 C CONTROL PASSES TO 480 WHERE AN ATTEMPT TO CHANGE THE STEPSIZE AND C ORDER IS MADE. C ---------------------------------------------------------------------- IF (IDOUB.EQ.2.AND.L.NE.LMAX) THEN LMP4=LMAX+4 DO 390 I = 1,N Y(I,LMP4) = Y(I,LL) 390 CONTINUE ENDIF IDOUB = IDOUB - 1 TRANGE=(TEND-TOLD-H)*H IF(TRANGE.LT.0.0D+0) THEN IDOUB = IDOUB + 2 GOTO 440 ENDIF JCHANG = 0 IF (IDOUB.EQ.0) THEN SAMPLE = .FALSE. ISAMP = ISAMP + 1 IF (ISAMP.EQ.4) THEN SAMPLE = .TRUE. ISAMP = 0 ENDIF C ********************************************************************** C NOW COMPUTE THE FACTORS PR1, PR2 AND PR3, BY WHICH C H COULD BE DIVIDED AT ORDER NQ-1, ORDER NQ AND ORDER NQ+1 C RESPECTIVELY. THE SMALLEST OF THESE IS DETERMINED AND THE NEW C ORDER CHOSEN ACCORDINGLY. IF THE ORDER IS TO BE INCREASED WE C MUST COMPUTE ONE MORE BACKWARD DIFFERENCE. C ********************************************************************** PR3 = 1.D+20 FAC = 1.5D+0 IF(IEMB.EQ.1) FAC = 1.8D+0 DO 400 I = 1,N AYI = DABS(Y(I,1)) IF(ITOL.EQ.1) THEN VHOLD = YMAX(I) ELSE IF(ITOL.EQ.2) THEN VHOLD = RTOL(1)*AYI + ATOL(1) ELSE IF(ITOL.EQ.3) THEN VHOLD = RTOL(1)*AYI + ATOL(I) ELSE IF(ITOL.EQ.4) THEN VHOLD = RTOL(I)*AYI + ATOL(1) ELSE IF(ITOL.EQ.5) THEN VHOLD = RTOL(I)*AYI + ATOL(I) ENDIF SCALE(I)=VHOLD 400 CONTINUE IF(NIND2.NE.0) THEN DO 4461 I=NIND1+1,NIND1+NIND2 SCALE(I) = SCALE(I)/HUSED 4461 CONTINUE ENDIF IF(NIND3.NE.0) THEN DO 4462 I=NIND1+NIND2+1,N SCALE(I) = SCALE(I)/(HUSED**2) 4462 CONTINUE ENDIF IF(L.NE.LMAX) THEN LMP4 = LMAX + 4 DUP = ZERO DO 401 I=1,N DUP = DUP + ((Y(I,LL)-Y(I,LMP4))/SCALE(I))**2 401 CONTINUE IF(ITOL .EQ. 1) DUP = DUP/(RTOL(1)**2) ENQ3 = 0.5D+0/DBLE(FLOAT(L+1)) PR3 = ((DUP/EUP)**ENQ3)*(FAC+0.2D+0) + 1.8D-6 ENDIF ENQ2 = 0.5D+0/DBLE(FLOAT(L)) D = ZERO c DDOWN = ZERO c DO 410 I = 1,N D = D + (Y(I,LL)/SCALE(I))**2 410 CONTINUE IF(ITOL .EQ.1) D = D/(RTOL(1)**2) PR2 = ((D/E)**ENQ2)*FAC + 1.6D-6 PR1 = 1.D+20 IF (NQ.GT.1) THEN DDDOWN=ZERO DDOWN = ZERO DO 57420 I = 1,N DDOWN = DDOWN + (Y(I,L)/SCALE(I))**2 DDDOWN = DDDOWN + (Y(I,L-1)/SCALE(I))**2 57420 CONTINUE IF(ITOL .EQ. 1) DDOWN = DDOWN/(RTOL(1)**2) ENQ1 = 0.5D+0/DBLE(FLOAT(NQ)) PR1 = ((DDOWN/EDN)**ENQ1)*(FAC+0.1D+0) + 1.7D-6 IF(NQ.GT.2) THEN ENQ0 = 0.5D+0/DBLE(FLOAT(NQ-1)) PR0 = ((DDDOWN/EDDN)**ENQ0)*(FAC+0.1D+0) + 1.7D-6 IF(PR0.GT.PR1) PR1 = PR0 IF(DDDOWN.LT.DDOWN) DDOWN = DDDOWN ENDIF ENDIF IF(L.EQ.LMAX) DUP = 0.0D+0 IF(NQ.LE.1) GOTO 6578 IF(DUP.GT.D.AND.D.GT.DDOWN) THEN PR2=1.0D+30 PR3=1.0D+30 ENDIF 6578 CONTINUE IF (PR2.LE.PR3) THEN IF (PR2.GT.PR1) THEN NEWQ = NQ - 1 RH = 1.0D+0/PR1 ELSE NEWQ = NQ RH = 1.0D+0/PR2 ENDIF ELSE IF (PR3.LT.PR1) THEN NEWQ = L RH = 1.0D+0/PR3 ELSE NEWQ = NQ - 1 RH = 1.0D+0/PR1 ENDIF IEMB=0 IF(RH.GT.1.0D+0.AND.RH.LT.1.1D+0) THEN IDOUB=10 NQ=NQUSED L=NQ+1 GOTO 440 ENDIF RH = DMIN1(RH,RMAX) CALL HCHOSE(RH,H,OVRIDE) IF ((JSINUP.LE.20).AND.(KFLAG.EQ.0).AND.(RH.LT.1.1D+0)) THEN C WE HAVE RUN INTO PROBLEMS IDOUB = 10 NQ = NQUSED L = NQ + 1 GO TO 440 ENDIF C ********************************************************************** C IF THERE IS A CHANGE IN ORDER, RESET NQ, L AND THE C COEFFICIENTS. IN ANY CASE H IS RESET AND THE C Y ARRAY IS RE-SCALED C ********************************************************************** IF(NIND3 .NE. 0) THEN IF(NQ.LE.2.AND.PR3.LT.1.0D+0) NEWQ=NQ+1 ENDIF IF (NEWQ.NE.NQ) THEN IF (NEWQ.GT.NQ) THEN C ADD AN EXTRA TERM TO THE HISTORY ARRAY DO 430 I = 1,N Y(I,LL) = Y(I,L) - YHOLD(I,L) 430 CONTINUE ENDIF NQ = NEWQ L = NQ + 1 C RESET ORDER,RECALCULATE ERROR BOUNDS CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) ENDIF RH = DMAX1(RH,HMIN/DABS(H)) RH = DMIN1(RH,HMAX/DABS(H),RMAX) CALL RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH RC = RC*RH IF(JSNOLD.GT.IBND) RC=ZERO IDOUB = L + 1 ENDIF 440 CONTINUE C ---------------------------------------------------------------------- C STORE THE Y ARRAY IN THE MATRIX YHOLD. STORE IN THE Y ARRAY THE C INFORMATION NECESSARY TO PERFORM AN INTERPOLATION TO FIND THE C SOLUTION AT THE SPECIFIED OUTPUT POINT IF APPROPRIATE. C ---------------------------------------------------------------------- CALL CPYARY(N*L,Y,YHOLD) NSTEP = NSTEP + 1 JSINUP = JSINUP + 1 JSNOLD = JSNOLD + 1 JSTART = NQUSED T = TOLD + HUSED HOLD = H KFAIL = 0 NEWPAR = 0 CFAIL = .FALSE. RETURN 450 CONTINUE FINISH = .FALSE. T=TOLD RMAX=2.0D+0 DO 460 J1=1,L DO 460 I=1,N Y(I,J1)=YHOLD(I,J1) 460 CONTINUE IF(DABS(H).LE.HMIN*1.00001D+0) THEN C C CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED C IF(NSTEP.EQ.0) THEN KFLAG=-1 ELSE KFLAG=-3 ENDIF C C TO SUPPRESS ERROR MESSAGES AT START AS H MAY C HAVE BEEN TOO LARGE ON THE FIRST STEP. C HOLD=H FINISH = .TRUE. ENDIF RH = RED IREDO=1 C C TRY AGAIN WITH UPDATED PARTIALS C 8000 IF (IERR .NE. 0) THEN write(*,1975) 1975 FORMAT (/,/,'IERR IS NON-ZERO BECAUSE OF AN ILLEGAL FUNCTION +CALL') h= h/2 IF(H.LT.EPSJAC/100.0D+0) THEN WRITE(6,9161) 9161 FORMAT (/,/,'STEPSIZE IS TOO SMALL') IDID = -7 RETURN ENDIF T= TOLD IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT TOUT WRITE (LOUT,*) T,TOUT,H CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) HO = H T0 = TOUT IDID = -5 RETURN ENDIF IERR = 0 jstart = -1 goto 30 endif C IF(IJUS.EQ.0) CALL HCHOSE(RH,H,OVRIDE) IF(.NOT.FINISH) THEN GO TO 40 ELSE RETURN ENDIF 9000 FORMAT (1X,' CORRECTOR HAS NOT CONVERGED') END C------------------- END OF SUBROUTINE STIFF -------------------------- SUBROUTINE RSCALE(N,L,RH,Y) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. INTEGER L,N C .. C .. ARRAY ARGUMENTS .. DIMENSION Y(N,12) C .. C .. LOCAL SCALARS .. INTEGER I,J,J1 C .. C .. LOCAL ARRAYS .. DIMENSION DI(8,8) C .. C .. DATA STATEMENTS .. C SUBROUTINE IS FOR RESCALING THE HISTORY ARRAY AFTER A CHANGE IN C STEPSIZE C C N ORDER OF THE PROBLEM C L NUMBER OF TERMS IN THE HISTORY ARRAY TO BE RESCALED C RH RATIO OF THE STEPSIZE CHANGE (I.E. RH = HNEW/HOLD) C Y() THE HISTORY ARRAY C DATA ZERO/0.0D+0/ C .. DI(2,2) = RH IF (L.GT.2) THEN TA = RH*RH DI(2,3) = RH* (1.0D+0-RH)/2.0D+0 DI(3,3) = TA IF (L.GT.3) THEN TB = TA*RH DI(2,4) = RH* ((RH-3.0D+0)*RH+2.0D+0)/6.0D+0 DI(3,4) = TA* (1.0D+0-RH) DI(4,4) = TB IF (L.GT.4) THEN TC = TB*RH DI(2,5) = - (((RH-6.0D+0)*RH+11.0D+0)*RH-6.0D+0)*RH/ + 24.0D+0 DI(3,5) = TA* ((7.0D+0*RH-18.0D+0)*RH+11.0D+0)/12.0D+0 DI(4,5) = 1.5D+0*TB* (1.0D+0-RH) DI(5,5) = TC IF (L.GT.5) THEN TD = TC*RH DI(2,6) = ((((RH-10.0D+0)*RH+35.0D+0)*RH-50.0D+0) + *RH+24.0D+0)*RH/120.0D+0 DI(3,6) = - (((3.0D+0*RH-14.0D+0)*RH+21.0D+0)*RH + -10.0D+0)*TA/12.0D+0 DI(4,6) = ((5.0D+0*RH-12.0D+0)*RH+7.0D+0)*TB/4.0D+0 DI(5,6) = 2.0D+0*TC* (1.0D+0-RH) DI(6,6) = TD IF (L.GT.6) THEN TE = TD*RH DI(2,7) = -RH* (RH-1.0D+0)* (RH-2.0D+0)* + (RH-3.0D+0)*(RH-4.0D+0)*(RH-5.0D+0)/ + 720.0D+0 DI(3,7) = TA* ((((62.0D+0*RH-450.0D+0)*RH+ + 1190.0D+0)*RH-1350.0D+0)*RH+548.0D+0) + /720.0D+0 DI(4,7) = TB* (((-18.0D+0*RH+75.0D+0)*RH + -102.0D+0)*RH+45.0D+0)/24.0D+0 DI(5,7) = TC* ((13.0D+0*RH-30.0D+0)*RH+17.0D+0) + /6.0D+0 DI(6,7) = 2.5D+0*TD* (1.0D+0-RH) DI(7,7) = TE IF (L.GT.7) THEN TF = TE*RH DI(2,8) = RH*(RH-1.0D+0)*(RH-2.0D+0)*(RH + -3.0D+0)*(RH-4.0D+0)*(RH-5.0D+0) + *(RH-6.0D+0)/5040.0D+0 DI(3,8) = TA* ((((((-126.0D+0*RH)+1302.0D+0)*RH- + 5250.0D+0)*RH+10290.0D+0)*RH-9744.0D+0 + )*RH+3528.0D+0)/5040.0D+0 DI(4,8) = TB* ((((43.0D+0*RH-270.0D+0)*RH+ + 625.0D+0)*RH-630.0D+0)*RH+232.0D+0) + /120.0D+0 DI(5,8) = TC* (((-10.0D+0*RH+39.0D+0)*RH- + 50.0D+0)*RH+21.0D+0)/6.0D+0 DI(6,8) = TD* ((20.0D+0*RH-45.0D+0)*RH+25.0D+0 + )/6.0D+0 DI(7,8) = 3.0D+0*TE* (1.0D+0-RH) DI(8,8) = TF ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF DO 30 I = 1,N DO 20 J = 2,L ZZ = ZERO DO 10 J1 = J,L ZZ = ZZ + DI(J,J1)*Y(I,J1) 10 CONTINUE Y(I,J) = ZZ 20 CONTINUE 30 CONTINUE RETURN END c---------------------------------------------------------------------- SUBROUTINE CPYARY(NELEM,SOURCE,TARGET) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C COPIES THE ARRAY SOURCE() INTO THE ARRAY TARGET() C C THIS SUBROUTINE COULD BE REPLACED BY THE BLAS ROUTINE SCOPY C (AFTER CHANGING THE ARGUMENT LIST APPROPRIATELY) C C .. SCALAR ARGUMENTS .. INTEGER NELEM C .. C .. ARRAY ARGUMENTS .. DIMENSION SOURCE(NELEM),TARGET(NELEM) C .. C .. LOCAL SCALARS .. INTEGER I C .. DO 10 I = 1,NELEM TARGET(I) = SOURCE(I) 10 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE HCHOSE(RH,H,OVRIDE) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON / STPSZE / HSTPSZ(2,14) LOGICAL OVRIDE C C FIRST MOVE ALL ELEMENTS DOWN ONE PLACE C IF (H.NE.HSTPSZ(2,1)) THEN DO 5 I=12,2,-1 I2=I-1 HSTPSZ(1,I)=HSTPSZ(1,I2) 5 HSTPSZ(2,I)=HSTPSZ(2,I2) C C NOW INSERT VALUE OF H USED BEFORE THIS CALL C HSTPSZ(1,2)=H/HSTPSZ(2,1) HSTPSZ(2,1)=H ENDIF C C NOW DECIDE ON THE NEW CHANGE C IF (RH.GT.1.0D+0) THEN OVRIDE=.FALSE. ELSE IF (HSTPSZ(1,2).LE.1.0D+0) THEN OVRIDE=.FALSE. ELSE IF ((RH*H).LE.HSTPSZ(2,2)) THEN OVRIDE=.FALSE. ELSE RH=HSTPSZ(2,2)/H OVRIDE=.TRUE. ENDIF HSTPSZ(1,1)=RH RETURN END C-------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. CHARACTER CMACH * .. * * Purpose * ======= * * DLAMCH determines double precision machine parameters. * * Arguments * ========= * * CMACH (input) CHARACTER*1 * Specifies the value to be returned by DLAMCH: * = 'E' or 'e', DLAMCH := eps * = 'S' or 's , DLAMCH := sfmin * = 'B' or 'b', DLAMCH := base * = 'P' or 'p', DLAMCH := eps*base * = 'N' or 'n', DLAMCH := t * = 'R' or 'r', DLAMCH := rnd * = 'M' or 'm', DLAMCH := emin * = 'U' or 'u', DLAMCH := rmin * = 'L' or 'l', DLAMCH := emax * = 'O' or 'o', DLAMCH := rmax * * where * * eps = relative machine precision * sfmin = safe minimum, such that 1/sfmin does not overflow * base = base of the machine * prec = eps*base * t = number of (base) digits in the mantissa * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise * emin = minimum exponent before (gradual) underflow * rmin = underflow threshold - base**(emin-1) * emax = largest exponent before overflow * rmax = overflow threshold - (base**emax)*(1-eps) * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, $ RND, SFMIN, SMALL, T * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLAMC2 * .. * .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, $ EMAX, RMAX, PREC * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) BASE = BETA T = IT IF( LRND ) THEN RND = ONE EPS = ( BASE**( 1-IT ) ) / 2 ELSE RND = ZERO EPS = BASE**( 1-IT ) END IF PREC = EPS*BASE EMIN = IMIN EMAX = IMAX SFMIN = RMIN SMALL = ONE / RMAX IF( SMALL.GE.SFMIN ) THEN * * Use SMALL plus a bit, to avoid the possibility of rounding * causing overflow when computing 1/sfmin. * SFMIN = SMALL*( ONE+EPS ) END IF END IF * IF( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF * DLAMCH = RMACH RETURN * * End of DLAMCH * END * ************************************************************************ * SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * DLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * The routine is based on the routine ENVRON by Malcolm and * incorporates suggestions by Gentleman and Marovich. See * * Malcolm M. A. (1972) Algorithms to reveal properties of * floating-point arithmetic. Comms. of the ACM, 15, 949-951. * * Gentleman W. M. and Marovich S. B. (1974) More on algorithms * that reveal properties of floating point arithmetic units. * Comms. of the ACM, 17, 276-277. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * Compute a = 2.0**m with the smallest positive integer m such * that * * fl( a + 1.0 ) = a. * A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 10 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 10 END IF *+ END WHILE * * Now compute b = 2.0**m with the smallest positive integer m * such that * * fl( a + b ) .gt. a. * B = 1 C = DLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = DLAMC3( A, B ) GO TO 20 END IF *+ END WHILE * * Now compute the base. a and c are neighbouring floating point * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so * their difference is beta. Adding 0.25 to c is to ensure that it * is truncated to beta and not ( beta - 1 ). * QTR = ONE / 4 SAVEC = C C = DLAMC3( C, -A ) LBETA = C + QTR * * Now determine whether rounding or chopping occurs, by adding a * bit less than beta/2 and a bit more than beta/2 to a. * B = LBETA F = DLAMC3( B / 2, -B / 100 ) C = DLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = DLAMC3( B / 2, B / 100 ) C = DLAMC3( F, A ) IF( ( LRND ) .AND. ( C.EQ.A ) ) $ LRND = .FALSE. * * Try and decide whether rounding is done in the IEEE 'round to * nearest' style. B/2 is half a unit in the last place of the two * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit * zero, and SAVEC is odd. Thus adding B/2 to A should not change * A, but adding B/2 to SAVEC should change SAVEC. * T1 = DLAMC3( B / 2, A ) T2 = DLAMC3( B / 2, SAVEC ) LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND * * Now find the mantissa, t. It should be the integer part of * log to the base beta of a, however it is safer to determine t * by powering. So we find t as the smallest positive integer for * which * * fl( beta**t + 1.0 ) = 1.0. * LT = 0 A = 1 C = 1 * *+ WHILE( C.EQ.ONE )LOOP 30 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 RETURN * * End of DLAMC1 * END * ************************************************************************ * SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T DOUBLE PRECISION EPS, RMAX, RMIN * .. * * Purpose * ======= * * DLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies whether proper rounding ( RND = .TRUE. ) or * chopping ( RND = .FALSE. ) occurs in addition. This may not * be a reliable guide to the way in which the machine performs * its arithmetic. * * EPS (output) DOUBLE PRECISION * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) DOUBLE PRECISION * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) DOUBLE PRECISION * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. External Subroutines .. EXTERNAL DLAMC1, DLAMC4, DLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN FIRST = .FALSE. ZERO = 0 ONE = 1 TWO = 2 * * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of * BETA, T, RND, EPS, EMIN and RMIN. * * Throughout this routine we use the function DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) * * Start to find EPS. * B = LBETA A = B**( -LT ) LEPS = A * * Try some tricks to see whether or not this is the correct EPS. * B = TWO / 3 HALF = ONE / 2 SIXTH = DLAMC3( B, -HALF ) THIRD = DLAMC3( SIXTH, SIXTH ) B = DLAMC3( THIRD, -HALF ) B = DLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = DLAMC3( HALF, -C ) B = DLAMC3( HALF, C ) C = DLAMC3( HALF, -B ) B = DLAMC3( HALF, C ) GO TO 10 END IF *+ END WHILE * IF( A.LT.LEPS ) $ LEPS = A * * Computation of EPS complete. * * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). * Keep dividing A by BETA until (gradual) underflow occurs. This * is detected when we cannot recover the previous A. * RBASE = ONE / LBETA SMALL = ONE DO 20 I = 1, 3 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = DLAMC3( ONE, SMALL ) CALL DLAMC4( NGPMIN, ONE, LBETA ) CALL DLAMC4( NGNMIN, -ONE, LBETA ) CALL DLAMC4( GPMIN, A, LBETA ) CALL DLAMC4( GNMIN, -A, LBETA ) IEEE = .FALSE. * IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN IF( NGPMIN.EQ.GPMIN ) THEN LEMIN = NGPMIN * ( Non twos-complement machines, no gradual underflow; * e.g., VAX ) ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN LEMIN = NGPMIN - 1 + LT IEEE = .TRUE. * ( Non twos-complement machines, with gradual underflow; * e.g., IEEE standard followers ) ELSE LEMIN = MIN( NGPMIN, GPMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) * ( Twos-complement machines, no gradual underflow; * e.g., CYBER 205 ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. $ ( GPMIN.EQ.GNMIN ) ) THEN IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT * ( Twos-complement machines with gradual underflow; * no known machine ) ELSE LEMIN = MIN( NGPMIN, NGNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF * ELSE LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) * ( A guess; no known machine ) IWARN = .TRUE. END IF *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )LEMIN END IF *** * * Assume IEEE arithmetic if we found denormalised numbers above, * or if arithmetic seems to round in the IEEE style, determined * in routine DLAMC1. A true IEEE machine should have both things * true; however, faulty machines may have one or the other. * IEEE = IEEE .OR. LIEEE1 * * Compute RMIN by successive division by BETA. We could compute * RMIN as BASE**( EMIN - 1 ), but some machines underflow during * this computation. * LRMIN = 1 DO 30 I = 1, 1 - LEMIN LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 30 CONTINUE * * Finally, call DLAMC5 to compute EMAX and RMAX. * CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) END IF * BETA = LBETA T = LT RND = LRND EPS = LEPS EMIN = LEMIN RMIN = LRMIN EMAX = LEMAX RMAX = LRMAX * RETURN * 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', $ ' EMIN = ', I8, / $ ' If, after inspection, the value EMIN looks', $ ' acceptable please comment out ', $ / ' the IF block as marked within the code of routine', $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) * * End of DLAMC2 * END * ************************************************************************ * DOUBLE PRECISION FUNCTION DLAMC3( A, B ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. DOUBLE PRECISION A, B * .. * * Purpose * ======= * * DLAMC3 is intended to force A and B to be stored prior to doing * the addition of A and B , for use in situations where optimizers * might hold one of these in a register. * * Arguments * ========= * * A, B (input) DOUBLE PRECISION * The values A and B. * * ===================================================================== * * .. Executable Statements .. * DLAMC3 = A + B * RETURN * * End of DLAMC3 * END * ************************************************************************ * SUBROUTINE DLAMC4( EMIN, START, BASE ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER BASE, EMIN DOUBLE PRECISION START * .. * * Purpose * ======= * * DLAMC4 is a service routine for DLAMC2. * * Arguments * ========= * * EMIN (output) EMIN * The minimum exponent before (gradual) underflow, computed by * setting A = START and dividing by BASE until the previous A * can not be recovered. * * START (input) DOUBLE PRECISION * The starting point for determining EMIN. * * BASE (input) INTEGER * The base of the machine. * * ===================================================================== * * .. Local Scalars .. INTEGER I DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Executable Statements .. * A = START ONE = 1 RBASE = ONE / BASE ZERO = 0 EMIN = 1 B1 = DLAMC3( A*RBASE, ZERO ) C1 = A C2 = A D1 = A D2 = A *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 10 CONTINUE IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. $ ( D2.EQ.A ) ) THEN EMIN = EMIN - 1 A = B1 B1 = DLAMC3( A / BASE, ZERO ) C1 = DLAMC3( B1*BASE, ZERO ) D1 = ZERO DO 20 I = 1, BASE D1 = D1 + B1 20 CONTINUE B2 = DLAMC3( A*RBASE, ZERO ) C2 = DLAMC3( B2 / RBASE, ZERO ) D2 = ZERO DO 30 I = 1, BASE D2 = D2 + B2 30 CONTINUE GO TO 10 END IF *+ END WHILE * RETURN * * End of DLAMC4 * END * ************************************************************************ * SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. LOGICAL IEEE INTEGER BETA, EMAX, EMIN, P DOUBLE PRECISION RMAX * .. * * Purpose * ======= * * DLAMC5 attempts to compute RMAX, the largest machine floating-point * number, without overflow. It assumes that EMAX + abs(EMIN) sum * approximately to a power of 2. It will fail on machines where this * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, * EMAX = 28718). It will also fail if the value supplied for EMIN is * too large (i.e. too close to zero), probably with overflow. * * Arguments * ========= * * BETA (input) INTEGER * The base of floating-point arithmetic. * * P (input) INTEGER * The number of base BETA digits in the mantissa of a * floating-point value. * * EMIN (input) INTEGER * The minimum exponent before (gradual) underflow. * * IEEE (input) LOGICAL * A logical flag specifying whether or not the arithmetic * system is thought to comply with the IEEE standard. * * EMAX (output) INTEGER * The largest exponent before overflow * * RMAX (output) DOUBLE PRECISION * The largest machine floating-point number. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP DOUBLE PRECISION OLDY, RECBAS, Y, Z * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. * .. Executable Statements .. * * First compute LEXP and UEXP, two powers of 2 that bound * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum * approximately to the bound that is closest to abs(EMIN). * (EMAX is the exponent of the required number RMAX). * LEXP = 1 EXBITS = 1 10 CONTINUE TRY = LEXP*2 IF( TRY.LE.( -EMIN ) ) THEN LEXP = TRY EXBITS = EXBITS + 1 GO TO 10 END IF IF( LEXP.EQ.-EMIN ) THEN UEXP = LEXP ELSE UEXP = TRY EXBITS = EXBITS + 1 END IF * * Now -LEXP is less than or equal to EMIN, and -UEXP is greater * than or equal to EMIN. EXBITS is the number of bits needed to * store the exponent. * IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN EXPSUM = 2*LEXP ELSE EXPSUM = 2*UEXP END IF * * EXPSUM is the exponent range, approximately equal to * EMAX - EMIN + 1 . * EMAX = EXPSUM + EMIN - 1 NBITS = 1 + EXBITS + P * * NBITS is the total number of bits needed to store a * floating-point number. * IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN * * Either there are an odd number of bits used to store a * floating-point number, which is unlikely, or some bits are * not used in the representation of numbers, which is possible, * (e.g. Cray machines) or the mantissa has an implicit bit, * (e.g. IEEE machines, Dec Vax machines), which is perhaps the * most likely. We have to assume the last alternative. * If this is true, then we need to reduce EMAX by one because * there must be some way of representing zero in an implicit-bit * system. On machines like Cray, we are reducing EMAX by one * unnecessarily. * EMAX = EMAX - 1 END IF * IF( IEEE ) THEN * * Assume we are on an IEEE machine which reserves one exponent * for infinity and NaN. * EMAX = EMAX - 1 END IF * * Now create RMAX, the largest machine number, which should * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . * * First compute 1.0 - BETA**(-P), being careful that the * result is less than 1.0 . * RECBAS = ONE / BETA Z = BETA - ONE Y = ZERO DO 20 I = 1, P Z = Z*RECBAS IF( Y.LT.ONE ) $ OLDY = Y Y = DLAMC3( Y, Z ) 20 CONTINUE IF( Y.GE.ONE ) $ Y = OLDY * * Now multiply by BETA**EMAX to get RMAX. * DO 30 I = 1, EMAX Y = DLAMC3( Y*BETA, ZERO ) 30 CONTINUE * RMAX = Y RETURN * * End of DLAMC5 * END LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END ****************************************************************** C######DATE 01 JAN 1984 COPYRIGHT UKAEA, HARWELL. C######ALIAS TD02AD C D AND S ARE STANDARD FORTRAN DOUBLE AND SINGLE VERSIONS. C... SUBROUTINE TD02AD(M,N,IRN,IP,F,H,X,Y,FF,HMAX,A,IG,W,Z, + IPAR,RPAR,IERR) C... COMMON/TD02DD/UMIN,UAIM,UMAX,EPS,EPS1,LP,ADJUST INTEGER IRN(*),IP(N+1),IG(2*N+1),IPAR(2) DOUBLE PRECISION HMAX(N),H(N),Z(*),TRUNC,ROUND,HM DOUBLE PRECISION X(N),FF(M),A(*),W(*),Y,RPAR(2), & DER,UAIM,EPS,UMAX,HJL,ONE,UMIN,XJ,EPS1,ZERO C EXTERNAL F LOGICAL REPEAT,ADJUST DATA ZERO/0.D0/,ONE/1.D0/ C C IG(J) POINTS TO THE FIRST ENTRY OF GROUP J NFUN=0 N1=N+1 MJ=M+J C... IG(N1+J),J=1,2,... ARE COLUMN NUMBERS FOR GROUP 1 THEN GROUP 2 ETC IF(IG(1).NE.0)GO TO 70 ICC=1 DO 10 J=1,N IG(J+1)=0 MJ=M+J 10 W(MJ )=ZERO C... W(M+J) IS ZERO IF COLUMN J HAS NOT BEEN INCLUDED IN A GROUP YET DO 60 NC=1,N1 IG(NC)=ICC DO 20 I=1,M 20 W(I)=ZERO C... W(I),I=1,M HOLDS THE BOOLEAN PATTERN OF THE UNION OF THE COLUMNS C... SO FAR INCLUDED IN GROUP NC DO 50 J=1,N MJ=M+J IF(W(MJ ).NE.ZERO)GO TO 50 K1=IP(J) K2=IP(J+1)-1 IF(K2.LT.K1)GO TO 50 DO 30 K=K1,K2 I=IRN(K) IF(W(I).NE.ZERO)GO TO 50 30 CONTINUE C... ACCEPT COLUMN DO 40 K=K1,K2 I=IRN(K) 40 W(I)=ONE I= N1+ICC IG(I)=J ICC=ICC+1 W(MJ)=ONE 50 CONTINUE IF(ICC.EQ.IG(NC))GO TO 70 60 CONTINUE C... PRESERVE X IN W(M+1),...,W(M+N). CHECK ALL H(J). 70 DO 82 J=1,N HM=-HMAX(1) IF(HM.LT.0.)HM=HMAX(J) XJ=DABS(X(J)) IF(ADJUST) 1 H(J)=(XJ+DMAX1(dsqrt(EPS)*XJ,DBLE(DMIN1(H(J),HM)), 2 EPS1*HM))-XJ MJ=M+J 82 W(MJ )=X(J) C... Z(J) HOLDS MAXIMUM OF ESTIMATED RATIO OF TRUNCATION ERROR TO C... ROUNDOFF ERROR IN COLUMN J,J=1,2,...,N. IP(J) IS NEGATED IF C... H(J) HAS DECREASED DURING THIS CALL OF TD02. DO 163 NIMP=1,60 DO 85 J=1,N 85 Z(J)=1. C... FIND INITIAL APPROXIMATE DERIVATIVES DO 105 NG=1,N L1=N1+IG(NG) L2=N1+IG(NG+1)-1 IF(L2.LT.L1)GO TO 110 DO 90 L=L1,L2 J=IG(L) IF(H(J))90,90,93 90 CONTINUE GO TO 105 93 DO 95 L=L1,L2 J=IG(L) H(J)=ABS(H(J)) 95 X(J)=X(J)+H(J) CALL F(N,Y,X,W,IPAR,RPAR,IERR) NFUN=NFUN+1 DO 100 L=L1,L2 J=IG(L) MJ=M+J X(J)=W(MJ) K1=IABS(IP(J)+0) K2=IABS(IP(J+1)+0)-1 IF(K2.LT.K1)GO TO 100 DO 98 K=K1,K2 I=IRN(K) A(K)=(W(I)-FF(I))/H(J) 98 continue 100 CONTINUE 105 CONTINUE 110 IF(.NOT.ADJUST)GO TO 180 C... ESTIMATE ERRORS AND IMPROVED STEP-LENGTHS DO 133 NG=1,N L1=N1+IG(NG) L2=N1+IG(NG+1)-1 IF(L2.LT.L1)GO TO 135 DO 113 L=L1,L2 J=IG(L) IF(H(J))113,113,117 113 CONTINUE GO TO 133 117 DO 120 L=L1,L2 J=IG(L) 120 X(J)=X(J)-H(J) CALL F(N,Y,X,W,IPAR,RPAR,IERR) NFUN=NFUN+1 DO 130 L=L1,L2 J=IG(L) MJ=M+J X(J)=W(MJ) K1=IABS(IP(J)+0) K2=IABS(IP(J+1)+0)-1 IF(K2.LT.K1)GO TO 130 DO 123 K=K1,K2 I=IRN(K) DER=(FF(I)-W(I))/H(J) ROUND=EPS*(0.5*(DABS(W(I))+DABS(A(K)*H(J)+FF(I))) & +DMAX1(DABS(A(K)),DABS(DER))*(DABS(X(J))+ & H(J)))/H(J) A(K)=(A(K)+DER)/2. TRUNC=DER-A(K) IF(ROUND.EQ.0.)ROUND=1. Z(J)=DMAX1(Z(J),ABS(TRUNC/ROUND)) 123 CONTINUE 130 CONTINUE 133 CONTINUE C... FIND IMPROVED STEPS AND DECIDE WHETHER FURTHER SWEEP IS NEEDED 135 REPEAT=.FALSE. DO 160 J=1,N IF(H(J).LT.0.)GO TO 160 HJL=H(J) XJ=DABS(X(J)) HM=-HMAX(1) IF(HM.LT.0.)HM=HMAX(J) IF(Z(J).LT.UMIN)GO TO 140 H(J)=(DMAX1(H(J)*DSQRT(UAIM/Z( J)),EPS*XJ,EPS1*HM)+XJ)-XJ IF(H(J).LT.HJL)IP(J)=-IABS(IP(J)+0) IF(Z(J).GT.UMAX)GO TO 150 138 H(J)=-H(J) GO TO 160 140 H(J)=(DMIN1(DSQRT(UAIM/Z( J))*H(J), DBLE(HM))+XJ)-XJ 150 IF(DABS(HJL/H(J)-1.).LE.0.001 )GO TO 138 IF(NIMP.GE.6 .AND. IP(J).LT.0)GO TO 138 REPEAT=.TRUE. 160 CONTINUE IF(.NOT.REPEAT)GO TO 166 163 CONTINUE 166 DO 170 J=1,N IP(J)=IABS(IP(J)+0) 170 H(J)= ABS(H(J)) stop 180 W(1)=NFUN ADJUST = .FALSE. I=2 W(I)=NC-1 RETURN END C---------------------------------------------------------------------- C... SUBROUTINE TD02BD(M,N,IRN,IP,F,H,X,T,FF,W,IA, + SCALE,IPAR,RPAR,IERR) C... COMMON/TD02DD/UMIN,UAIM,UMAX,EPS,EPS1,LP,ADJUST INTEGER IRN(*),IP(N+1),IPAR(*) DOUBLE PRECISION H(N),DXJ,SCALE(N),RPAR(*) DOUBLE PRECISION X(N),FF(M),W(2*N),T,XJ,UMIN,UAIM,UMAX,EPS,EPS1 C EXTERNAL F LOGICAL ADJUST C... K=1 IP(1)=1 DO 200 J=1,N XJ=X(J) DXJ = DSIGN(SCALE(J),XJ) X(J)=XJ+DXJ CALL F(N,T,X,W,IPAR,RPAR,IERR) X(J)=XJ DO 210 I=1,M IF(DABS((W(I)-FF(I))/DXJ).LE.0.d0) GO TO 210 IF(K.GT.IA) THEN WRITE(LP,*) 'NZMAX SHOULD BE AT LEAST',K GO TO 223 ENDIF IRN(K)=I K=K+1 210 CONTINUE IP(1+J) = K 200 continue W(1) = J -1 RETURN 223 WRITE(LP,226) 226 FORMAT(49H0ERROR RETURN FROM TD02BD BECAUSE IA IS TOO SMALL) IP(1)=-1 RETURN END C------------------------------------------------------------------- C.... BLOCK DATA TD02ED COMMON/TD02DD/UMIN,UAIM,UMAX,EPS,EPS1,LP,ADJUST DOUBLE PRECISION UMIN,UAIM,UMAX,EPS,EPS1 LOGICAL ADJUST DATA UMIN/10.D0/, UAIM/1.D2/, UMAX/1.D3 /, EPS/25.D-17/, & LP/6/,EPS1/25.D-17/,ADJUST/.FALSE./ END C... C----------------------------------------------------------------------- C... SUBROUTINE SPARSDET(N,T,Y0,HO,SCALE,HDUMMY,YDOT,IP,NZMAX,IRNT, & F,NZ,JCOL,IROW,W1,IPAR,RPAR,IERR) C... IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION Y0(N),HO,SCALE(N),HDUMMY(N),YDOT(N) DOUBLE PRECISION W1(2*N),RPAR(*),FAC INTEGER IP(N+1),JCOL(*),IROW(*),IPAR(*),IRNT(NZMAX) C... C... EXTERNAL SUBROUTINES EXTERNAL F,TD02BD C... C... INTRINSIC FUNCTIONS INTRINSIC DABS,DMAX1 C... C... PERTURB INITIAL CONDITION, Y, FOR STRUCTURE DETERMINATION C... DO I= 1,N FAC = 1.0D0+1.0D0/(DBLE(I)+1.0D0) Y0(I) = Y0(I) + FAC*DSIGN(SCALE(I),Y0(I)) ENDDO C... C... EVALUATION OF THE SPARSITY PATTERN C... DO K = 1,N HDUMMY(K) = K*HO ENDDO C... C... COMPUTE THE SPARSITY STRUCTURE FROM RESULT OF N+1 CALL TO F C... CALL F(N,T,Y0,YDOT,IPAR,RAPR,IERR) CALL TD02BD(N,N,IRNT,IP,F,HDUMMY,Y0,T,YDOT,W1,NZMAX,SCALE, + IPAR,RPAR,IERR) C... NZ = IP(N+1) -1 DO I=1,NZ IROW(I) = IRNT(I) ENDDO C DO K =1,N K1 =IP(K) K2 =IP(K+1) -1 DO I=K1,K2 JCOL(I) = K ENDDO ENDDO RETURN END C***********************************************************************