c 19/2/2000 c----------------------------------------------------------------------------- c----------------------------------------------------------------------------- c c MEBDF VERSION THAT USES YSMP (YALE SPARSE MATRIX PACKAGE) c c----------------------------------------------------------------------------- c----------------------------------------------------------------------------- SUBROUTINE MEBDFSO(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT, + LWORK, WORK,LIWORK,IWORK,MAXDER,ITOL,RTOL,ATOL,F,PDERV, + NSP,NONZ,IPAR,RPAR,IERR) c **************************************************************************** c **************************************************************************** c Written by T.J. Abdulla and J.R. Cash, 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, + JK1,JK2,JK3,JK4,JK5,JK6,JK7,JK8 C LOGICAL ORDER C COMMON BLOCKS COMMON /ORDERING/ORDER C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL OVDRIV,F,PDERV C .. C .. SAVE STATEMENT .. SAVE I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13, + JK1,JK2,JK3,JK4,JK5,JK6,JK7,JK8 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 + NONZ I12 = I11 + NONZ I13 = I12 + NSP C C.....INTEGER WORKSPACE C JK1 = N+15 JK2 = JK1 + NONZ JK3 = JK2 + N JK4 = JK3 + N JK5 = JK4 + NONZ JK6 = JK5 + N JK7 = JK6 + N JK8 = JK7 + N UROUND=DLAMCH('Epsilon') WORK(1)=UROUND EPSJAC=SQRT(WORK(1)) C IF (LWORK.LT.(I13-1)) THEN IDID = -11 WRITE (LOUT,9000) I13 - 1 ENDIF IF (LIWORK.LT.JK8-1) THEN IDID = -12 WRITE (LOUT,9010) JK8 - 1 ENDIF ENDIF IF (IDID.LT.0) RETURN ORDER = .FALSE. ENDIF C C THE DIMENSION OF THE REAL WORKSPACE, WORK, HAS TO BE AT LEAST C (33*N + 2*NONZ + NSP + 3). WHILE THE DIMENSION OF THE INTEGER c WORKSPACE HAS TO BE AT LEAST (6*N + 2*NONZ + 15), WHERE NONZ IS C THE NUMBER OF THE NONZERO ENTRIES IN THE JACOBIAN MATRIX (IF THIS C IS UNKNOWN, USE AN ESTIMATE) AND NSP SHOULD BE LARGER THAN C 8*N+2+2*NONZ (WORKSPACE FOR THE SPARSE SOLVER) 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),WORK(I12),IWORK(13), + IWORK(JK1),IWORK(JK2),IWORK(JK3),WORK(I12),MAXDER,ITOL, + RTOL,ATOL,F,PDERV,IWORK(1),IWORK(2),IWORK(3),IWORK(4), + IWORK(5),IWORK(6),IWORK(7),IWORK(8),IWORK(9),IWORK(10), + IWORK(11),IWORK(12),UROUND,WORK(2),EPSJAC,NSP,NONZ, + IWORK(JK4),IWORK(JK5),IWORK(JK6),IWORK(JK7),IPAR,RPAR,IERR) 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(NONZ) C PWCOPY(NONZ) 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,SAVF,PW,PWCOPY, + RSP,IA,JA,P,IP,ISP,MAXDER,ITOL,RTOL,ATOL,F,PDERV,NQUSED, + NSTEP,NFAIL,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,MAXSTP, + NNZ,UROUND,HUSED,EPSJAC,NSP,NONZ,IGP,JGP,INCL,JDONE, + IPAR,RPAR,IERR) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C THIS IS THE FEBRUARY 23rd 2000 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 6. S.C. EISENSTAT, M.C. GURSKY, M.H. SCHULTZ, AND, A.H. SHERMAN, C YALE SPARSE MATRIX PACKAGE II, THE NONSYEMMETRIC CODES, RESEARCH C REPORT NO. 114, DEPT. OF COMPUTER SCIENCES, YALE UNIVERSITY, 1977. C----------------------------------------------------------------------------- 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 OR MF = 26 ARE C ALLOWED. THESE ARE EXTENDED BACKWARD DIFFERENTIATION C FORMULAE USING THE CHORD METHOD WITH ANALYTIC JACOBIAN C FOR MF = 25 AND NUMERICAL JACOBIAN FOR MF = 26. C WITH MF = 25 THE USER NEED TO SPECIFY SUBROUTINE PDERV. 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 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 NONZ THE NUMBER OF NONZERO ELEMENTS IN THE JACOBIAN MATRIX. C IF THIS UNKNOWN, USE AN ESTIMATE. 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 SPARSDET( - ) CONSTRUCT THE SPARSITY STRUCUTRE OF THE JACOBIAN MATRIX C AND COLUMN GROUPING FOR THE NUMERICAL JACOBIAN. C C AND THE FOLLOWING YALE SPARSE ROUTINES C OVDRV( - ) CONSTRUCT A REORDERING OF THE ROWS AND COLUMNS OF C THE ITERATION MATRIX BY MINIMUM DEGREE ALGORITHM C CDRV( - ) PERFORMS REORDERING, SYMBOLIC FACTORIZATION, NUMERICAL C FACTORIZATION OR LINEAR SYSTEM SOLUTION OPERATIONS C DEPENDING ON A PATH ARGUMENT. 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,J,PD,IPAR,RPAR,IERR) C COMPUTES THE NONZERO ENTRIES IN THE Jth COLUMN OF THE C JACOBIAN MATRIX AND STORES IT IN PD, I.E., C PD(J) = DF(I)/DY(J) FOR ALL RELEVANT VALUES OF I. C PDERV IS CALLED ONLY IF MITER = 5. OTHERWISE A DUMMY C ROUTINE CAN BE SUBSTITUTED. C C C THE DIMENSIONS OF YMAX,ERROR,SAVE1,SAVE2 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(1)) LAST ORDER SUCCESSFULLY USED C NSTEP (=IWORK(2)) NUMBER OF SUCCESSFUL STEPS TAKEN SO FAR C NFAIL (=IWORK(3)) NUMBER OF FAILED STEPS C NFE (=IWORK(4)) NUMBER OF FUNCTION EVALUATIONS SO FAR C NJE (=IWORK(5)) NUMBER OF JACOBIAN EVALUATIONS SO FAR C NDEC (=IWORK(6)) NUMBER OF LU DECOMPOSITIONS SO FAR C NBSOL (=IWORK(7)) NUMBER OF 'BACKSOLVES' SO FAR C NPSET (=IWORK(8)) NUMBER OF TIMES A NEW COEFFICIENT MATRIX HAS BEEN C FORMED SO FAR C NCOSET (=IWORK(9)) NUMBER OF TIMES THE ORDER OF THE METHOD USED HAS C BEEN CHANGED SO FAR C MAXORD (=IWORK(10)) THE MAXIMUM ORDER USED SO FAR IN THE INTEGRATION C MAXSTP (=IWORK(11)) THE MAXIMUM ALLOWED NUMBER OF STEPS SET BY THE USER C NNZ (=IWORK(12)) NUMBER OF NONZER ELEMENTS IN THE JACOBIAN MATRIX 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),RSP(NSP),IGP(NONZ),JGP(N), + INCL(N),JDONE(N),IPAR(*),RPAR(*) C .. C .. LOCAL SCALARS .. INTEGER I,KGO,NHCUT,NSP,IA(N+1),JA(NONZ),IP(N),P(N),ISP(NSP) C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL INTERP,STIFF,F,PDERV,SPARSDET C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1 C .. C .. COMMON BLOCKS .. INTEGER IDOUB,ISAMP,JSTART,KFLAG,L,MAXORD,NBSOL,NCOSET,NDEC, + NEWPAR,NFE,NJE,NPSET,NQUSED,NRENEW,NSTEP C SAVE T,H,HMIN,HMAX,KFLAG,JSTART LOGICAL ORDER COMMON/ORDERING/ORDER 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(LOUT,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 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 ENDIF ENDIF C C THE NEXT ROUTINE CONSTRUCT THE SPARSITY STRUCTURE OF THE ITERATION C MATRIX (ROW, COLUMN) INDICES. THIS NEED TO BE DONE ONLY ONCE C PER INTEGRATION C IF (.NOT. ORDER ) 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 CALL SPARSDET(N,T,Y0,PWCOPY,SCALE,IA,JA,SAVF,P,IP,NSP, + ISP,RSP,NONZ,NGP,JGP,IGP,MF,INCL,JDONE,LOUT,F,NNZ, + IPAR,RPAR,IERR) ENDIF 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,NQUSED,NSTEP, + NFAIL,NFE,NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,UROUND, + EPSJAC,HUSED,RSP,IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGP,JGP, + IGP,LOUT,IPAR,RPAR,IERR) C IF (IERR .NE. 0)THEN WRITE(LOUT,9150) RETURN 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 C 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',/) 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',/,/) 9150 FORMAT (/,/,'IERR IS NON-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,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP 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, + NRENEW,YMAX,SAVE1,SAVE2,PW,PWCOPY,WRKSPC,ITOL,RTOL,ATOL, + NPSET,NJE,NFE,NDEC,RSP,IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGRP, + JGP,IGP,LOUT,IPAR,RPAR,IERR) 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), & YMAX(N),SAVE1(N),RTOL(1),ATOL(1),RSP(NSP),FTEM(1),SAVF(N), & IPAR(*),RPAR(*) INTEGER P(N),IP(N),IA(N+1),JA(NONZ),ISP(NSP),IGP(NGRP+1),JGP(N) C LOGICAL ORDER C .. EXTERNAL SUBROUTINES .. EXTERNAL F,PDERV,ODRV,CDRV C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSQRT C .. C .. COMMON BLOCKS .. COMMON /ORDERING/ORDER INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP c c NPSET = NPSET + 1 IF (NRENEW.EQ.0) THEN IF (MITER.EQ.5 .OR. MITER .EQ. 6) THEN KK=1 KMIN = IA(1) DO 275 J=1,N KMAX=IA(1+J)-1 DO 270 K=KMIN,KMAX I =JA(K) PIJ=-PWCOPY(K) IF ( I .EQ. J) PIJ = PIJ +1.0D0/CON PW(K) = PIJ 270 CONTINUE KMIN=KMAX+1 275 CONTINUE ENDIF GO TO 131 ENDIF NJE = NJE + 1 70 CONTINUE C C.... FORM THE ITERATION MATRIX P =I/CON - J C.... MITER = 5, CALL PDERV ROUTINE AND FORM P C... IF (MITER .EQ. 5) THEN KMIN = IA(1) DO 130 J = 1,N KMAX = IA(1+J) -1 DO 110 I = 1,N 110 PWCOPY(I) = 0.0D0 CALL PDERV(N,T,Y,J,PWCOPY,IPAR,RPAR,IERR) DO 120 K = KMIN,KMAX II = JA(K) PW(K) = -PWCOPY(II) IF (II .EQ. J) PW(K) = PW(K) + 1.0D0/CON 120 CONTINUE KMIN = KMAX + 1 130 CONTINUE ELSE C.... C.... NUMERICAL APRROXIMATION TO THE SPARSE JACOBIAN, MAKE NGRP CALS C.... TO F TO APPROXIMATE THE JACOBIAN AND FORM THE MATRIX P C.... DO 35 I=1,N IF(ITOL.EQ.2) THEN YMAX(I)=dabs(Y(I,1))*RTOL(1)+ATOL(1) ELSE IF (ITOL.EQ.3) THEN YMAX(I)=dabs(Y(I,1))*RTOL(1)+ATOL(I) ELSE IF(ITOL.EQ.4) THEN YMAX(I)=dabs(Y(I,1))*RTOL(I)+ATOL(1) ELSE IF(ITOL.EQ.5) THEN YMAX(I)=dabs(Y(I,1))*RTOL(I)+ATOL(I) END IF 35 CONTINUE D = 0.D0 DO 41 I = 1,N D = D + (SAVE2(I)*YMAX(I))**2 41 CONTINUE IF (ITOL .EQ. 1) D= D*RTOL(1)**2 D = DSQRT(D/DFLOAT(N)) R0 = DABS(H)*D*(1.0D+03)*UROUND*DFLOAT(N) IF (R0.EQ.0.D0) R0 = 1.D0 JMIN = IGP(1) DO 240 NG =1,NGRP JMAX= IGP(1+NG)-1 DO 210 J=JMIN, JMAX JJ= JGP(J) SAVE1(JJ) =Y(JJ,1) YI = Y(JJ,1) R = DMAX1(EPSJAC*DABS(YI),R0/YMAX(JJ)) 210 Y(JJ,1)=Y(JJ,1)+R CALL F(N,T,Y,WRKSPC,IPAR,RPAR,IERR) DO 230 J=JMIN,JMAX JJ = JGP(J) Y(JJ,1)=SAVE1(JJ) YJJ = Y(JJ,1) R = DMAX1(EPSJAC*DABS(YJJ),R0/YMAX(JJ)) D = CON/R K1 = IA(JJ) K2 = IA(JJ+1) -1 DO 220 K=K1,K2 I = JA(K) TEMPRY = WRKSPC(I) - SAVE2(I) PWCOPY(k)=TEMPRY/R PW(k) = -PWCOPY(k) IF (I .EQ. JJ ) PW(K) = PW(K) + 1.0D0/CON 220 CONTINUE 230 CONTINUE JMIN = JMAX + 1 240 CONTINUE NFE = NFE + NGRP ENDIF 131 CONTINUE C.... C.... NUMERICAL FACTORIZATION OF THE ITERATION MATRIX P C.... DO 132 I=1,N 132 WRKSPC(I) = 0.0D0 CALL CDRV(N,P,P,IP,IA,JA,PW,WRKSPC,WRKSPC,NSP,ISP, + RSP,IEXTRA,2,IFLAG) C IF (IFLAG .NE. 0) THEN WRITE(LOUT,1333) IFALG 1333 FORMAT('FATAL ERROR IN THE CDRV ROUTINE (SPARSE + NUMERICAL FACTORIZATION) IFLAG RETURN IS ',I3) IER = 1 ORDER = .FALSE. ENDIF NDEC = NDEC + 1 C RETURN END C ---------------------- END OF SUBROUTINE PSET --------------------- 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) 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),IPAR(*),RPAR(*) 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----------------------------------------------------------------------------- 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,RSP,IA,JA,P,IP,NSP, + ISP,NONZ,LOUT,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(nonz),SAVE1(N),SAVE2(N),Y(N,12), & YMAX(N),RTOL(1),ATOL(1),SCALE(N),RSP(NSP),RPAR(*) INTEGER P(N),IP(N),ISP(NSP),IA(N+1),JA(NONZ),IPAR(*) INTEGER MAXORD,NBSOL,NCOSET,NDEC,NFE,NJE,NPSET,NQUSED,NSTEP C C .. EXTERNAL SUBROUTINES .. EXTERNAL F,CDRV C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DMAX1,DMIN1 C .. C .. DATA STATEMENTS .. DATA ZERO/0.0D+0/ C .. 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(LMB.EQ.1) GOTO 25 DO 10 I = 1,N SAVE1(I) = (-SAVE1(I)+HBETA*SAVE2(I)+ARH(I))/QQQ 10 CONTINUE IFLAG = 0 IF ( MF .EQ. 25 .OR. MF .EQ. 26 ) THEN CALL CDRV(N,P,P,IP,IA,JA,PW,SAVE1,SAVE1,NSP,ISP, + RSP,IEXTRA,4,IFLAG) IF (IFLAG .NE. 0) THEN WRITE(LOUT,2222) IFLAG 2222 FORMAT('FALTAL ERROR IN THE CDRV ROUTINE + (SPARSE SOLVER) IFLAG RETURN IS ' , I3) STOP ENDIF 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 DO 40 I = 1,N SAVE1(I)=(-SAVE1(I)+HBETA*SAVE2(I)+ARH(I))/qqq 40 CONTINUE C C IF WE ARE HERE THEN PARTIALS ARE O.K. C IFLAG = 0 IF( MF.EQ. 25 .OR. MF .EQ. 26) THEN CALL CDRV(N,P,P,IP,IA,JA,PW,SAVE1,SAVE1,NSP,ISP, + RSP,IEXTRA,4,IFLAG) IF (IFLAG .NE. 0) THEN WRITE(LOUT,3333) IFLAG 3333 FORMAT('FALTAL ERROR IN THE CDRV ROUTINE + (SPARSE SOLVER) IFLAG RETURN IS ', I3) STOP ENDIF 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----------------------------------------------------------------------------- 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,NQUSED,NSTEP,NFAIL,NFE, + NJE,NDEC,NBSOL,NPSET,NCOSET,MAXORD,UROUND,EPSJAC,HUSED,RSP, + IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGP,JGP,IGP,IPAR,RPAR,IERR) 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),RSP(NSP),IPAR(*),RPAR(*) INTEGER P(N),IP(N),IA(N+1),JA(NONZ),ISP(NSP),JGP(N),IGP(NONZ) C .. C .. LOCAL SCALARS .. INTEGER I,IER,IITER,IITER2,IREDO,ITST,J,J1,J2,KFAIL,LL,LMAX, + M3STEP,MAXDER,METH,MFOLD,MITER,NEWQ, NSP C 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,CDRV 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,NT 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 .. 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 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 ENDIF 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 ENDIF 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 ENDIF ENDIF TCRAT1 = ZERO TCRAT2 = ZERO IF (CFAIL) THEN NRENEW = 1 NEWPAR = 1 JSINUP = -1 JNEWIM = .TRUE. ENDIF CFAIL = .FALSE. JSNOLD = 0 MQ1TMP = MEQC1 MQ2TMP = MEQC2 CALL PSET(Y,N,H,T,UROUND,EPSJAC,QI,MITER,IER,F,PDERV,NRENEW,YMAX, + SAVE1,SAVE2,PW,PWCOPY,ERROR,ITOL,RTOL,ATOL,NPSET,NJE,NFE, + NDEC,RSP,IA,JA,P,IP,NSP,ISP,NONZ,SAVF,NGP,JGP,IGP,LOUT, + IPAR,RPAR,IERR) C IF (IERR .NE. 0 ) GOTO 8000 QQQ=QI 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, + RSP,IA,JA,P,IP,NSP,ISP,NONZ,LOUT,IPAR,RPAR,IERR) 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, + RSP,IA,JA,P,IP,NSP,ISP,NONZ,LOUT,IPAR,RPAR,IERR) 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 DO 145 I=1,N SAVE2(I) = Y(I,1)*QQ Y(I,1) = SAVE1(I) 145 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 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, + RSP,IA,JA,P,IP,NSP,ISP,NONZ,LOUT,IPAR,RPAR,IERR) 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 DO 233 I=1,N Y(I,LMP2) = Y(I,LMP3)*QQ Y(I,LMP3) = SAVE1(I) 233 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 DO 310 I=1,N SAVE1(I) = (-Y(I,1) + QI*SAVE2(I) + ARH(I))/qqq 310 CONTINUE C IF (MF .EQ. 25 .OR. MF .EQ. 26) THEN IFLAG = 0 CALL CDRV(N,P,P,IP,IA,JA,PW,SAVE1,SAVE1,NSP,ISP, + RSP,IEXTRA,4,IFLAG) IF (IFLAG .NE. 0) THEN WRITE(LOUT,4444) IFLAG 4444 FORMAT('FALTAL ERROR IN CDRV ROUTINE + (SPARSE SOLVER) IFLAG RETURN IS ', I3) STOP ENDIF NBSOL = NBSOL + 1 ENDIF 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 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.4) THEN c WRITE (LOUT,9000) 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 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(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 DDOWN=ZERO 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 (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 ) RETURN 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----------------------------------------------------------------------------- C THE NEXT ROUTINE PREFORMS PREPROCESSING RELATED TO THE SPARSE LINEAR C SYSTEMS. C THE OPERATIONS THAT ARE PERFORMED HERE ARE C * COMPUTE SPARSENESS STRUCTURE OF THE JACOBIAN MATRIX ACCORDING TO MF, C * COMPUTE GROUPING OF COLUMN INDICES (MF = 26), C * COMPUTE A NEW ORDERING OF ROWS AND COLUMNS OF THE MATRIX, C * REORDER JA CORRESPONDING TO THE NEW ORDERING, C * PERFORM A SYMBOLIC LU FACTORIZATION OF THE MATRIX. C----------------------------------------------------------------------------- SUBROUTINE SPARSDET(N,T,Y0,PWCOPY,SCALE,IA,JA,SAVF,P,IP,NSP, + ISP,RSP,NONZ,NGP,JGP,IGP,MF,INCL,JDONE,LOUT,F,NNZ, + IPAR,RPAR,IERR) C ..ARRAY ARGUMENTS .. DOUBLE PRECISION PWCOPY(NONZ),Y0(N),RSP(NSP),SAVF(N),YJ DOUBLE PRECISION FAC,SCALE(N),FTEM(1),DQ,RPAR(*) INTEGER P(N),IP(N),IA(N+1),JA(NONZ),ISP(NSP), JGP(N),IGP(NONZ), + INCL(N),JDONE(N),IPAR(*) C LOGICAL ORDER C .. EXTERNAL SUBROUTINES .. EXTERNAL F,PDERV,ODRV,CDRV C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1 C .. C .. COMMON BLOCKS .. COMMON /ORDERING/ORDER C... C... PERTURB THE INITIAL CONDITIONS, Y0, FOR STRUCTURE DETERMINATION. C... DO I=1,N FAC= 1.0D0+1.0D0/(DFLOAT(I)+1.D0) Y0(I) =Y0(I) +FAC*DSIGN(SCALE(I),Y0(I)) ENDDO C... C... COMPUTE THE SPARSITY STRUCTURE FROM THE USER SUPPLIES JACOBIAN C... ROUTINE PDERV. C... IF (MF .EQ. 25) THEN KK = 1 IA(1) =1 DO 111 J =1,N JA(KK) = J KK = KK+ 1 DO K=1,N SAVF(K) = 0.0D0 ENDDO CALL PDERV(N,T,Y0,J,SAVF,IPAR,RPAR,IERR) DO 112 II=1,N IF (DABS(SAVF(II)) .LE. 0.0D0 ) GOTO 112 IF (II .EQ. J) GOTO 112 JA(KK) = II KK=KK+1 112 CONTINUE IA(1+J) = KK 111 CONTINUE ELSE C.... C.... MF = 26 COMPUTE SPARSITY STRUCTURE FROM RESULTS OF N+1 CALLS TO F C.... 100 KK = 1 IA(1) = 1 CALL F(N,T,Y0,SAVF,IPAR,RPAR,IERR) DO 120 J = 1,N JA(KK) = J KK = KK + 1 YJ = Y0(J) DYJ = DSIGN(SCALE(J),YJ) Y0(J) = YJ + DYJ CALL F(N,T,Y0,PWCOPY,IPAR,RPAR,IERR) Y0(J) = YJ DO 110 II = 1,N DQ = (PWCOPY(II) - SAVF(II))/DYJ IF (DABS(DQ) .LE. 0.D0) GOTO 110 IF (II .EQ. J) GOTO 110 JA(KK) = II KK = KK + 1 110 CONTINUE IA(1+J) = KK 120 CONTINUE ENDIF C..... NNZ = IA(1+N)-1 IF (NNZ .GT. NONZ) THEN WRITE(LOUT,*) 'NUMBER OF NONZEROS ',NNZ, 'NOT',NONZ STOP ENDIF C.... C.... COMPUTE COLUMN GROUPING FOR SPARSE NUMERICAL JACOBIAN I.E. MF = 26 C.... IF ( MF .EQ. 26) THEN MAXG = N+1 CALL JGROUP(N,IA,JA,MAXG,NGP,IGP,JGP,INCL,JDONE,IER,NNZ) ENDIF C.... C.... COMPUTE NEW ORDERING OF NONZERO ENTERIS IN THE ITERATION MATRIX C.... DO I=1,N P(I) = I ENDDO CALL ODRV(N,IA,JA,PWCOPY,P,IP,NSP,ISP,1,IFLAG) IF (IFLAG .NE. 0) THEN WRITE(LOUT,*)'FATAL ERROR IN THE ODRV ROUTINE (SPARSE ORDERING) + IFALG RETURN IS= ',IFLAG STOP ENDIF C.... C.... REORDER AND DO SYMBOLIC LU FACTORIZATION OF THE ITERATION MATRIX C.... DO I=1, NNZ PWCOPY(I) =0.0D0 ENDDO CALL CDRV(N,P,P,IP,IA,JA,PWCOPY,FTEM,FTEM,NSP,ISP, + RSP,IEXTRA,5,IFLAG) IF (IFLAG .EQ. 0) THEN ORDER = .TRUE. ELSE WRITE(LOUT,1111) IFLAG 1111 FORMAT('FATAL ERROR IN THE CDRV ROUTINE (SYMBOLIC + FACTORIZATION) IFLAG RETURN IS ', I3) STOP ENDIF C RETURN END C----------------------------------------------------------------------------- SUBROUTINE JGROUP(N,IA,JA,MAXG,NGRP,IGP,JGP,INCL,JDONE,IER,NNZ) INTEGER N, IA, JA, MAXG, NGRP, IGP, JGP, INCL, JDONE, IER DIMENSION IA(1), JA(1), IGP(1), JGP(N), INCL(N), JDONE(N) C----------------------------------------------------------------------------- C THIS SUNROUTINE CONSTRUCTS GROUPING OF THE COLUMN INDICES OF THE C JACOBIAN MATRIX, USED IN THE NUMERICAL EVALUATION OF THE JACOBIAN BY C FINITE DIFFERENCES. C C INPUT.. C N = THE ORDER OF THE MATRIX. C IA,JA = SPARSE STRUCURE DESCRIPTORS OF THE MATRIX BY ROWS. C MAXG = LENGTH OF AVAILABLE STORAGE IN THE IGP ARRAY. C C OUTPUT.. C NGRP = NUMBER OF GROUPS. C JGP = ARRAY OF LENGTH N CONTAINING THE COLUMN INDICES BY GROUPS. C IGP = POINTER ARRAY OF LENGTH NGRP + 1 TO THE LOCATIONS IN JGP OF C THE BEGINNING OF EACH GROUP. C IER = ERROR INDICATOR. IER = 0 IF NO ERROR OCCURRED, OR 1 IF C MAXG WAS INSUFFICIENT. C C INCL AND JDONE ARE WORKING ARRAYS OF LENGTH N. C----------------------------------------------------------------------------- INTEGER I, J, K, KMIN, KMAX, NCOL, NG C IER = 0 DO 10 JK = 1,N 10 JDONE(JK) = 0 NCOL = 1 DO 60 NG = 1,MAXG IGP(NG) = NCOL DO 20 I = 1,N 20 INCL(I) = 0 DO 50 J = 1,N C... C... REJECT COLUMN J IF IT IS ALREADY IN A GROUP. C... IF (JDONE(J) .EQ. 1) GOTO 50 KMIN = IA(J) KMAX = IA(J+1) - 1 DO 30 K = KMIN, KMAX C... C... REJECT COLUMN J IF IT OVERLAPS ANY COLUMN ALREADY IN THIS GROUP. C... I = JA(K) IF (INCL(I) .EQ. 1) GOTO 50 30 CONTINUE C... C... ACCEPT COLUMN J INTO GROUP NG. C... JGP(NCOL) = J NCOL = NCOL + 1 JDONE(J) = 1 DO 40 K = KMIN,KMAX I = JA(K) INCL(I) = 1 40 CONTINUE 50 CONTINUE C... C... STOP IF THIS GROUP IS EMPTY (GROUPING IS COMPLETE). C... IF (NCOL .EQ. IGP(NG)) GOTO 70 60 CONTINUE C... C... ERROR RETURN IF NOT ALL COLUMNS WERE CHOSEN (MAXG TOO SMALL). C... IF (NCOL.LE. N) GOTO 80 NG = MAXG 70 NGRP = NG - 1 RETURN 80 IER = 1 RETURN END C----------------------- END OF SUBROUTINE JGROUP ---------------------- 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----------------------------------------------------------------------------