Mercurial > hg > octave-avbm
changeset 3987:258c1d15ad78
[project @ 2002-07-11 19:33:35 by tenny]
changed real to double precision
author | tenny |
---|---|
date | Thu, 11 Jul 2002 19:33:35 +0000 |
parents | 7cb85d5c7aad |
children | 9a73d0e7e737 |
files | libcruft/odessa/odessa.f libcruft/odessa/odessa_cfode.f libcruft/odessa/odessa_ewset.f libcruft/odessa/odessa_intdy.f libcruft/odessa/odessa_prepd.f libcruft/odessa/odessa_prepj.f libcruft/odessa/odessa_solsy.f libcruft/odessa/odessa_stesa.f libcruft/odessa/odessa_stode.f libcruft/odessa/odessa_vnorm.f |
diffstat | 10 files changed, 86 insertions(+), 77 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/odessa/odessa.f +++ b/libcruft/odessa/odessa.f @@ -1468,7 +1468,7 @@ DO 80 I = 1,N 80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1) C RELOAD WM(1) = RWORK(LWM), SINCE LWM MAY HAVE CHANGED. --------------- - 90 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) + 90 IF (MITER .GT. 0) RWORK(LWM) = DSQRT(UROUND) GO TO 200 C----------------------------------------------------------------------- C BLOCK C. @@ -1486,7 +1486,7 @@ IF (H0 .NE. ZERO .AND. (T + H0 - TCRIT)*H0 .GT. ZERO) 1 H0 = TCRIT - T 105 JSTART = 0 - IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) + IF (MITER .GT. 0) RWORK(LWM) = DSQRT(UROUND) NHNIL = 0 NST = 0 NJE = 0 @@ -1547,29 +1547,29 @@ C THE SIGN OF H0 IS INFERRED FROM THE INITIAL VALUES OF TOUT AND T. C----------------------------------------------------------------------- 125 IF (H0 .NE. ZERO) GO TO 180 - TDIST = ABS(TOUT - T) - W0 = MAX(ABS(T),ABS(TOUT)) + TDIST = DABS(TOUT - T) + W0 = DMAX1(DABS(T),DABS(TOUT)) IF (TDIST .LT. TWO*UROUND*W0) GO TO 622 TOL = RTOL(1) IF (ITOL .LE. 2) GO TO 140 DO 130 I = 1,N - 130 TOL = MAX(TOL,RTOL(I)) + 130 TOL = DMAX1(TOL,RTOL(I)) 140 IF (TOL .GT. ZERO) GO TO 160 ATOLI = ATOL(1) DO 150 I = 1,N IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) - AYI = ABS(Y(I)) - IF (AYI .NE. ZERO) TOL = MAX(TOL,ATOLI/AYI) + AYI = DABS(Y(I)) + IF (AYI .NE. ZERO) TOL = DMAX1(TOL,ATOLI/AYI) 150 CONTINUE - 160 TOL = MAX(TOL,100.0D0*UROUND) - TOL = MIN(TOL,0.001D0) + 160 TOL = DMAX1(TOL,100.0D0*UROUND) + TOL = DMIN1(TOL,0.001D0) SUM = ODESSA_VNORM (N, RWORK(LF0), RWORK(LEWT)) SUM = ONE/(TOL*W0*W0) + TOL*SUM**2 - H0 = ONE/SQRT(SUM) + H0 = ONE/DSQRT(SUM) H0 = MIN(H0,TDIST) - H0 = SIGN(H0,TOUT-T) + H0 = DSIGN(H0,TOUT-T) C ADJUST H0 IF NECESSARY TO MEET HMAX BOUND. --------------------------- - 180 RH = ABS(H0)*HMXI + 180 RH = DABS(H0)*HMXI IF (RH .GT. ONE) H0 = H0/RH C LOAD H WITH H0 AND SCALE YH(*,2) BY H0. ------------------------------ H = H0 @@ -1602,8 +1602,8 @@ GO TO 420 240 TCRIT = RWORK(1) IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624 - 245 HMX = ABS(TN) + ABS(H) - IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX + 245 HMX = DABS(TN) + DABS(H) + IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX IF (IHIT) GO TO 400 TNEXT = TN + H*(ONE + FOUR*UROUND) IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250 @@ -1677,8 +1677,8 @@ CALL ODESSA_INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) T = TOUT GO TO 420 - 345 HMX = ABS(TN) + ABS(H) - IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX + 345 HMX = DABS(TN) + DABS(H) + IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX IF (IHIT) GO TO 400 TNEXT = TN + H*(ONE + FOUR*UROUND) IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250 @@ -1686,8 +1686,8 @@ JSTART = -2 GO TO 250 C ITASK = 5. SEE IF TCRIT WAS REACHED AND JUMP TO EXIT. --------------- - 350 HMX = ABS(TN) + ABS(H) - IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX + 350 HMX = DABS(TN) + DABS(H) + IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX C----------------------------------------------------------------------- C BLOCK G. C THE FOLLOWING BLOCK HANDLES ALL SUCCESSFUL RETURNS FROM ODESSA. @@ -1770,7 +1770,7 @@ 560 BIG = ZERO IMXER = 1 DO 570 I = 1,NYH - SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1)) + SIZE = DABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1)) IF (BIG .GE. SIZE) GO TO 570 BIG = SIZE IMXER = I
--- a/libcruft/odessa/odessa_cfode.f +++ b/libcruft/odessa/odessa_cfode.f @@ -47,9 +47,9 @@ C INITIALLY, P(X) = 1. C----------------------------------------------------------------------- RQ1FAC = RQFAC - RQFAC = RQFAC/REAL(NQ) + RQFAC = RQFAC/DBLE(NQ) NQM1 = NQ - 1 - FNQM1 = REAL(NQM1) + FNQM1 = DBLE(NQM1) NQP1 = NQ + 1 C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ---------------------------------- PC(NQ) = ZERO @@ -63,17 +63,17 @@ TSIGN = ONE DO 120 I = 2,NQ TSIGN = -TSIGN - PINT = PINT + TSIGN*PC(I)/REAL(I) - 120 XPIN = XPIN + TSIGN*PC(I)/REAL(I+1) + PINT = PINT + TSIGN*PC(I)/DBLE(I) + 120 XPIN = XPIN + TSIGN*PC(I)/DBLE(I+1) C STORE COEFFICIENTS IN ELCO AND TESCO. -------------------------------- ELCO(1,NQ) = PINT*RQ1FAC ELCO(2,NQ) = ONE DO 130 I = 2,NQ - 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/REAL(I) + 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/DBLE(I) AGAMQ = RQFAC*XPIN RAGQ = ONE/AGAMQ TESCO(2,NQ) = RAGQ - IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/REAL(NQP1) + IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/DBLE(NQP1) TESCO(3,NQM1) = RAGQ 140 CONTINUE RETURN @@ -86,7 +86,7 @@ C P(X) = (X+1)*(X+2)*...*(X+NQ). C INITIALLY, P(X) = 1. C----------------------------------------------------------------------- - FNQ = REAL(NQ) + FNQ = DBLE(NQ) NQP1 = NQ + 1 C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------ PC(NQP1) = ZERO @@ -99,8 +99,8 @@ 220 ELCO(I,NQ) = PC(I)/PC(2) ELCO(2,NQ) = ONE TESCO(1,NQ) = RQ1FAC - TESCO(2,NQ) = REAL(NQP1)/ELCO(1,NQ) - TESCO(3,NQ) = REAL(NQ+2)/ELCO(1,NQ) + TESCO(2,NQ) = DBLE(NQP1)/ELCO(1,NQ) + TESCO(3,NQ) = DBLE(NQ+2)/ELCO(1,NQ) RQ1FAC = RQ1FAC/FNQ 230 CONTINUE RETURN
--- a/libcruft/odessa/odessa_ewset.f +++ b/libcruft/odessa/odessa_ewset.f @@ -12,7 +12,7 @@ DO 10 I = 1,N IF (ITOL .GE. 3) RTOLI = RTOL(I) IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I) - EWT(I) = RTOLI*ABS(YCUR(I)) + ATOLI + EWT(I) = RTOLI*DABS(YCUR(I)) + ATOLI 10 CONTINUE RETURN C----------------------- END OF SUBROUTINE ODESSA_EWSET -----------------------
--- a/libcruft/odessa/odessa_intdy.f +++ b/libcruft/odessa/odessa_intdy.f @@ -37,7 +37,7 @@ JJ1 = L - K DO 10 JJ = JJ1,NQ 10 IC = IC*JJ - 15 C = REAL(IC) + 15 C = DBLE(IC) DO 20 I = 1,NYH 20 DKY(I) = C*YH(I,L) IF (K .EQ. NQ) GO TO 55 @@ -50,7 +50,7 @@ JJ1 = JP1 - K DO 30 JJ = JJ1,J 30 IC = IC*JJ - 35 C = REAL(IC) + 35 C = DBLE(IC) DO 40 I = 1,NYH 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) 50 CONTINUE
--- a/libcruft/odessa/odessa_prepd.f +++ b/libcruft/odessa/odessa_prepd.f @@ -36,7 +36,7 @@ GO TO (100, 200), IDF1 C IDF = 0, CALL F TO APPROXIMATE DFDP. --------------------------------- 100 RPAR = PAR(JPAR) - R = MAX(SRUR*ABS(RPAR),SRUR) + R = DMAX1(SRUR*DABS(RPAR),SRUR) PAR(JPAR) = RPAR + R FAC = 1.0D0/R CALL F (NEQ, TN, Y, PAR, FTEM) @@ -52,3 +52,12 @@ RETURN C -------------------- END OF SUBROUTINE ODESSA_PREPDF ------------------------ END + + + + + + + + +
--- a/libcruft/odessa/odessa_prepj.f +++ b/libcruft/odessa/odessa_prepj.f @@ -63,13 +63,13 @@ GO TO 240 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. -------------------- 200 FAC = ODESSA_VNORM (N, SAVF, EWT) - R0 = 1000.0D0*ABS(H)*UROUND*REAL(N)*FAC + R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC IF (R0 .EQ. ZERO) R0 = ONE SRUR = WM(1) J1 = 2 DO 230 J = 1,N YJ = Y(J) - R = MAX(SRUR*ABS(YJ),R0/EWT(J)) + R = DMAX1(SRUR*DABS(YJ),R0/EWT(J)) Y(J) = Y(J) + R FAC = -HL0/R CALL F (NEQ, TN, Y, PAR, FTEM) @@ -100,8 +100,8 @@ R0 = H*SAVF(I) - YH(I,2) DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I)) WM(I+2) = 1.0D0 - IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320 - IF (ABS(DI) .EQ. ZERO) GO TO 330 + IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320 + IF (DABS(DI) .EQ. ZERO) GO TO 330 WM(I+2) = 0.1D0*R0/DI 320 CONTINUE RETURN @@ -126,26 +126,26 @@ 500 ML = IWM(1) MU = IWM(2) MBAND = ML + MU + 1 - MBA = MIN(MBAND,N) + MBA = MIN0(MBAND,N) MEBAND = MBAND + ML MEB1 = MEBAND - 1 SRUR = WM(1) FAC = ODESSA_VNORM (N, SAVF, EWT) - R0 = 1000.0D0*ABS(H)*UROUND*REAL(N)*FAC + R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC IF (R0 .EQ. ZERO) R0 = ONE DO 560 J = 1,MBA DO 530 I = J,N,MBAND YI = Y(I) - R = MAX(SRUR*ABS(YI),R0/EWT(I)) + R = DMAX1(SRUR*DABS(YI),R0/EWT(I)) 530 Y(I) = Y(I) + R CALL F (NEQ, TN, Y, PAR, FTEM) DO 550 JJ = J,N,MBAND Y(JJ) = YH(JJ,1) YJJ = Y(JJ) - R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ)) + R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ)) FAC = -HL0/R - I1 = MAX(JJ-MU,1) - I2 = MIN(JJ+ML,N) + I1 = MAX0(JJ-MU,1) + I2 = MIN0(JJ+ML,N) II = JJ*MEB1 - ML + 2 DO 540 I = I1,I2 540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC
--- a/libcruft/odessa/odessa_solsy.f +++ b/libcruft/odessa/odessa_solsy.f @@ -43,7 +43,7 @@ R = HL0/PHL0 DO 320 I = 1,N DI = ONE - R*(ONE - ONE/WM(I+2)) - IF (ABS(DI) .EQ. ZERO) GO TO 390 + IF (DABS(DI) .EQ. ZERO) GO TO 390 320 WM(I+2) = ONE/DI 330 DO 340 I = 1,N 340 X(I) = WM(I+2)*X(I)
--- a/libcruft/odessa/odessa_stesa.f +++ b/libcruft/odessa/odessa_stesa.f @@ -110,8 +110,8 @@ IF (L .EQ. LMAX) GO TO 70 DO 60 I= 1,N 60 Y(I,J) = ACOR(I,J) - YH(I,J,LMAX) - DUPS = MAX(DUPS,ODESSA_VNORM(N,Y(1,J),EWT(1,J))*TI3) - 70 DSMS = MAX(DSMS,ERR) + DUPS = DMAX1(DUPS,ODESSA_VNORM(N,Y(1,J),EWT(1,J))*TI3) + 70 DSMS = DMAX1(DSMS,ERR) 100 CONTINUE RETURN C-----------------------------------------------------------------------
--- a/libcruft/odessa/odessa_stode.f +++ b/libcruft/odessa/odessa_stode.f @@ -164,14 +164,14 @@ NQNYH = NQ*NYH RC = RC*EL(1)/EL0 EL0 = EL(1) - CONIT = 0.5D0/REAL(NQ+2) + CONIT = 0.5D0/DBLE(NQ+2) DDN = ODESSA_VNORM (N, SAVF, EWT)/TESCO(1,L) - EXDN = ONE/REAL(L) + EXDN = ONE/DBLE(L) RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0) - RH = MIN(RHDN,ONE) + RH = DMIN1(RHDN,ONE) IREDO = 3 IF (H .EQ. HOLD) GO TO 170 - RH = MIN(RH,ABS(H/HOLD)) + RH = DMIN1(RH,DABS(H/HOLD)) H = HOLD GO TO 175 C----------------------------------------------------------------------- @@ -185,7 +185,7 @@ NQNYH = NQ*NYH RC = RC*EL(1)/EL0 EL0 = EL(1) - CONIT = 0.5D0/REAL(NQ+2) + CONIT = 0.5D0/DBLE(NQ+2) GO TO (160, 170, 200), IRET C----------------------------------------------------------------------- C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST @@ -198,9 +198,9 @@ H = HOLD IREDO = 3 GO TO 175 - 170 RH = MAX(RH,HMIN/ABS(H)) - 175 RH = MIN(RH,RMAX) - RH = RH/MAX(ONE,ABS(H)*HMXI*RH) + 170 RH = DMAX1(RH,HMIN/DABS(H)) + 175 RH = DMIN1(RH,RMAX) + RH = RH/DMAX1(ONE,DABS(H)*HMXI*RH) R = ONE DO 180 J = 2,L R = R*RH @@ -219,7 +219,7 @@ C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS FOR ISOPT = 0, C AND AT LEAST ONCE EVERY STEP FOR ISOPT = 1. C----------------------------------------------------------------------- - 200 IF (ABS(RC-ONE) .GT. CCMAX) IPUP = MITER + 200 IF (DABS(RC-ONE) .GT. CCMAX) IPUP = MITER IF (NST .GE. NSLP+MSBP) IPUP = MITER TN = TN + H I1 = NQNYH + 1 @@ -288,8 +288,8 @@ 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----------------------------------------------------------------------- - 400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP) - DCON = DEL*MIN(ONE,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) + 400 IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP) + DCON = DEL*DMIN1(ONE,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) IF (DCON .LE. ONE) GO TO 450 M = M + 1 IF (M .EQ. MAXCOR) GO TO 410 @@ -320,7 +320,7 @@ 440 YH1(I) = YH1(I) - YH1(I+NYH) 445 CONTINUE IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 - IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670 + IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670 IF (NCF .EQ. MXNCF) GO TO 670 RH = 0.25D0 IPUP = MITER @@ -393,7 +393,7 @@ 510 YH1(I) = YH1(I) - YH1(I+NYH) 515 CONTINUE RMAX = 2.0D0 - IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660 + IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660 IF (KFLAG .LE. -3) GO TO 640 IREDO = 2 RHUP = ZERO @@ -415,23 +415,23 @@ DO 530 I = 1,N 530 SAVF(I) = ACOR(I) - YH(I,LMAX) DUP = ODESSA_VNORM (N, SAVF, EWT)/TESCO(3,NQ) - DUP = MAX(DUP,DUPS) - EXUP = ONE/REAL(L+1) + DUP = DMAX1(DUP,DUPS) + EXUP = ONE/DBLE(L+1) RHUP = ONE/(1.4D0*DUP**EXUP + 0.0000014D0) - 540 EXSM = ONE/REAL(L) - DSM = MAX(DSM,DSMS) + 540 EXSM = ONE/DBLE(L) + DSM = DMAX1(DSM,DSMS) RHSM = ONE/(1.2D0*DSM**EXSM + 0.0000012D0) RHDN = ZERO IF (NQ .EQ. 1) GO TO 560 JPOINT = 1 DO 550 J = 1,NSV DDN = ODESSA_VNORM (N, YH(JPOINT,L), EWT(JPOINT))/TESCO(1,NQ) - DDNS = MAX(DDNS,DDN) + DDNS = DMAX1(DDNS,DDN) JPOINT = JPOINT + N 550 CONTINUE DDN = DDNS DDNS = ZERO - EXDN = ONE/REAL(NQ) + EXDN = ONE/DBLE(NQ) RHDN = ONE/(1.3D0*DDN**EXDN + 0.0000013D0) 560 IF (RHSM .GE. RHUP) GO TO 570 IF (RHUP .GT. RHDN) GO TO 590 @@ -447,14 +447,14 @@ 590 NEWQ = L RH = RHUP IF (RH .LT. 1.1D0) GO TO 610 - R = EL(L)/REAL(L) + R = EL(L)/DBLE(L) DO 600 I = 1,NYH 600 YH(I,NEWQ+1) = ACOR(I)*R GO TO 630 610 IALTH = 3 GO TO 700 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610 - IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) + IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0) C----------------------------------------------------------------------- C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED. @@ -476,7 +476,7 @@ C----------------------------------------------------------------------- 640 IF (KFLAG .EQ. -10) GO TO 660 RH = 0.1D0 - RH = MAX(HMIN/ABS(H),RH) + RH = DMAX1(HMIN/DABS(H),RH) H = H*RH DO 645 I = 1,NYH 645 Y(I) = YH(I,1)
--- a/libcruft/odessa/odessa_vnorm.f +++ b/libcruft/odessa/odessa_vnorm.f @@ -34,38 +34,38 @@ I = 1 20 SX = V(I)*W(I) GO TO (30,40,70,80),NEXT -30 IF (ABS(SX).GT.CUTLO) GO TO 110 +30 IF (DABS(SX).GT.CUTLO) GO TO 110 NEXT = 2 XMAX = ZERO 40 IF (SX.EQ.ZERO) GO TO 130 - IF (ABS(SX).GT.CUTLO) GO TO 110 + IF (DABS(SX).GT.CUTLO) GO TO 110 NEXT = 3 GO TO 60 50 I=J NEXT = 4 SUM = (SUM/SX)/SX -60 XMAX = ABS(SX) +60 XMAX = DABS(SX) GO TO 90 -70 IF(ABS(SX).GT.CUTLO) GO TO 100 -80 IF(ABS(SX).LE.XMAX) GO TO 90 +70 IF(DABS(SX).GT.CUTLO) GO TO 100 +80 IF(DABS(SX).LE.XMAX) GO TO 90 SUM = ONE + SUM * (XMAX/SX)**2 - XMAX = ABS(SX) + XMAX = DABS(SX) GO TO 130 90 SUM = SUM + (SX/XMAX)**2 GO TO 130 100 SUM = (SUM*XMAX)*XMAX -110 HITEST = CUTHI/REAL(N) +110 HITEST = CUTHI/DBLE(N) DO 120 J = I,N SX = V(J)*W(J) - IF(ABS(SX).GE.HITEST) GO TO 50 + IF(DABS(SX).GE.HITEST) GO TO 50 SUM = SUM + SX**2 120 CONTINUE - ODESSA_VNORM = SQRT(SUM) + ODESSA_VNORM = DSQRT(SUM) GO TO 140 130 CONTINUE I = I + 1 IF (I.LE.N) GO TO 20 - ODESSA_VNORM = XMAX * SQRT(SUM) + ODESSA_VNORM = XMAX * DSQRT(SUM) 140 CONTINUE RETURN C----------------------- END OF FUNCTION ODESSA_VNORM -------------------------