--- cash.orig/mebdfdae.f 2007-12-15 15:37:46.000000000 -0500 +++ cash/mebdfdae.f 2014-03-02 16:20:04.683064613 -0500 @@ -53,13 +53,13 @@ C C NOVEMBER 6th 1998: FIRST RELEASE C -C OVDRIV +C A_OVDRIV C A PACKAGE FOR THE SOLUTION OF THE INITIAL VALUE PROBLEM C FOR SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS C DY/DT = F(Y,T), Y=(Y(1),Y(2),Y(3), . . . ,Y(N)) C AND LINEARLY IMPLICIT DIFFERENTIAL ALGEBRAIC EQUATIONS C M(DY/DT) = F(Y,T) -C SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE +C SUBROUTINE A_OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE C C REFERENCES C @@ -79,8 +79,8 @@ C SPRINGER 1996, page 267. 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 A_OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND +C IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR A_STIFF. C C THE INPUT PARAMETERS ARE .. C N = THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. @@ -165,7 +165,7 @@ 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 COMPUTED IN A_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. @@ -251,20 +251,20 @@ C IN ADDITION TO OVDRIVE, THE FOLLOWING ROUTINES ARE PROVIDED C IN THE PACKAGE.. C -C INTERP( - ) INTERPOLATES TO GET THE OUTPUT VALUES +C A_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 A_STIFF( - ) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A C SINGLE STEP AND ASSOCIATED ERROR CONTROL. -C COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION +C A_COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION C SCHEMES FOR USE IN THE CORE INTEGRATOR. -C PSET( - ) COMPUTES AND PROCESSES THE JACOBIAN +C A_PSET( - ) COMPUTES AND PROCESSES THE JACOBIAN C MATRIX J = DF/DY -C DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. -C SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER DEC +C A_DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. +C A_SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER A_DEC C HAS BEEN CALLED FOR THE MATRIX A -C DGBFA ( - ) FACTORS A DOUBLE PRECISION BAND MATRIX BY +C A_DGBFA ( - ) FACTORS A DOUBLE PRECISION BAND MATRIX BY C ELIMINATION. -C DGBSL ( - ) SOLVES A BANDED LINEAR SYSTEM A*x=b +C A_DGBSL ( - ) SOLVES A BANDED LINEAR SYSTEM A*x=b C C ALSO SUPPLIED ARE THE BLAS ROUTINES C @@ -338,7 +338,7 @@ 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 FOR USE INSIDE THE INTEGRATOR A_STIFF * C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C .. SCALAR ARGUMENTS .. @@ -353,7 +353,7 @@ C COMMON BLOCKS C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL OVDRIV,F,PDERV,MAS + EXTERNAL A_OVDRIV,F,PDERV,MAS C .. C .. SAVE STATEMENT .. SAVE I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12 @@ -401,7 +401,7 @@ c WORKSPACE HAS TO BE AT LEAST N+14. c - CALL OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,WORK(3),WORK(I1), + CALL A_OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,WORK(3),WORK(I1), + WORK(I2),WORK(I3),WORK(I4),WORK(I5),WORK(I6),WORK(I7), + WORK(I8),WORK(I9),WORK(I10),WORK(I11),IWORK(15),MBND,MASBND, + IWORK(1),IWORK(2),IWORK(3),MAXDER,ITOL,RTOL,ATOL,RPAR,IPAR, @@ -431,7 +431,7 @@ + ' WITH N = ',I6) END - SUBROUTINE OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,Y,YHOLD, + SUBROUTINE A_OVDRIV(N,T0,HO,Y0,TOUT,TEND,MF,IDID,LOUT,Y,YHOLD, + YNHOLD,YMAX,ERRORS,SAVE1,SAVE2,SCALE,ARH,PW,PWCOPY, + AM,IPIV,MBND,MASBND,NIND1,NIND2,NIND3,MAXDER,ITOL, + RTOL,ATOL,RPAR,IPAR,F,PDERV,MAS,NQUSED,NSTEP,NFAIL, @@ -457,7 +457,7 @@ INTEGER I,KGO,NHCUT C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL INTERP,STIFF,F,PDERV,MAS + EXTERNAL A_INTERP,A_STIFF,F,PDERV,MAS C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1 @@ -475,7 +475,7 @@ 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) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG T0 = TOUT HO = H @@ -493,7 +493,7 @@ 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) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT HO = H IDID = KFLAG @@ -520,7 +520,7 @@ 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) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) HO = H T0 = TOUT IDID = -5 @@ -534,7 +534,7 @@ 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) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG T0 = TOUT HO = H @@ -667,7 +667,7 @@ 20 IF ((T+H).EQ.T) THEN WRITE (LOUT,9000) END IF - CALL STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND,MASBND, + CALL A_STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND,MASBND, + NIND1,NIND2,NIND3,T,TOUT,TEND,Y,N, + YMAX,ERRORS,SAVE1,SAVE2,SCALE,PW,PWCOPY,AM,YHOLD, + YNHOLD,ARH,IPIV,LOUT,MAXDER,ITOL,RTOL,ATOL,RPAR,IPAR,F, @@ -679,7 +679,7 @@ ENDIF KGO = 1 - KFLAG IF (KGO.EQ.1) THEN -C NORMAL RETURN FROM STIFF +C NORMAL RETURN FROM A_STIFF GO TO 30 ELSE IF (KGO.EQ.2) THEN @@ -756,7 +756,7 @@ 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) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT HO = H IDID = KFLAG @@ -773,7 +773,7 @@ ELSE IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT, SO INTERPOLATE - CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG HO = H T0 = TOUT @@ -812,14 +812,14 @@ ELSE C HAVE PASSED TOUT SO INTERPOLATE - CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) + CALL A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT IDID = KFLAG END IF HO = H IF(KFLAG.NE.0) IDID = KFLAG RETURN -C -------------------------- END OF SUBROUTINE OVDRIV ----------------- +C -------------------------- END OF SUBROUTINE A_OVDRIV ----------------- 9000 FORMAT (' WARNING.. T + H = T ON NEXT STEP.') 9010 FORMAT (/,/,' KFLAG = -2 FROM INTEGRATOR AT T = ',E16.8,' H =', + E16.8,/, @@ -853,7 +853,7 @@ END - SUBROUTINE INTERP(N,JSTART,H,T,Y,TOUT,Y0) + SUBROUTINE A_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. @@ -880,14 +880,14 @@ 20 CONTINUE 30 CONTINUE RETURN -C -------------- END OF SUBROUTINE INTERP --------------------------- +C -------------- END OF SUBROUTINE A_INTERP --------------------------- END - SUBROUTINE COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + SUBROUTINE A_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 A_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 @@ -1017,10 +1017,10 @@ 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 --------------------- +C --------------------- END OF SUBROUTINE A_COSET --------------------- END - SUBROUTINE PSET(Y,N,H,T,UROUND,EPSJAC,CON,MITER,MBND, + SUBROUTINE A_PSET(Y,N,H,T,UROUND,EPSJAC,CON,MITER,MBND, + MASBND,NIND1,NIND2,NIND3,IER,F,PDERV,MAS, + NRENEW,YMAX,SAVE1,SAVE2,PW,PWCOPY,AM,WRKSPC,IPIV, + ITOL,RTOL,ATOL,NPSET,NJE,NFE,NDEC,IPAR,RPAR,IERR) @@ -1029,7 +1029,7 @@ IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------- -C PSET IS CALLED BY STIFF TO COMPUTE AND PROCESS THE MATRIX +C A_PSET IS CALLED BY A_STIFF TO COMPUTE AND PROCESS THE MATRIX C M/(H*EL(1)) - J WHERE J IS AN APPROXIMATION TO THE RELEVANT JACOBIAN C AND M IS THE MASS MATRIX. THIS MATRIX IS THEN SUBJECTED TO LU C DECOMPOSITION IN PREPARATION FOR LATER SOLUTION OF LINEAR SYSTEMS @@ -1037,7 +1037,7 @@ C MATRIX J IS FOUND BY THE USER-SUPPLIED ROUTINE PDERV IF MITER=1 C OR 3 OR BY FINITE DIFFERENCING IF MITER = 2 OR 4. C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH -C PSET USES THE FOLLOWING .. +C A_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 @@ -1060,7 +1060,7 @@ INTEGER I,J,J1,JJKK,FOUR,FIVE C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL DEC,F,PDERV,DGBFA,MAS + EXTERNAL A_DEC,F,PDERV,A_DGBFA,MAS C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSQRT @@ -1267,7 +1267,7 @@ II = II + MBND(4) 75 CONTINUE ENDIF - CALL DGBFA(PW,MBND(4),N,ML,MU,IPIV,IER) + CALL A_DGBFA(PW,MBND(4),N,ML,MU,IPIV,IER) NDEC = NDEC + 1 ELSE IF(MASBND(1).EQ.0) THEN @@ -1278,13 +1278,13 @@ J = J + NP1 80 CONTINUE ENDIF - CALL DEC(N,N,PW,IPIV,IER) + CALL A_DEC(N,N,PW,IPIV,IER) NDEC = NDEC + 1 ENDIF RETURN -C ---------------------- END OF SUBROUTINE PSET --------------------- +C ---------------------- END OF SUBROUTINE A_PSET --------------------- END - SUBROUTINE DEC(N,NDIM,A,IP,IER) + SUBROUTINE A_DEC(N,NDIM,A,IP,IER) IMPLICIT DOUBLE PRECISION(A-H,O-Z) @@ -1301,9 +1301,9 @@ C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0. C IER = 0 IF MATRIX IS NON-SINGULAR, OR K IF FOUND TO BE SINGULAR C AT STAGE K. -C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. +C USE A_SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. C DETERM(A) = IP(N)*A(1,1)*A(2,2)* . . . *A(N,N). -C IF IP(N) = 0, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. +C IF IP(N) = 0, A IS SINGULAR, A_SOL WILL DIVIDE BY ZERO. C C REFERENCE. C C.B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, C.A.C.M @@ -1362,9 +1362,9 @@ 80 IER = K IP(N) = 0 RETURN -C --------------------- END OF SUBROUTINE DEC ---------------------- +C --------------------- END OF SUBROUTINE A_DEC ---------------------- END - SUBROUTINE SOL(N,NDIM,A,B,IP) + SUBROUTINE A_SOL(N,NDIM,A,B,IP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) @@ -1386,10 +1386,10 @@ C INPUT .. C N = ORDER OF MATRIX. C NDIM = DECLARED DIMENSION OF MATRIX A. -C A = TRIANGULARISED MATRIX OBTAINED FROM DEC. +C A = TRIANGULARISED MATRIX OBTAINED FROM A_DEC. C B = RIGHT HAND SIDE VECTOR. -C IP = PIVOT VECTOR OBTAINED FROM DEC. -C DO NOT USE IF DEC HAS SET IER .NE. 0 +C IP = PIVOT VECTOR OBTAINED FROM A_DEC. +C DO NOT USE IF A_DEC HAS SET IER .NE. 0 C OUTPUT.. C B = SOLUTION VECTOR, X. C ------------------------------------------------------------------ @@ -1416,15 +1416,15 @@ 40 CONTINUE 50 B(1) = B(1)/A(1,1) RETURN -C ------------------------- END OF SUBROUTINE SOL ------------------ +C ------------------------- END OF SUBROUTINE A_SOL ------------------ END - subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info) + subroutine a_dgbfa(abd,lda,n,ml,mu,ipvt,info) integer lda,n,ml,mu,ipvt(1),info double precision abd(lda,1) c -c dgbfa factors a double precision band matrix by elimination. +c a_dgbfa factors a double precision band matrix by elimination. c -c dgbfa is usually called by dgbco, but it can be called +c a_dgbfa is usually called by dgbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry @@ -1466,7 +1466,7 @@ c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does -c indicate that dgbsl will divide by zero if +c indicate that a_dgbsl will divide by zero if c called. use rcond in dgbco for a reliable c indication of singularity. c @@ -1593,151 +1593,18 @@ return end c - subroutine daxpy(n,da,dx,incx,dy,incy) -c -c constant times a vector plus a vector. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1),da - integer i,incx,incy,ix,iy,m,mp1,n -c - if(n.le.0)return - if (da .eq. 0.0d0) return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dy(iy) = dy(iy) + da*dx(ix) - ix = ix + incx - iy = iy + incy - 10 continue - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,4) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dy(i) = dy(i) + da*dx(i) - 30 continue - if( n .lt. 4 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,4 - dy(i) = dy(i) + da*dx(i) - dy(i + 1) = dy(i + 1) + da*dx(i + 1) - dy(i + 2) = dy(i + 2) + da*dx(i + 2) - dy(i + 3) = dy(i + 3) + da*dx(i + 3) - 50 continue - return - end -c - subroutine dscal(n,da,dx,incx) -c -c scales a vector by a constant. -c uses unrolled loops for increment equal to one. -c jack dongarra, linpack, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. -c - double precision da,dx(1) - integer i,incx,ix,m,mp1,n -c - if(n.le.0)return - if(incx.eq.1)go to 20 -c -c code for increment not equal to 1 -c - ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - do 10 i = 1,n - dx(ix) = da*dx(ix) - ix = ix + incx - 10 continue - return -c -c code for increment equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dx(i) = da*dx(i) - 30 continue - if( n .lt. 5 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - dx(i) = da*dx(i) - dx(i + 1) = da*dx(i + 1) - dx(i + 2) = da*dx(i + 2) - dx(i + 3) = da*dx(i + 3) - dx(i + 4) = da*dx(i + 4) - 50 continue - return - end -c - integer function idamax(n,dx,incx) -c -c finds the index of element having max. absolute value. -c jack dongarra, linpack, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. -c - double precision dx(1),dmax - integer i,incx,ix,n -c - idamax = 0 - if( n .lt. 1 ) return - idamax = 1 - if(n.eq.1)return - if(incx.eq.1)go to 20 -c -c code for increment not equal to 1 -c - ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - dmax = dabs(dx(ix)) - ix = ix + incx - do 10 i = 2,n - if(dabs(dx(ix)).le.dmax) go to 5 - idamax = i - dmax = dabs(dx(ix)) - 5 ix = ix + incx - 10 continue - return -c -c code for increment equal to 1 -c - 20 dmax = dabs(dx(1)) - do 30 i = 2,n - if(dabs(dx(i)).le.dmax) go to 30 - idamax = i - dmax = dabs(dx(i)) - 30 continue - return - end - - subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job) + subroutine a_dgbsl(abd,lda,n,ml,mu,ipvt,b,job) integer lda,n,ml,mu,ipvt(*),job double precision abd(lda,*),b(*) c -c dgbsl solves the double precision band system +c a_dgbsl solves the double precision band system c a * x = b or trans(a) * x = b -c using the factors computed by dgbco or dgbfa. +c using the factors computed by dgbco or a_dgbfa. c c on entry c c abd double precision(lda, n) -c the output from dgbco or dgbfa. +c the output from dgbco or a_dgbfa. c c lda integer c the leading dimension of the array abd . @@ -1752,7 +1619,7 @@ c number of diagonals above the main diagonal. c c ipvt integer(n) -c the pivot vector from dgbco or dgbfa. +c the pivot vector from dgbco or a_dgbfa. c c b double precision(n) c the right hand side vector. @@ -1773,14 +1640,14 @@ c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgbco has set rcond .gt. 0.0 -c or dgbfa has set info .eq. 0 . +c or a_dgbfa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p -c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) +c call a_dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . @@ -1862,62 +1729,13 @@ return end c - double precision function ddot(n,dx,incx,dy,incy) -c -c forms the dot product of two vectors. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1),dtemp - integer i,incx,incy,ix,iy,m,mp1,n -c - ddot = 0.0d0 - dtemp = 0.0d0 - if(n.le.0)return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dtemp = dtemp + dx(ix)*dy(iy) - ix = ix + incx - iy = iy + incy - 10 continue - ddot = dtemp - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dtemp = dtemp + dx(i)*dy(i) - 30 continue - if( n .lt. 5 ) go to 60 - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + - * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) - 50 continue - 60 ddot = dtemp - return - end - - SUBROUTINE ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) + SUBROUTINE A_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 IN A_STIFF . C C *************************************************** C .. SCALAR ARGUMENTS .. @@ -1950,7 +1768,7 @@ C ** ERROR ASSOCIATED WITH METHOD OF ORDER TWO LOWER. RETURN END - SUBROUTINE PRDICT(T,H,Y,L,N,YPRIME,NFE,IPAR,RPAR,F,IERR) + SUBROUTINE A_PRDICT(T,H,Y,L,N,YPRIME,NFE,IPAR,RPAR,F,IERR) @@ -1987,10 +1805,10 @@ RETURN END - SUBROUTINE ITRAT2(QQQ,Y,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE,M,WORKED - + ,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND,AM,MASBND,NIND1, - + NIND2,NIND3,IPIV,LMB,ITOL,RTOL,ATOL,IPAR,RPAR,HUSED,NBSOL, - + NFE,NQUSED,F,IERR) + SUBROUTINE A_ITRAT2(QQQ,Y,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE,M, + + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND,AM,MASBND, + + NIND1,NIND2,NIND3,IPIV,LMB,ITOL,RTOL,ATOL,IPAR,RPAR,HUSED, + + NBSOL,NFE,NQUSED,F,IERR) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. @@ -2006,7 +1824,7 @@ INTEGER I C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL F,SOL,DGBSL + EXTERNAL F,A_SOL,A_DGBSL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DMAX1,DMIN1 @@ -2077,10 +1895,10 @@ 8812 CONTINUE ENDIF IF(MF.GE.23) THEN - CALL DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) + CALL A_DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) NBSOL = NBSOL + 1 ELSE - CALL SOL(N,N,PW,SAVE1,IPIV) + CALL A_SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 ENDIF D = ZERO @@ -2131,10 +1949,10 @@ C IF WE ARE HERE THEN PARTIALS ARE O.K. C IF( MF.GE. 23) THEN - CALL DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) + CALL A_DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) NBSOL=NBSOL + 1 ELSE - CALL SOL(N,N,PW,SAVE1,IPIV) + CALL A_SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 ENDIF C @@ -2180,7 +1998,7 @@ END - SUBROUTINE STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND, + SUBROUTINE A_STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND, + MASBND,NIND1,NIND2,NIND3,T,TOUT,TEND,Y,N, + YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,PWCOPY,AM,YHOLD, + YNHOLD,ARH,IPIV,LOUT,MAXDER,ITOL,RTOL,ATOL,RPAR,IPAR,F, @@ -2191,13 +2009,13 @@ IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------ -C THE SUBROUTINE STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN +C THE SUBROUTINE A_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 COMMUNICATION WITH A_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 MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE A_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. @@ -2207,7 +2025,7 @@ 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 RTOL,ATOL THE ERROR BOUNDS. SEE DESCRIPTION IN A_OVDRIV. 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.. @@ -2242,7 +2060,7 @@ 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 NRENEW FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE A_PSET. C IF NRENEW > 0 THEN FORM A NEW JACOBIAN BEFORE C COMPUTING THE COEFFICIENT MATRIX FOR C THE NEWTON-RAPHSON ITERATION @@ -2271,8 +2089,8 @@ DIMENSION EL(10),ELST(10),TQ(5) C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL COSET,CPYARY,ERRORS,F,HCHOSE,ITRAT2, - + PRDICT,PSET,RSCALE,SOL,DGBSL,PDERV,MAS + EXTERNAL A_COSET,A_CPYARY,A_ERRORS,F,A_HCHOSE,A_ITRAT2, + + A_PRDICT,A_PSET,A_RSCALE,A_SOL,A_DGBSL,PDERV,MAS C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DMIN1 @@ -2378,7 +2196,7 @@ 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) + CALL A_COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) @@ -2389,20 +2207,20 @@ 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) + CALL A_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) + CALL A_CPYARY(N*L,Y,YHOLD) QI = H*EL(1) QQ = ONE/QI - CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,IPAR,RPAR,F,IERR) + CALL A_PRDICT(T,H,Y,L,N,SAVE2,NFE,IPAR,RPAR,F,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) + 30 CALL A_CPYARY(N*L,YHOLD,Y) IF (MF.NE.MFOLD) THEN METH = MF/10 MITER = MF - 10*METH @@ -2445,7 +2263,7 @@ C ********************************************* 40 RH = DMAX1(RH,HMIN/DABS(H)) 50 RH = DMIN1(RH,HMAX/DABS(H),RMAX) - CALL RSCALE(N,L,RH,Y) + CALL A_RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH @@ -2462,7 +2280,7 @@ END IF IDOUB = L + 1 - CALL CPYARY(N*L,Y,YHOLD) + CALL A_CPYARY(N*L,Y,YHOLD) 60 IF (DABS(RC-ONE).GT.UPBND) IWEVAL = MITER HUSED = H @@ -2487,7 +2305,7 @@ 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,IPAR,RPAR,F,IERR) + CALL A_PRDICT(T,H,Y,L,N,SAVE2,NFE,IPAR,RPAR,F,IERR) IF(IERR.NE.0) GOTO 8000 ELSE @@ -2507,7 +2325,7 @@ 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 THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN A_PSET. C THE PROCESSED MATRIX IS STORED IN PW C ------------------------------------------------------------------- IWEVAL = 0 @@ -2573,13 +2391,13 @@ JSNOLD = 0 MQ1TMP = MEQC1 MQ2TMP = MEQC2 - CALL PSET(Y,N,H,T,UROUND,EPSJAC,QI,MITER,MBND,MASBND, + CALL A_PSET(Y,N,H,T,UROUND,EPSJAC,QI,MITER,MBND,MASBND, + NIND1,NIND2,NIND3,IER,F,PDERV,MAS,NRENEW,YMAX,SAVE1,SAVE2, + PW,PWCOPY,AM,ERROR,IPIV,ITOL,RTOL,ATOL,NPSET,NJE,NFE,NDEC,IPAR + ,RPAR,IERR) IF(IERR.NE.0) GOTO 8000 QQQ=QI -C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET +C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY A_PSET IF (IER.NE.0) THEN C IF IER>0 THEN WE HAVE HAD A SINGULARITY IN THE ITERATION MATRIX IJUS=1 @@ -2603,7 +2421,7 @@ 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, + CALL A_ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1,WORKED,YMAX, + ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND,AM,MASBND, + NIND1,NIND2,NIND3,IPIV,1,ITOL,RTOL,ATOL,IPAR,RPAR,HUSED,NBSOL, + NFE,NQUSED,F,IERR) @@ -2611,7 +2429,7 @@ ITST = 2 ELSE - CALL ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1,WORKED,YMAX, + CALL A_ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1,WORKED,YMAX, + ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND,AM,MASBND, +NIND1,NIND2,NIND3,IPIV,0,ITOL,RTOL,ATOL,IPAR,RPAR,HUSED,NBSOL, + NFE,NQUSED,F,IERR) @@ -2752,7 +2570,7 @@ ARH(I) = ARH(I) + EL(JP1)*Y(I,J1) 200 CONTINUE 210 CONTINUE - CALL PRDICT(T,H,Y,L,N,SAVE2,NFE,IPAR,RPAR,F,IERR) + CALL A_PRDICT(T,H,Y,L,N,SAVE2,NFE,IPAR,RPAR,F,IERR) IF(IERR.NE.0) GOTO 8000 DO 220 I = 1,N SAVE1(I) = Y(I,1) @@ -2763,7 +2581,7 @@ 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, + CALL A_ITRAT2(QQQ,Y,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2,WORKED,YMAX, + ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND,AM,MASBND, +NIND1,NIND2,NIND3,IPIV,1,ITOL,RTOL,ATOL,IPAR,RPAR,HUSED,NBSOL, + NFE,NQUSED,F,IERR) @@ -2872,10 +2690,10 @@ 3111 CONTINUE ENDIF IF (MF.GE. 23) THEN - CALL DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) + CALL A_DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) NBSOL=NBSOL+1 ELSE - CALL SOL(N,N,PW,SAVE1,IPIV) + CALL A_SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 ENDIF DO 321 I=1,N @@ -2971,7 +2789,7 @@ 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) + CALL A_CPYARY(N*L,YHOLD,Y) RMAX = 2.0D+0 IF (DABS(H).LE.HMIN*1.00001D+0) THEN C @@ -3000,10 +2818,10 @@ NQ=NEWQ RH=ONE/(PLFAIL*DBLE(FLOAT(-KFAIL))) L=NQ+1 - CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + CALL A_COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) RC=RC*EL(1)/OLDLO OLDLO=EL(1) - CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) + CALL A_ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) ELSE NEWQ = NQ RH = ONE/ (PRFAIL*DBLE(FLOAT(-KFAIL))) @@ -3029,7 +2847,7 @@ C ********************************* JCHANG = 1 RH = DMAX1(HMIN/DABS(H),0.1D+0) - CALL HCHOSE(RH,H,OVRIDE) + CALL A_HCHOSE(RH,H,OVRIDE) H = H*RH CALL F(N,T,YHOLD,SAVE1,IPAR,RPAR,IERR) IF(IERR.NE.0) GOTO 8000 @@ -3048,11 +2866,11 @@ NQ = 1 L = 2 C RESET ORDER, RECALCULATE ERROR BOUNDS - CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + CALL A_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) + CALL A_ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) C NOW JUMP TO NORMAL CONTINUATION POINT GO TO 60 C ********************************************************************** @@ -3216,7 +3034,7 @@ GOTO 440 ENDIF RH = DMIN1(RH,RMAX) - CALL HCHOSE(RH,H,OVRIDE) + CALL A_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 @@ -3244,17 +3062,17 @@ NQ = NEWQ L = NQ + 1 C RESET ORDER,RECALCULATE ERROR BOUNDS - CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + CALL A_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) + CALL A_ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) END IF RH = DMAX1(RH,HMIN/DABS(H)) RH = DMIN1(RH,HMAX/DABS(H),RMAX) - CALL RSCALE(N,L,RH,Y) + CALL A_RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH @@ -3271,7 +3089,7 @@ 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) + CALL A_CPYARY(N*L,Y,YHOLD) NSTEP = NSTEP + 1 JSINUP = JSINUP + 1 JSNOLD = JSNOLD + 1 @@ -3312,17 +3130,17 @@ C TRY AGAIN WITH UPDATED PARTIALS C 8000 IF(IERR.NE.0) RETURN - IF(IJUS.EQ.0) CALL HCHOSE(RH,H,OVRIDE) + IF(IJUS.EQ.0) CALL A_HCHOSE(RH,H,OVRIDE) IF(.NOT.FINISH) THEN GO TO 40 ELSE RETURN END IF -C ------------------- END OF SUBROUTINE STIFF -------------------------- +C ------------------- END OF SUBROUTINE A_STIFF -------------------------- 9000 FORMAT (1X,' CORRECTOR HAS NOT CONVERGED') END - SUBROUTINE RSCALE(N,L,RH,Y) + SUBROUTINE A_RSCALE(N,L,RH,Y) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. @@ -3432,7 +3250,7 @@ RETURN END - SUBROUTINE CPYARY(NELEM,SOURCE,TARGET) + SUBROUTINE A_CPYARY(NELEM,SOURCE,TARGET) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C COPIES THE ARRAY SOURCE() INTO THE ARRAY TARGET() @@ -3455,7 +3273,7 @@ RETURN END - SUBROUTINE HCHOSE(RH,H,OVRIDE) + SUBROUTINE A_HCHOSE(RH,H,OVRIDE) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON / STPSZE / HSTPSZ(2,14) LOGICAL OVRIDE @@ -3492,947 +3310,3 @@ C ************************************************************ C END - 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 --- cash.orig/mebdfi.f 2007-12-15 15:37:46.000000000 -0500 +++ cash/mebdfi.f 2014-03-02 16:22:33.208828923 -0500 @@ -58,11 +58,11 @@ C C SEPTEMBER 20th 1999: FIRST RELEASE C -C OVDRIV +C I_OVDRIV C A PACKAGE FOR THE SOLUTION OF THE INITIAL VALUE PROBLEM C FOR SYSTEMS OF IMPLICIT DIFFERENTIAL ALGEBRAIC EQUATIONS c G(t,Y,Y')=0, Y=(Y(1),Y(2),Y(3),.....,Y(N)). -C SUBROUTINE OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE. +C SUBROUTINE I_OVDRIV IS A DRIVER ROUTINE FOR THIS PACKAGE. C C REFERENCES C @@ -82,7 +82,7 @@ C SPRINGER 1996, page 267. C C ---------------------------------------------------------------- -C OVDRIV IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND +C I_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 .. @@ -158,7 +158,7 @@ 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 COMPUTED IN I_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. @@ -242,23 +242,23 @@ C -12 INSUFFICIENT INTEGER WORKSPACE FOR THE INTEGRATION C C -C IN ADDITION TO OVDRIVE, THE FOLLOWING ROUTINES ARE PROVIDED +C IN ADDITION TO I_OVDRIVE, THE FOLLOWING ROUTINES ARE PROVIDED C IN THE PACKAGE.. C -C INTERP( - ) INTERPOLATES TO GET THE OUTPUT VALUES +C I_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 I_STIFF( - ) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A C SINGLE STEP AND ASSOCIATED ERROR CONTROL. -C COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION +C I_COSET( - ) SETS COEFFICIENTS FOR BACKWARD DIFFERENTIATION C SCHEMES FOR USE IN THE CORE INTEGRATOR. -C PSET( - ) COMPUTES AND PROCESSES THE NEWTON ITERATION +C I_PSET( - ) COMPUTES AND PROCESSES THE NEWTON ITERATION C MATRIX DG/DY + (1/(H*BETA))DG/DY' -C DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. -C SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER DEC +C I_DEC( - ) PERFORMS AN LU DECOMPOSITION ON A MATRIX. +C I_SOL( - ) SOLVES LINEAR SYSTEMS A*X = B AFTER I_DEC C HAS BEEN CALLED FOR THE MATRIX A -C DGBFA ( - ) FACTORS A DOUBLE PRECISION BAND MATRIX BY +C I_DGBFA ( - ) FACTORS A DOUBLE PRECISION BAND MATRIX BY C ELIMINATION. -C DGBSL ( - ) SOLVES A BANDED LINEAR SYSTEM A*x=b +C I_DGBSL ( - ) SOLVES A BANDED LINEAR SYSTEM A*x=b C C ALSO SUPPLIED ARE THE BLAS ROUTINES C @@ -330,7 +330,7 @@ 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 FOR USE INSIDE THE INTEGRATOR I_STIFF * C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C .. SCALAR ARGUMENTS .. @@ -346,7 +346,7 @@ C COMMON BLOCKS C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL OVDRIV,PDERV,RESID + EXTERNAL I_OVDRIV,PDERV,RESID C .. C .. SAVE STATEMENT .. SAVE I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11 @@ -397,7 +397,7 @@ c THE ERROR FLAG IS INITIALISED c - CALL OVDRIV(N,T0,HO,Y0,YPRIME,TOUT,TEND,MF,IDID,LOUT,WORK(3), + CALL I_OVDRIV(N,T0,HO,Y0,YPRIME,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),IWORK(15), + MBND,IWORK(1),IWORK(2),IWORK(3),MAXDER,ITOL,RTOL,ATOL,RPAR, @@ -428,7 +428,7 @@ END C-------------------------------------------------------------------------- C - SUBROUTINE OVDRIV(N,T0,HO,Y0,YPRIME,TOUT,TEND,MF,IDID,LOUT,Y, + SUBROUTINE I_OVDRIV(N,T0,HO,Y0,YPRIME,TOUT,TEND,MF,IDID,LOUT,Y, + YHOLD,YNHOLD,YMAX,ERRORS,SAVE1,SAVE2,SCALE,ARH,PW,PWCOPY, + IPIV,MBND,NIND1,NIND2,NIND3,MAXDER,ITOL,RTOL,ATOL,RPAR, + IPAR,PDERV,RESID,NQUSED,NSTEP,NFAIL,NRE,NJE,NDEC,NBSOL, @@ -452,7 +452,7 @@ INTEGER I,KGO,NHCUT C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL INTERP,STIFF,PDERV,RESID + EXTERNAL I_INTERP,I_STIFF,PDERV,RESID C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1 @@ -468,7 +468,7 @@ 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) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG T0 = TOUT HO = H @@ -486,7 +486,7 @@ 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) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT HO = H IDID = KFLAG @@ -513,7 +513,7 @@ 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) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) HO = H T0 = TOUT IDID = -5 @@ -527,7 +527,7 @@ 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) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG T0 = TOUT HO = H @@ -660,7 +660,7 @@ 20 IF ((T+H).EQ.T) THEN WRITE (LOUT,9000) END IF - CALL STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND, + CALL I_STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND, + NIND1,NIND2,NIND3,T,TOUT,TEND,Y,YPRIME,N, + YMAX,ERRORS,SAVE1,SAVE2,SCALE,PW,PWCOPY,YHOLD, + YNHOLD,ARH,IPIV,LOUT,MAXDER,ITOL,RTOL,ATOL,RPAR,IPAR, @@ -672,7 +672,7 @@ C ENDIF KGO = 1 - KFLAG IF (KGO.EQ.1) THEN -C NORMAL RETURN FROM STIFF +C NORMAL RETURN FROM I_STIFF GO TO 30 ELSE IF (KGO.EQ.2) THEN @@ -708,7 +708,7 @@ 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 IF INTERPOLATION IS NOT DESIRED, THE CALL TO I_INTERP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520. C -------------------------------------------------------------------- IF(NSTEP.GT.MAXSTP) THEN @@ -749,7 +749,7 @@ 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) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT HO = H IDID = KFLAG @@ -766,7 +766,7 @@ ELSE IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT, SO INTERPOLATE - CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IDID = KFLAG HO = H T0 = TOUT @@ -805,14 +805,14 @@ ELSE C HAVE PASSED TOUT SO INTERPOLATE - CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) T0 = TOUT IDID = KFLAG END IF HO = H IF(KFLAG.NE.0) IDID = KFLAG RETURN -C -------------------------- END OF SUBROUTINE OVDRIV ----------------- +C -------------------------- END OF SUBROUTINE I_OVDRIV ----------------- 9000 FORMAT (' WARNING.. T + H = T ON NEXT STEP.') 9010 FORMAT (/,/,' KFLAG = -2 FROM INTEGRATOR AT T = ',E16.8,' H =', + E16.8,/, @@ -848,7 +848,7 @@ END C-------------------------------------------------------------------------- C - SUBROUTINE INTERP(N,JSTART,H,T,Y,TOUT,Y0) + SUBROUTINE I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. @@ -875,15 +875,15 @@ 20 CONTINUE 30 CONTINUE RETURN -C -------------- END OF SUBROUTINE INTERP --------------------------- +C -------------- END OF SUBROUTINE I_INTERP --------------------------- END C - SUBROUTINE COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + SUBROUTINE I_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 I_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 @@ -1013,24 +1013,24 @@ 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 --------------------- +C --------------------- END OF SUBROUTINE I_COSET --------------------- END - SUBROUTINE PSET(Y,YPRIME,N,H,T,UROUND,EPSJAC,CON,MITER,MBND, + SUBROUTINE I_PSET(Y,YPRIME,N,H,T,UROUND,EPSJAC,CON,MITER,MBND, + NIND1,NIND2,NIND3,IER,PDERV,RESID,NRENEW,YMAX,SAVE1,SAVE2, + SAVE3,PW,PWCOPY,WRKSPC,IPIV,ITOL,RTOL,ATOL,NPSET,NJE,NRE, + NDEC,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_PSET IS CALLED BY I_STIFF TO COMPUTE AND PROCESS THE MATRIX C PD=DG/DY + (1/CON)DG/DY'. 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. THE C MATRIX PD IS FOUND BY THE USER-SUPPLIED ROUTINE PDERV IF MITER=1 C OR 3 OR BY FINITE DIFFERENCING IF MITER = 2 OR 4. C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH -C PSET USES THE FOLLOWING .. +C I_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 @@ -1052,7 +1052,7 @@ INTEGER I,J,J1,JJKK C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL DEC,PDERV,DGBFA,RESID + EXTERNAL I_DEC,PDERV,I_DGBFA,RESID C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DSQRT @@ -1192,17 +1192,17 @@ NRE=NRE+ MIN(MBND(3),N) C 70 IF (MITER.GT.2) THEN - CALL DGBFA(PW,MBND(4),N,ML,MU,IPIV,IER) + CALL I_DGBFA(PW,MBND(4),N,ML,MU,IPIV,IER) NDEC = NDEC + 1 ELSE - CALL DEC(N,N,PW,IPIV,IER) + CALL I_DEC(N,N,PW,IPIV,IER) NDEC = NDEC + 1 ENDIF RETURN -C ---------------------- END OF SUBROUTINE PSET --------------------- +C ---------------------- END OF SUBROUTINE I_PSET --------------------- END C - SUBROUTINE DEC(N,NDIM,A,IP,IER) + SUBROUTINE I_DEC(N,NDIM,A,IP,IER) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------- @@ -1218,9 +1218,9 @@ C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR 0. C IER = 0 IF MATRIX IS NON-SINGULAR, OR K IF FOUND TO BE SINGULAR C AT STAGE K. -C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. +C USE I_SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. C DETERM(A) = IP(N)*A(1,1)*A(2,2)* . . . *A(N,N). -C IF IP(N) = 0, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. +C IF IP(N) = 0, A IS SINGULAR, I_SOL WILL DIVIDE BY ZERO. C C REFERENCE. C C.B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, C.A.C.M @@ -1279,10 +1279,10 @@ 80 IER = K IP(N) = 0 RETURN -C--------------------- END OF SUBROUTINE DEC ---------------------- +C--------------------- END OF SUBROUTINE I_DEC ---------------------- END C - SUBROUTINE SOL(N,NDIM,A,B,IP) + SUBROUTINE I_SOL(N,NDIM,A,B,IP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. @@ -1305,8 +1305,8 @@ C NDIM = DECLARED DIMENSION OF MATRIX A. C A = TRIANGULARISED MATRIX OBTAINED FROM DEC. C B = RIGHT HAND SIDE VECTOR. -C IP = PIVOT VECTOR OBTAINED FROM DEC. -C DO NOT USE IF DEC HAS SET IER .NE. 0 +C IP = PIVOT VECTOR OBTAINED FROM I_DEC. +C DO NOT USE IF I_DEC HAS SET IER .NE. 0 C OUTPUT.. C B = SOLUTION VECTOR, X. C ------------------------------------------------------------------ @@ -1333,16 +1333,16 @@ 40 CONTINUE 50 B(1) = B(1)/A(1,1) RETURN -C------------------------- END OF SUBROUTINE SOL ------------------ +C------------------------- END OF SUBROUTINE I_SOL ------------------ END C - subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info) + subroutine i_dgbfa(abd,lda,n,ml,mu,ipvt,info) integer lda,n,ml,mu,ipvt(1),info double precision abd(lda,1) c -c dgbfa factors a double precision band matrix by elimination. +c i_dgbfa factors a double precision band matrix by elimination. c -c dgbfa is usually called by dgbco, but it can be called +c i_dgbfa is usually called by dgbco, but it can be called c directly with a saving in time if rcond is not needed. c c on entry @@ -1384,7 +1384,7 @@ c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does -c indicate that dgbsl will divide by zero if +c indicate that i_dgbsl will divide by zero if c called. use rcond in dgbco for a reliable c indication of singularity. c @@ -1511,151 +1511,18 @@ return end C-------------------------------------------------------------------------- - subroutine daxpy(n,da,dx,incx,dy,incy) -c -c constant times a vector plus a vector. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1),da - integer i,incx,incy,ix,iy,m,mp1,n -c - if(n.le.0)return - if (da .eq. 0.0d0) return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dy(iy) = dy(iy) + da*dx(ix) - ix = ix + incx - iy = iy + incy - 10 continue - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,4) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dy(i) = dy(i) + da*dx(i) - 30 continue - if( n .lt. 4 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,4 - dy(i) = dy(i) + da*dx(i) - dy(i + 1) = dy(i + 1) + da*dx(i + 1) - dy(i + 2) = dy(i + 2) + da*dx(i + 2) - dy(i + 3) = dy(i + 3) + da*dx(i + 3) - 50 continue - return - end -C--------------------------------------------------------------------------- - subroutine dscal(n,da,dx,incx) -c -c scales a vector by a constant. -c uses unrolled loops for increment equal to one. -c jack dongarra, linpack, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. -c - double precision da,dx(1) - integer i,incx,ix,m,mp1,n -c - if(n.le.0)return - if(incx.eq.1)go to 20 -c -c code for increment not equal to 1 -c - ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - do 10 i = 1,n - dx(ix) = da*dx(ix) - ix = ix + incx - 10 continue - return -c -c code for increment equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dx(i) = da*dx(i) - 30 continue - if( n .lt. 5 ) return - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - dx(i) = da*dx(i) - dx(i + 1) = da*dx(i + 1) - dx(i + 2) = da*dx(i + 2) - dx(i + 3) = da*dx(i + 3) - dx(i + 4) = da*dx(i + 4) - 50 continue - return - end -C-------------------------------------------------------------------------- - integer function idamax(n,dx,incx) -c -c finds the index of element having max. absolute value. -c jack dongarra, linpack, 3/11/78. -c modified to correct problem with negative increment, 8/21/90. -c - double precision dx(1),dmax - integer i,incx,ix,n -c - idamax = 0 - if( n .lt. 1 ) return - idamax = 1 - if(n.eq.1)return - if(incx.eq.1)go to 20 -c -c code for increment not equal to 1 -c - ix = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - dmax = dabs(dx(ix)) - ix = ix + incx - do 10 i = 2,n - if(dabs(dx(ix)).le.dmax) go to 5 - idamax = i - dmax = dabs(dx(ix)) - 5 ix = ix + incx - 10 continue - return -c -c code for increment equal to 1 -c - 20 dmax = dabs(dx(1)) - do 30 i = 2,n - if(dabs(dx(i)).le.dmax) go to 30 - idamax = i - dmax = dabs(dx(i)) - 30 continue - return - end -C-------------------------------------------------------------------------- - subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job) + subroutine i_dgbsl(abd,lda,n,ml,mu,ipvt,b,job) integer lda,n,ml,mu,ipvt(*),job double precision abd(lda,*),b(*) c -c dgbsl solves the double precision band system +c i_dgbsl solves the double precision band system c a * x = b or trans(a) * x = b -c using the factors computed by dgbco or dgbfa. +c using the factors computed by dgbco or i_dgbfa. c c on entry c c abd double precision(lda, n) -c the output from dgbco or dgbfa. +c the output from dgbco or i_dgbfa. c c lda integer c the leading dimension of the array abd . @@ -1670,7 +1537,7 @@ c number of diagonals above the main diagonal. c c ipvt integer(n) -c the pivot vector from dgbco or dgbfa. +c the pivot vector from dgbco or i_dgbfa. c c b double precision(n) c the right hand side vector. @@ -1691,14 +1558,14 @@ c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgbco has set rcond .gt. 0.0 -c or dgbfa has set info .eq. 0 . +c or i_dgbfa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p -c call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) +c call i_dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . @@ -1780,64 +1647,14 @@ return end C--------------------------------------------------------------------------- - double precision function ddot(n,dx,incx,dy,incy) -c -c forms the dot product of two vectors. -c uses unrolled loops for increments equal to one. -c jack dongarra, linpack, 3/11/78. -c - double precision dx(1),dy(1),dtemp - integer i,incx,incy,ix,iy,m,mp1,n -c - ddot = 0.0d0 - dtemp = 0.0d0 - if(n.le.0)return - if(incx.eq.1.and.incy.eq.1)go to 20 -c -c code for unequal increments or equal increments -c not equal to 1 -c - ix = 1 - iy = 1 - if(incx.lt.0)ix = (-n+1)*incx + 1 - if(incy.lt.0)iy = (-n+1)*incy + 1 - do 10 i = 1,n - dtemp = dtemp + dx(ix)*dy(iy) - ix = ix + incx - iy = iy + incy - 10 continue - ddot = dtemp - return -c -c code for both increments equal to 1 -c -c -c clean-up loop -c - 20 m = mod(n,5) - if( m .eq. 0 ) go to 40 - do 30 i = 1,m - dtemp = dtemp + dx(i)*dy(i) - 30 continue - if( n .lt. 5 ) go to 60 - 40 mp1 = m + 1 - do 50 i = mp1,n,5 - dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + - * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + - * dx(i + 4)*dy(i + 4) - 50 continue - 60 ddot = dtemp - return - end -C--------------------------------------------------------------------------- - SUBROUTINE ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) + SUBROUTINE I_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 IN I_STIFF . C C *************************************************** C .. SCALAR ARGUMENTS .. @@ -1872,7 +1689,7 @@ END C-------------------------------------------------------------------------- - SUBROUTINE PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) + SUBROUTINE I_PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ********************************************************************** @@ -1903,7 +1720,7 @@ END C------------------------------------------------------------------------ - SUBROUTINE ITRAT2(QQQ,Y,YPRIME,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE + SUBROUTINE I_ITRAT2(QQQ,Y,YPRIME,N,T,HBETA,ERRBND,ARH,CRATE,TCRATE + ,M,WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND,NIND1, + NIND2,NIND3,IPIV,LMB,ITOL,RTOL,ATOL,IPAR,RPAR,HUSED,NBSOL, + NRE,NQUSED,RESID,IERR) @@ -1922,7 +1739,7 @@ INTEGER I C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL SOL,DGBSL,RESID + EXTERNAL I_SOL,I_DGBSL,RESID C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DMAX1,DMIN1 @@ -1963,10 +1780,10 @@ C call resid(n,t,y,save2,yprime,ipar,rpar,ierr) IF(MF.GE.23) THEN - CALL DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE2,0) + CALL I_DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE2,0) NBSOL = NBSOL + 1 ELSE - CALL SOL(N,N,PW,SAVE2,IPIV) + CALL I_SOL(N,N,PW,SAVE2,IPIV) NBSOL = NBSOL + 1 ENDIF D = ZERO @@ -1992,10 +1809,10 @@ C IF WE ARE HERE THEN PARTIALS ARE O.K. C IF( MF.GE. 23) THEN - CALL DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE2,0) + CALL I_DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE2,0) NBSOL=NBSOL + 1 ELSE - CALL SOL(N,N,PW,SAVE2,IPIV) + CALL I_SOL(N,N,PW,SAVE2,IPIV) NBSOL = NBSOL + 1 ENDIF C @@ -2043,7 +1860,7 @@ END C-------------------------------------------------------------------------- - SUBROUTINE STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND, + SUBROUTINE I_STIFF(H,HMAX,HMIN,JSTART,KFLAG,MF,MBND, + NIND1,NIND2,NIND3,T,TOUT,TEND,Y,YPRIME,N, + YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,PWCOPY,YHOLD, + YNHOLD,ARH,IPIV,LOUT,MAXDER,ITOL,RTOL,ATOL,RPAR,IPAR, @@ -2052,13 +1869,13 @@ IMPLICIT DOUBLE PRECISION(A-H,O-Z) C ------------------------------------------------------------------ -C THE SUBROUTINE STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN +C THE SUBROUTINE I_STIFF PERFORMS ONE STEP OF THE INTEGRATION OF AN C INITIAL VALUE PROBLEM FOR A SYSTEM OF C IMPLICIT DIFFERENTIAL ALGEBRAIC EQUATIONS. -C COMMUNICATION WITH STIFF IS DONE WITH THE FOLLOWING VARIABLES.. +C COMMUNICATION WITH I_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 MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE I_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. @@ -2068,7 +1885,7 @@ 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 RTOL,ATOL THE ERROR BOUNDS. SEE DESCRIPTION IN I_OVDRIV. 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.. @@ -2103,7 +1920,7 @@ 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 NRENEW FLAG THAT IS USED IN COMMUNICATION WITH SUBROUTINE I_PSET. C IF NRENEW > 0 THEN FORM A NEW JACOBIAN BEFORE C COMPUTING THE COEFFICIENT MATRIX FOR C THE NEWTON-RAPHSON ITERATION @@ -2130,10 +1947,11 @@ C .. C .. LOCAL ARRAYS .. DIMENSION EL(10),ELST(10),TQ(5) + DIMENSION Y0(N) C .. C .. EXTERNAL SUBROUTINES .. - EXTERNAL COSET,CPYARY,ERRORS,HCHOSE,ITRAT2, - + PRDICT,PSET,RSCALE,SOL,DGBSL,PDERV,RESID + EXTERNAL I_COSET,I_CPYARY,I_ERRORS,I_HCHOSE,I_ITRAT2, + + I_PRDICT,I_PSET,I_RSCALE,I_SOL,I_DGBSL,PDERV,RESID C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC DABS,DMAX1,DMIN1 @@ -2225,14 +2043,14 @@ 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 AND BND MUST BE RESET. E IS A COMPARISON FOR I_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) + CALL I_COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDLO OLDLO = EL(1) @@ -2243,14 +2061,14 @@ 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) + CALL I_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) + CALL I_CPYARY(N*L,Y,YHOLD) QI = H*EL(1) QQ = ONE/QI - CALL PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) + CALL I_PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) IF(IERR.NE.0) THEN H=H/2 IERR = 0 @@ -2263,7 +2081,7 @@ C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> C DIFFERENT PARAMETERS ON THIS CALL < C <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - 30 CALL CPYARY(N*L,YHOLD,Y) + 30 CALL I_CPYARY(N*L,YHOLD,Y) IF (MF.NE.MFOLD) THEN METH = MF/10 MITER = MF - 10*METH @@ -2306,7 +2124,7 @@ C ********************************************* 40 RH = DMAX1(RH,HMIN/DABS(H)) 50 RH = DMIN1(RH,HMAX/DABS(H),RMAX) - CALL RSCALE(N,L,RH,Y) + CALL I_RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH @@ -2323,7 +2141,7 @@ END IF IDOUB = L + 1 - CALL CPYARY(N*L,Y,YHOLD) + CALL I_CPYARY(N*L,Y,YHOLD) 60 IF (DABS(RC-ONE).GT.UPBND) IWEVAL = MITER HUSED = H @@ -2348,7 +2166,7 @@ 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,IPAR,RPAR,IERR) + CALL I_PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) IF(IERR.NE.0) GOTO 8000 DO 95 I=1,N YPRIME(I)=(Y(I,1)-ARH(I))/QI @@ -2371,7 +2189,7 @@ 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 THAT THIS HAS BEEN DONE. P IS COMPUTED AND PROCESSED IN I_PSET. C THE PROCESSED MATRIX IS STORED IN PW C ------------------------------------------------------------------- IWEVAL = 0 @@ -2436,14 +2254,14 @@ JSNOLD = 0 MQ1TMP = MEQC1 MQ2TMP = MEQC2 - CALL PSET(Y,YPRIME,N,H,T,UROUND,EPSJAC,QI,MITER,MBND, + CALL I_PSET(Y,YPRIME,N,H,T,UROUND,EPSJAC,QI,MITER,MBND, + NIND1,NIND2,NIND3,IER,PDERV,RESID,NRENEW,YMAX,SAVE1,SAVE2, + SCALE,PW,PWCOPY,ERROR,IPIV,ITOL,RTOL,ATOL,NPSET,NJE,NRE,NDEC + ,IPAR,RPAR,IERR) IF(IERR.NE.0) GOTO 8000 QQQ=QI C -C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY PSET +C NOTE THAT ERROR() IS JUST BEING USED AS A WORKSPACE BY I_PSET IF (IER.NE.0) THEN C IF IER>0 THEN WE HAVE HAD A SINGULARITY IN THE ITERATION MATRIX IJUS=1 @@ -2467,14 +2285,14 @@ C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C ********************************************************************** IF (.NOT.SAMPLE) THEN - CALL ITRAT2(QQQ,Y,YPRIME,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1, + CALL I_ITRAT2(QQQ,Y,YPRIME,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1, + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND, + NIND1,NIND2,NIND3,IPIV,1,ITOL,RTOL,ATOL,IPAR,RPAR, + HUSED,NBSOL,NRE,NQUSED,resid,IERR) IF(IERR.NE.0) GOTO 8000 ELSE - CALL ITRAT2(QQQ,Y,YPRIME,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1, + CALL I_ITRAT2(QQQ,Y,YPRIME,N,T,QI,BND,ARH,CRATE1,TCRAT1,M1, + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND, + NIND1,NIND2,NIND3,IPIV,0,ITOL,RTOL,ATOL,IPAR,RPAR, + HUSED,NBSOL,NRE,NQUSED,resid,IERR) @@ -2589,7 +2407,7 @@ ARH(I) = ARH(I) + EL(JP1)*Y(I,J1) 200 CONTINUE 210 CONTINUE - CALL PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) + CALL I_PRDICT(T,H,Y,L,N,IPAR,RPAR,IERR) IF(IERR.NE.0) GOTO 8000 DO 215 I=1,N YPRIME(I)=(Y(I,1)-ARH(I))/QQQ @@ -2603,7 +2421,7 @@ C FOR NOW WILL ASSUME THAT WE DO NOT WISH TO SAMPLE C AT THE N+2 STEP POINT C - CALL ITRAT2(QQQ,Y,YPRIME,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2, + CALL I_ITRAT2(QQQ,Y,YPRIME,N,T,QI,BND,ARH,CRATE2,TCRAT2,M2, + WORKED,YMAX,ERROR,SAVE1,SAVE2,SCALE,PW,MF,MBND, + NIND1,NIND2,NIND3,IPIV,1,ITOL,RTOL,ATOL,IPAR,RPAR, + HUSED,NBSOL,NRE,NQUSED,resid,IERR) @@ -2661,10 +2479,10 @@ NRE=NRE+1 C IF (MF.GE. 23) THEN - CALL DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) + CALL I_DGBSL(PW,MBND(4),N,MBND(1),MBND(2),IPIV,SAVE1,0) NBSOL=NBSOL+1 ELSE - CALL SOL(N,N,PW,SAVE1,IPIV) + CALL I_SOL(N,N,PW,SAVE1,IPIV) NBSOL = NBSOL + 1 ENDIF DO 321 I=1,N @@ -2758,7 +2576,7 @@ 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) + CALL I_CPYARY(N*L,YHOLD,Y) RMAX = 2.0D+0 IF (DABS(H).LE.HMIN*1.00001D+0) THEN C @@ -2787,10 +2605,10 @@ NQ=NEWQ RH=ONE/(PLFAIL*DBLE(FLOAT(-KFAIL))) L=NQ+1 - CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + CALL I_COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) RC=RC*EL(1)/OLDLO OLDLO=EL(1) - CALL ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) + CALL I_ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) ELSE NEWQ = NQ RH = ONE/ (PRFAIL*DBLE(FLOAT(-KFAIL))) @@ -2816,7 +2634,7 @@ C ********************************* JCHANG = 1 RH = DMAX1(HMIN/DABS(H),0.1D+0) - CALL HCHOSE(RH,H,OVRIDE) + CALL I_HCHOSE(RH,H,OVRIDE) H = H*RH DO 350 I = 1,N Y(I,1) = YHOLD(I,1) @@ -2832,11 +2650,11 @@ NQ = 1 L = 2 C RESET ORDER, RECALCULATE ERROR BOUNDS - CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + CALL I_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) + CALL I_ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) C NOW JUMP TO NORMAL CONTINUATION POINT GO TO 60 C ********************************************************************** @@ -3003,7 +2821,7 @@ GOTO 440 ENDIF RH = DMIN1(RH,RMAX) - CALL HCHOSE(RH,H,OVRIDE) + CALL I_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 @@ -3031,16 +2849,16 @@ NQ = NEWQ L = NQ + 1 C RESET ORDER,RECALCULATE ERROR BOUNDS - CALL COSET(NQ,EL,ELST,TQ,NCOSET,MAXORD) + CALL I_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) + CALL I_ERRORS(N,TQ,EDN,E,EUP,BND,EDDN) END IF RH = DMAX1(RH,HMIN/DABS(H)) RH = DMIN1(RH,HMAX/DABS(H),RMAX) - CALL RSCALE(N,L,RH,Y) + CALL I_RSCALE(N,L,RH,Y) RMAX = 10.0D+0 JCHANG = 1 H = H*RH @@ -3057,7 +2875,7 @@ 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) + CALL I_CPYARY(N*L,Y,YHOLD) NSTEP = NSTEP + 1 JSINUP = JSINUP + 1 JSNOLD = JSNOLD + 1 @@ -3112,7 +2930,7 @@ IF ((T-TOUT)*H.GE.0.0D+0) THEN C HAVE OVERSHOT TOUT WRITE (LOUT,*) T,TOUT,H - CALL INTERP(N,JSTART,H,T,Y,TOUT,Y0) + CALL I_INTERP(N,JSTART,H,T,Y,TOUT,Y0) HO = H T0 = TOUT IDID = -5 @@ -3123,7 +2941,7 @@ goto 30 endif c - IF(IJUS.EQ.0) CALL HCHOSE(RH,H,OVRIDE) + IF(IJUS.EQ.0) CALL I_HCHOSE(RH,H,OVRIDE) IF(.NOT.FINISH) THEN GO TO 40 ELSE @@ -3132,9 +2950,9 @@ 9000 FORMAT (1X,' CORRECTOR HAS NOT CONVERGED') END -C ------------------- END OF SUBROUTINE STIFF -------------------------- +C ------------------- END OF SUBROUTINE I_STIFF -------------------------- - SUBROUTINE RSCALE(N,L,RH,Y) + SUBROUTINE I_RSCALE(N,L,RH,Y) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C .. SCALAR ARGUMENTS .. @@ -3246,7 +3064,7 @@ END C--------------------------------------------------------------------------- - SUBROUTINE CPYARY(NELEM,SOURCE,TARGET) + SUBROUTINE I_CPYARY(NELEM,SOURCE,TARGET) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C COPIES THE ARRAY SOURCE() INTO THE ARRAY TARGET() @@ -3271,7 +3089,7 @@ END C---------------------------------------------------------------------------- - SUBROUTINE HCHOSE(RH,H,OVRIDE) + SUBROUTINE I_HCHOSE(RH,H,OVRIDE) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON / STPSZE / HSTPSZ(2,14) LOGICAL OVRIDE @@ -3306,953 +3124,3 @@ RETURN END -C -C ************************************************************ -C - DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* October 31, 1992 -* -* .. Scalar Arguments .. - CHARACTER CMACH -* .. -* -* Purpose -* ======= -* -* DLAMCH determines double precision machine parameters. -* -* Arguments -* ========= -* -* CMACH (input) CHARACTER*1 -* Specifies the value to be returned by DLAMCH: -* = 'E' or 'e', DLAMCH := eps -* = 'S' or 's , DLAMCH := sfmin -* = 'B' or 'b', DLAMCH := base -* = 'P' or 'p', DLAMCH := eps*base -* = 'N' or 'n', DLAMCH := t -* = 'R' or 'r', DLAMCH := rnd -* = 'M' or 'm', DLAMCH := emin -* = 'U' or 'u', DLAMCH := rmin -* = 'L' or 'l', DLAMCH := emax -* = 'O' or 'o', DLAMCH := rmax -* -* where -* -* eps = relative machine precision -* sfmin = safe minimum, such that 1/sfmin does not overflow -* base = base of the machine -* prec = eps*base -* t = number of (base) digits in the mantissa -* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise -* emin = minimum exponent before (gradual) underflow -* rmin = underflow threshold - base**(emin-1) -* emax = largest exponent before overflow -* rmax = overflow threshold - (base**emax)*(1-eps) -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) -* .. -* .. Local Scalars .. - LOGICAL FIRST, LRND - INTEGER BETA, IMAX, IMIN, IT - DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, - $ RND, SFMIN, SMALL, T -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME -* .. -* .. External Subroutines .. - EXTERNAL DLAMC2 -* .. -* .. Save statement .. - SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, - $ EMAX, RMAX, PREC -* .. -* .. Data statements .. - DATA FIRST / .TRUE. / -* .. -* .. Executable Statements .. -* - IF( FIRST ) THEN - FIRST = .FALSE. - CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) - BASE = BETA - T = IT - IF( LRND ) THEN - RND = ONE - EPS = ( BASE**( 1-IT ) ) / 2 - ELSE - RND = ZERO - EPS = BASE**( 1-IT ) - END IF - PREC = EPS*BASE - EMIN = IMIN - EMAX = IMAX - SFMIN = RMIN - SMALL = ONE / RMAX - IF( SMALL.GE.SFMIN ) THEN -* -* Use SMALL plus a bit, to avoid the possibility of rounding -* causing overflow when computing 1/sfmin. -* - SFMIN = SMALL*( ONE+EPS ) - END IF - END IF -* - IF( LSAME( CMACH, 'E' ) ) THEN - RMACH = EPS - ELSE IF( LSAME( CMACH, 'S' ) ) THEN - RMACH = SFMIN - ELSE IF( LSAME( CMACH, 'B' ) ) THEN - RMACH = BASE - ELSE IF( LSAME( CMACH, 'P' ) ) THEN - RMACH = PREC - ELSE IF( LSAME( CMACH, 'N' ) ) THEN - RMACH = T - ELSE IF( LSAME( CMACH, 'R' ) ) THEN - RMACH = RND - ELSE IF( LSAME( CMACH, 'M' ) ) THEN - RMACH = EMIN - ELSE IF( LSAME( CMACH, 'U' ) ) THEN - RMACH = RMIN - ELSE IF( LSAME( CMACH, 'L' ) ) THEN - RMACH = EMAX - ELSE IF( LSAME( CMACH, 'O' ) ) THEN - RMACH = RMAX - END IF -* - DLAMCH = RMACH - RETURN -* -* End of DLAMCH -* - END -* -************************************************************************ -* - SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* October 31, 1992 -* -* .. Scalar Arguments .. - LOGICAL IEEE1, RND - INTEGER BETA, T -* .. -* -* Purpose -* ======= -* -* DLAMC1 determines the machine parameters given by BETA, T, RND, and -* IEEE1. -* -* Arguments -* ========= -* -* BETA (output) INTEGER -* The base of the machine. -* -* T (output) INTEGER -* The number of ( BETA ) digits in the mantissa. -* -* RND (output) LOGICAL -* Specifies whether proper rounding ( RND = .TRUE. ) or -* chopping ( RND = .FALSE. ) occurs in addition. This may not -* be a reliable guide to the way in which the machine performs -* its arithmetic. -* -* IEEE1 (output) LOGICAL -* Specifies whether rounding appears to be done in the IEEE -* 'round to nearest' style. -* -* Further Details -* =============== -* -* The routine is based on the routine ENVRON by Malcolm and -* incorporates suggestions by Gentleman and Marovich. See -* -* Malcolm M. A. (1972) Algorithms to reveal properties of -* floating-point arithmetic. Comms. of the ACM, 15, 949-951. -* -* Gentleman W. M. and Marovich S. B. (1974) More on algorithms -* that reveal properties of floating point arithmetic units. -* Comms. of the ACM, 17, 276-277. -* -* ===================================================================== -* -* .. Local Scalars .. - LOGICAL FIRST, LIEEE1, LRND - INTEGER LBETA, LT - DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMC3 - EXTERNAL DLAMC3 -* .. -* .. Save statement .. - SAVE FIRST, LIEEE1, LBETA, LRND, LT -* .. -* .. Data statements .. - DATA FIRST / .TRUE. / -* .. -* .. Executable Statements .. -* - IF( FIRST ) THEN - FIRST = .FALSE. - ONE = 1 -* -* LBETA, LIEEE1, LT and LRND are the local values of BETA, -* IEEE1, T and RND. -* -* Throughout this routine we use the function DLAMC3 to ensure -* that relevant values are stored and not held in registers, or -* are not affected by optimizers. -* -* Compute a = 2.0**m with the smallest positive integer m such -* that -* -* fl( a + 1.0 ) = a. -* - A = 1 - C = 1 -* -*+ WHILE( C.EQ.ONE )LOOP - 10 CONTINUE - IF( C.EQ.ONE ) THEN - A = 2*A - C = DLAMC3( A, ONE ) - C = DLAMC3( C, -A ) - GO TO 10 - END IF -*+ END WHILE -* -* Now compute b = 2.0**m with the smallest positive integer m -* such that -* -* fl( a + b ) .gt. a. -* - B = 1 - C = DLAMC3( A, B ) -* -*+ WHILE( C.EQ.A )LOOP - 20 CONTINUE - IF( C.EQ.A ) THEN - B = 2*B - C = DLAMC3( A, B ) - GO TO 20 - END IF -*+ END WHILE -* -* Now compute the base. a and c are neighbouring floating point -* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so -* their difference is beta. Adding 0.25 to c is to ensure that it -* is truncated to beta and not ( beta - 1 ). -* - QTR = ONE / 4 - SAVEC = C - C = DLAMC3( C, -A ) - LBETA = C + QTR -* -* Now determine whether rounding or chopping occurs, by adding a -* bit less than beta/2 and a bit more than beta/2 to a. -* - B = LBETA - F = DLAMC3( B / 2, -B / 100 ) - C = DLAMC3( F, A ) - IF( C.EQ.A ) THEN - LRND = .TRUE. - ELSE - LRND = .FALSE. - END IF - F = DLAMC3( B / 2, B / 100 ) - C = DLAMC3( F, A ) - IF( ( LRND ) .AND. ( C.EQ.A ) ) - $ LRND = .FALSE. -* -* Try and decide whether rounding is done in the IEEE 'round to -* nearest' style. B/2 is half a unit in the last place of the two -* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit -* zero, and SAVEC is odd. Thus adding B/2 to A should not change -* A, but adding B/2 to SAVEC should change SAVEC. -* - T1 = DLAMC3( B / 2, A ) - T2 = DLAMC3( B / 2, SAVEC ) - LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND -* -* Now find the mantissa, t. It should be the integer part of -* log to the base beta of a, however it is safer to determine t -* by powering. So we find t as the smallest positive integer for -* which -* -* fl( beta**t + 1.0 ) = 1.0. -* - LT = 0 - A = 1 - C = 1 -* -*+ WHILE( C.EQ.ONE )LOOP - 30 CONTINUE - IF( C.EQ.ONE ) THEN - LT = LT + 1 - A = A*LBETA - C = DLAMC3( A, ONE ) - C = DLAMC3( C, -A ) - GO TO 30 - END IF -*+ END WHILE -* - END IF -* - BETA = LBETA - T = LT - RND = LRND - IEEE1 = LIEEE1 - RETURN -* -* End of DLAMC1 -* - END -* -************************************************************************ -* - SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* October 31, 1992 -* -* .. Scalar Arguments .. - LOGICAL RND - INTEGER BETA, EMAX, EMIN, T - DOUBLE PRECISION EPS, RMAX, RMIN -* .. -* -* Purpose -* ======= -* -* DLAMC2 determines the machine parameters specified in its argument -* list. -* -* Arguments -* ========= -* -* BETA (output) INTEGER -* The base of the machine. -* -* T (output) INTEGER -* The number of ( BETA ) digits in the mantissa. -* -* RND (output) LOGICAL -* Specifies whether proper rounding ( RND = .TRUE. ) or -* chopping ( RND = .FALSE. ) occurs in addition. This may not -* be a reliable guide to the way in which the machine performs -* its arithmetic. -* -* EPS (output) DOUBLE PRECISION -* The smallest positive number such that -* -* fl( 1.0 - EPS ) .LT. 1.0, -* -* where fl denotes the computed value. -* -* EMIN (output) INTEGER -* The minimum exponent before (gradual) underflow occurs. -* -* RMIN (output) DOUBLE PRECISION -* The smallest normalized number for the machine, given by -* BASE**( EMIN - 1 ), where BASE is the floating point value -* of BETA. -* -* EMAX (output) INTEGER -* The maximum exponent before overflow occurs. -* -* RMAX (output) DOUBLE PRECISION -* The largest positive number for the machine, given by -* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point -* value of BETA. -* -* Further Details -* =============== -* -* The computation of EPS is based on a routine PARANOIA by -* W. Kahan of the University of California at Berkeley. -* -* ===================================================================== -* -* .. Local Scalars .. - LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND - INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, - $ NGNMIN, NGPMIN - DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, - $ SIXTH, SMALL, THIRD, TWO, ZERO -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMC3 - EXTERNAL DLAMC3 -* .. -* .. External Subroutines .. - EXTERNAL DLAMC1, DLAMC4, DLAMC5 -* .. -* .. Intrinsic Functions .. - INTRINSIC ABS, MAX, MIN -* .. -* .. Save statement .. - SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, - $ LRMIN, LT -* .. -* .. Data statements .. - DATA FIRST / .TRUE. / , IWARN / .FALSE. / -* .. -* .. Executable Statements .. -* - IF( FIRST ) THEN - FIRST = .FALSE. - ZERO = 0 - ONE = 1 - TWO = 2 -* -* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of -* BETA, T, RND, EPS, EMIN and RMIN. -* -* Throughout this routine we use the function DLAMC3 to ensure -* that relevant values are stored and not held in registers, or -* are not affected by optimizers. -* -* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. -* - CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) -* -* Start to find EPS. -* - B = LBETA - A = B**( -LT ) - LEPS = A -* -* Try some tricks to see whether or not this is the correct EPS. -* - B = TWO / 3 - HALF = ONE / 2 - SIXTH = DLAMC3( B, -HALF ) - THIRD = DLAMC3( SIXTH, SIXTH ) - B = DLAMC3( THIRD, -HALF ) - B = DLAMC3( B, SIXTH ) - B = ABS( B ) - IF( B.LT.LEPS ) - $ B = LEPS -* - LEPS = 1 -* -*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP - 10 CONTINUE - IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN - LEPS = B - C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) - C = DLAMC3( HALF, -C ) - B = DLAMC3( HALF, C ) - C = DLAMC3( HALF, -B ) - B = DLAMC3( HALF, C ) - GO TO 10 - END IF -*+ END WHILE -* - IF( A.LT.LEPS ) - $ LEPS = A -* -* Computation of EPS complete. -* -* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). -* Keep dividing A by BETA until (gradual) underflow occurs. This -* is detected when we cannot recover the previous A. -* - RBASE = ONE / LBETA - SMALL = ONE - DO 20 I = 1, 3 - SMALL = DLAMC3( SMALL*RBASE, ZERO ) - 20 CONTINUE - A = DLAMC3( ONE, SMALL ) - CALL DLAMC4( NGPMIN, ONE, LBETA ) - CALL DLAMC4( NGNMIN, -ONE, LBETA ) - CALL DLAMC4( GPMIN, A, LBETA ) - CALL DLAMC4( GNMIN, -A, LBETA ) - IEEE = .FALSE. -* - IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN - IF( NGPMIN.EQ.GPMIN ) THEN - LEMIN = NGPMIN -* ( Non twos-complement machines, no gradual underflow; -* e.g., VAX ) - ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN - LEMIN = NGPMIN - 1 + LT - IEEE = .TRUE. -* ( Non twos-complement machines, with gradual underflow; -* e.g., IEEE standard followers ) - ELSE - LEMIN = MIN( NGPMIN, GPMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -* - ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN - IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN - LEMIN = MAX( NGPMIN, NGNMIN ) -* ( Twos-complement machines, no gradual underflow; -* e.g., CYBER 205 ) - ELSE - LEMIN = MIN( NGPMIN, NGNMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -* - ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. - $ ( GPMIN.EQ.GNMIN ) ) THEN - IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN - LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT -* ( Twos-complement machines with gradual underflow; -* no known machine ) - ELSE - LEMIN = MIN( NGPMIN, NGNMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -* - ELSE - LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) -* ( A guess; no known machine ) - IWARN = .TRUE. - END IF -*** -* Comment out this if block if EMIN is ok - IF( IWARN ) THEN - FIRST = .TRUE. - WRITE( 6, FMT = 9999 )LEMIN - END IF -*** -* -* Assume IEEE arithmetic if we found denormalised numbers above, -* or if arithmetic seems to round in the IEEE style, determined -* in routine DLAMC1. A true IEEE machine should have both things -* true; however, faulty machines may have one or the other. -* - IEEE = IEEE .OR. LIEEE1 -* -* Compute RMIN by successive division by BETA. We could compute -* RMIN as BASE**( EMIN - 1 ), but some machines underflow during -* this computation. -* - LRMIN = 1 - DO 30 I = 1, 1 - LEMIN - LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) - 30 CONTINUE -* -* Finally, call DLAMC5 to compute EMAX and RMAX. -* - CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) - END IF -* - BETA = LBETA - T = LT - RND = LRND - EPS = LEPS - EMIN = LEMIN - RMIN = LRMIN - EMAX = LEMAX - RMAX = LRMAX -* - RETURN -* - 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', - $ ' EMIN = ', I8, / - $ ' If, after inspection, the value EMIN looks', - $ ' acceptable please comment out ', - $ / ' the IF block as marked within the code of routine', - $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) -* -* End of DLAMC2 -* - END -* -************************************************************************ -* - DOUBLE PRECISION FUNCTION DLAMC3( A, B ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* October 31, 1992 -* -* .. Scalar Arguments .. - DOUBLE PRECISION A, B -* .. -* -* Purpose -* ======= -* -* DLAMC3 is intended to force A and B to be stored prior to doing -* the addition of A and B , for use in situations where optimizers -* might hold one of these in a register. -* -* Arguments -* ========= -* -* A, B (input) DOUBLE PRECISION -* The values A and B. -* -* ===================================================================== -* -* .. Executable Statements .. -* - DLAMC3 = A + B -* - RETURN -* -* End of DLAMC3 -* - END -* -************************************************************************ -* - SUBROUTINE DLAMC4( EMIN, START, BASE ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* October 31, 1992 -* -* .. Scalar Arguments .. - INTEGER BASE, EMIN - DOUBLE PRECISION START -* .. -* -* Purpose -* ======= -* -* DLAMC4 is a service routine for DLAMC2. -* -* Arguments -* ========= -* -* EMIN (output) EMIN -* The minimum exponent before (gradual) underflow, computed by -* setting A = START and dividing by BASE until the previous A -* can not be recovered. -* -* START (input) DOUBLE PRECISION -* The starting point for determining EMIN. -* -* BASE (input) INTEGER -* The base of the machine. -* -* ===================================================================== -* -* .. Local Scalars .. - INTEGER I - DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMC3 - EXTERNAL DLAMC3 -* .. -* .. Executable Statements .. -* - A = START - ONE = 1 - RBASE = ONE / BASE - ZERO = 0 - EMIN = 1 - B1 = DLAMC3( A*RBASE, ZERO ) - C1 = A - C2 = A - D1 = A - D2 = A -*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. -* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP - 10 CONTINUE - IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. - $ ( D2.EQ.A ) ) THEN - EMIN = EMIN - 1 - A = B1 - B1 = DLAMC3( A / BASE, ZERO ) - C1 = DLAMC3( B1*BASE, ZERO ) - D1 = ZERO - DO 20 I = 1, BASE - D1 = D1 + B1 - 20 CONTINUE - B2 = DLAMC3( A*RBASE, ZERO ) - C2 = DLAMC3( B2 / RBASE, ZERO ) - D2 = ZERO - DO 30 I = 1, BASE - D2 = D2 + B2 - 30 CONTINUE - GO TO 10 - END IF -*+ END WHILE -* - RETURN -* -* End of DLAMC4 -* - END -* -************************************************************************ -* - SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* October 31, 1992 -* -* .. Scalar Arguments .. - LOGICAL IEEE - INTEGER BETA, EMAX, EMIN, P - DOUBLE PRECISION RMAX -* .. -* -* Purpose -* ======= -* -* DLAMC5 attempts to compute RMAX, the largest machine floating-point -* number, without overflow. It assumes that EMAX + abs(EMIN) sum -* approximately to a power of 2. It will fail on machines where this -* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, -* EMAX = 28718). It will also fail if the value supplied for EMIN is -* too large (i.e. too close to zero), probably with overflow. -* -* Arguments -* ========= -* -* BETA (input) INTEGER -* The base of floating-point arithmetic. -* -* P (input) INTEGER -* The number of base BETA digits in the mantissa of a -* floating-point value. -* -* EMIN (input) INTEGER -* The minimum exponent before (gradual) underflow. -* -* IEEE (input) LOGICAL -* A logical flag specifying whether or not the arithmetic -* system is thought to comply with the IEEE standard. -* -* EMAX (output) INTEGER -* The largest exponent before overflow -* -* RMAX (output) DOUBLE PRECISION -* The largest machine floating-point number. -* -* ===================================================================== -* -* .. Parameters .. - DOUBLE PRECISION ZERO, ONE - PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) -* .. -* .. Local Scalars .. - INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP - DOUBLE PRECISION OLDY, RECBAS, Y, Z -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMC3 - EXTERNAL DLAMC3 -* .. -* .. Intrinsic Functions .. - INTRINSIC MOD -* .. -* .. Executable Statements .. -* -* First compute LEXP and UEXP, two powers of 2 that bound -* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum -* approximately to the bound that is closest to abs(EMIN). -* (EMAX is the exponent of the required number RMAX). -* - LEXP = 1 - EXBITS = 1 - 10 CONTINUE - TRY = LEXP*2 - IF( TRY.LE.( -EMIN ) ) THEN - LEXP = TRY - EXBITS = EXBITS + 1 - GO TO 10 - END IF - IF( LEXP.EQ.-EMIN ) THEN - UEXP = LEXP - ELSE - UEXP = TRY - EXBITS = EXBITS + 1 - END IF -* -* Now -LEXP is less than or equal to EMIN, and -UEXP is greater -* than or equal to EMIN. EXBITS is the number of bits needed to -* store the exponent. -* - IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN - EXPSUM = 2*LEXP - ELSE - EXPSUM = 2*UEXP - END IF -* -* EXPSUM is the exponent range, approximately equal to -* EMAX - EMIN + 1 . -* - EMAX = EXPSUM + EMIN - 1 - NBITS = 1 + EXBITS + P -* -* NBITS is the total number of bits needed to store a -* floating-point number. -* - IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN -* -* Either there are an odd number of bits used to store a -* floating-point number, which is unlikely, or some bits are -* not used in the representation of numbers, which is possible, -* (e.g. Cray machines) or the mantissa has an implicit bit, -* (e.g. IEEE machines, Dec Vax machines), which is perhaps the -* most likely. We have to assume the last alternative. -* If this is true, then we need to reduce EMAX by one because -* there must be some way of representing zero in an implicit-bit -* system. On machines like Cray, we are reducing EMAX by one -* unnecessarily. -* - EMAX = EMAX - 1 - END IF -* - IF( IEEE ) THEN -* -* Assume we are on an IEEE machine which reserves one exponent -* for infinity and NaN. -* - EMAX = EMAX - 1 - END IF -* -* Now create RMAX, the largest machine number, which should -* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . -* -* First compute 1.0 - BETA**(-P), being careful that the -* result is less than 1.0 . -* - RECBAS = ONE / BETA - Z = BETA - ONE - Y = ZERO - DO 20 I = 1, P - Z = Z*RECBAS - IF( Y.LT.ONE ) - $ OLDY = Y - Y = DLAMC3( Y, Z ) - 20 CONTINUE - IF( Y.GE.ONE ) - $ Y = OLDY -* -* Now multiply by BETA**EMAX to get RMAX. -* - DO 30 I = 1, EMAX - Y = DLAMC3( Y*BETA, ZERO ) - 30 CONTINUE -* - RMAX = Y - RETURN -* -* End of DLAMC5 -* - END - LOGICAL FUNCTION LSAME( CA, CB ) -* -* -- LAPACK auxiliary routine (version 2.0) -- -* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., -* Courant Institute, Argonne National Lab, and Rice University -* September 30, 1994 -* -* .. Scalar Arguments .. - CHARACTER CA, CB -* .. -* -* Purpose -* ======= -* -* LSAME returns .TRUE. if CA is the same letter as CB regardless of -* case. -* -* Arguments -* ========= -* -* CA (input) CHARACTER*1 -* CB (input) CHARACTER*1 -* CA and CB specify the single characters to be compared. -* -* ===================================================================== -* -* .. Intrinsic Functions .. - INTRINSIC ICHAR -* .. -* .. Local Scalars .. - INTEGER INTA, INTB, ZCODE -* .. -* .. Executable Statements .. -* -* Test if the characters are equal -* - LSAME = CA.EQ.CB - IF( LSAME ) - $ RETURN -* -* Now test for equivalence if both characters are alphabetic. -* - ZCODE = ICHAR( 'Z' ) -* -* Use 'Z' rather than 'A' so that ASCII can be detected on Prime -* machines, on which ICHAR returns a value with bit 8 set. -* ICHAR('A') on Prime machines returns 193 which is the same as -* ICHAR('A') on an EBCDIC machine. -* - INTA = ICHAR( CA ) - INTB = ICHAR( CB ) -* - IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN -* -* ASCII is assumed - ZCODE is the ASCII code of either lower or -* upper case 'Z'. -* - IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 - IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 -* - ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN -* -* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or -* upper case 'Z'. -* - IF( INTA.GE.129 .AND. INTA.LE.137 .OR. - $ INTA.GE.145 .AND. INTA.LE.153 .OR. - $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 - IF( INTB.GE.129 .AND. INTB.LE.137 .OR. - $ INTB.GE.145 .AND. INTB.LE.153 .OR. - $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 -* - ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN -* -* ASCII is assumed, on Prime machines - ZCODE is the ASCII code -* plus 128 of either lower or upper case 'Z'. -* - IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 - IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 - END IF - LSAME = INTA.EQ.INTB -* -* RETURN -* -* End of LSAME -* - END - -C---------------------------------------------------------------------------- -