C$Id: cacm423.f,v 1.4 1997/11/26 11:04:17 walter Rel $
CFPP$ NOCONCUR F
C
C CAUTION: this is *only* for use by PSIDE
C
C For PSIDE's linear algebra we chose to use LAPACK.
C However, if you do *not* have available on your system both
C
C - a machine tuned LAPACK, and
C - a machine optimized BLAS
C
C we suggest you use these replacement routines instead of the portable
C LAPACK implementation (available from http://www.netlib.org/lapack/).
C Of course the latter will work, however, it may give much worse
C performance.
C
C DGBTRF/DGBTRS and DGETRF/DGETRS are not the original LAPACK routines
C they are wrappers to mimic LAPACK using CACM423 DEC/SOL and DECB/SOLB
C this works well for PSIDE because we only factorize
C
C - square systems M = N
C - non transposed matrices TRANS = 'N'
C - single right hand side NRHS = 1
C
C the supplied DLAMCH and LSAME are the original LAPACK routines
C
SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO )
INTEGER INFO, KL, KU, LDAB, M, N
INTEGER IPIV( * )
DOUBLE PRECISION AB( LDAB, * )
CALL DECB(N,LDAB,AB,KL,KU,IPIV,INFO)
RETURN
END
SUBROUTINE DGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB,
+ INFO )
CHARACTER TRANS
INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS
INTEGER IPIV( * )
DOUBLE PRECISION AB( LDAB, * ), B( LDB, * )
CALL SOLB(N,LDAB,AB,KL,KU,B,IPIV)
RETURN
END
SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
INTEGER INFO, LDA, M, N
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * )
CALL DEC(N,LDA,A,IPIV,INFO)
RETURN
END
SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
CHARACTER TRANS
INTEGER INFO, LDA, LDB, N, NRHS
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
CALL SOL(N,LDA,A,B,IPIV)
RETURN
END
SUBROUTINE DEC (N, NDIM, A, IP, IER)
INTEGER N, NDIM, IP(N), IER
DOUBLE PRECISION A(NDIM,N)
C-----------------------------------------------------------------------
C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C INPUT..
C N = ORDER OF MATRIX.
C NDIM = DECLARED DIMENSION OF ARRAY A .
C A = MATRIX TO BE TRIANGULARIZED.
C OUTPUT..
C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
C SINGULAR AT STAGE K.
C USE 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)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
C
C REFERENCE..
C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
C C.A.C.M. 15 (1972), P. 274.
C-----------------------------------------------------------------------
INTEGER I, J, K, KP1, M, NM1
DOUBLE PRECISION T
IER = 0
IP(N) = 1
IF (N .EQ. 1) GO TO 70
NM1 = N - 1
DO 60 K = 1,NM1
KP1 = K + 1
M = K
DO 10 I = KP1,N
IF (ABS(A(I,K)) .GT. ABS(A(M,K))) M = I
10 CONTINUE
IP(K) = M
T = A(M,K)
IF (M .EQ. K) GO TO 20
IP(N) = -IP(N)
A(M,K) = A(K,K)
A(K,K) = T
20 CONTINUE
IF (T .EQ. 0.D0) GO TO 80
T = 1.D0/T
DO 30 I = KP1,N
30 A(I,K) = -A(I,K)*T
DO 50 J = KP1,N
T = A(M,J)
A(M,J) = A(K,J)
A(K,J) = T
IF (T .EQ. 0.D0) GO TO 45
DO 40 I = KP1,N
40 A(I,J) = A(I,J) + A(I,K)*T
45 CONTINUE
50 CONTINUE
60 CONTINUE
70 K = N
IF (A(N,N) .EQ. 0.D0) GO TO 80
RETURN
80 IER = K
IP(N) = 0
RETURN
C----------------------- END OF SUBROUTINE DEC -------------------------
END
SUBROUTINE SOL (N, NDIM, A, B, IP)
INTEGER N, NDIM, IP(N)
DOUBLE PRECISION A(NDIM,N), B(N)
C-----------------------------------------------------------------------
C SOLUTION OF LINEAR SYSTEM, A*X = B .
C INPUT..
C N = ORDER OF MATRIX.
C NDIM = DECLARED DIMENSION OF ARRAY A .
C A = TRIANGULARIZED 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 OUTPUT..
C B = SOLUTION VECTOR, X .
C-----------------------------------------------------------------------
INTEGER I, K, KB, KM1, KP1, M, NM1
DOUBLE PRECISION T
IF (N .EQ. 1) GO TO 50
NM1 = N - 1
DO 20 K = 1,NM1
KP1 = K + 1
M = IP(K)
T = B(M)
B(M) = B(K)
B(K) = T
DO 10 I = KP1,N
10 B(I) = B(I) + A(I,K)*T
20 CONTINUE
DO 40 KB = 1,NM1
KM1 = N - KB
K = KM1 + 1
B(K) = B(K)/A(K,K)
T = -B(K)
DO 30 I = 1,KM1
30 B(I) = B(I) + A(I,K)*T
40 CONTINUE
50 B(1) = B(1)/A(1,1)
RETURN
C----------------------- END OF SUBROUTINE SOL -------------------------
END
SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER)
INTEGER N, NDIM, ML, MU, IP(N), IER
DOUBLE PRECISION A(NDIM,N)
C-----------------------------------------------------------------------
C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
C MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
C INPUT..
C N ORDER OF THE ORIGINAL MATRIX A.
C NDIM DECLARED DIMENSION OF ARRAY A.
C A CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS
C OF THE MATRIX ARE STORED IN THE COLUMNS OF A AND
C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
C ML+1 THROUGH 2*ML+MU+1 OF A.
C ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C OUTPUT..
C A AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
C IP INDEX VECTOR OF PIVOT INDICES.
C IP(N) (-1)**(NUMBER OF INTERCHANGES) OR O .
C IER = 0 IF MATRIX A IS NONSINGULAR, OR = K IF FOUND TO BE
C SINGULAR AT STAGE K.
C USE SOLB TO OBTAIN SOLUTION OF LINEAR SYSTEM.
C DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N) WITH MD=ML+MU+1.
C IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
C
C REFERENCE..
C THIS IS A MODIFICATION OF
C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
C C.A.C.M. 15 (1972), P. 274.
C-----------------------------------------------------------------------
INTEGER I, IJK, J, JK, JU, K, KP1, M, MD, MD1, MDL, MM, NM1
DOUBLE PRECISION T
IER = 0
IP(N) = 1
MD = ML + MU + 1
MD1 = MD + 1
JU = 0
IF (ML .EQ. 0) GO TO 70
IF (N .EQ. 1) GO TO 70
IF (N .LT. MU+2) GO TO 7
DO 5 J = MU+2,N
DO 5 I = 1,ML
5 A(I,J) = 0.D0
7 NM1 = N - 1
DO 60 K = 1,NM1
KP1 = K + 1
M = MD
MDL = MIN(ML,N-K) + MD
DO 10 I = MD1,MDL
IF (ABS(A(I,K)) .GT. ABS(A(M,K))) M = I
10 CONTINUE
IP(K) = M + K - MD
T = A(M,K)
IF (M .EQ. MD) GO TO 20
IP(N) = -IP(N)
A(M,K) = A(MD,K)
A(MD,K) = T
20 CONTINUE
IF (T .EQ. 0.D0) GO TO 80
T = 1.D0/T
DO 30 I = MD1,MDL
30 A(I,K) = -A(I,K)*T
JU = MIN(MAX(JU,MU+IP(K)),N)
MM = MD
IF (JU .LT. KP1) GO TO 55
DO 50 J = KP1,JU
M = M - 1
MM = MM - 1
T = A(M,J)
IF (M .EQ. MM) GO TO 35
A(M,J) = A(MM,J)
A(MM,J) = T
35 CONTINUE
IF (T .EQ. 0.D0) GO TO 45
JK = J - K
DO 40 I = MD1,MDL
IJK = I - JK
40 A(IJK,J) = A(IJK,J) + A(I,K)*T
45 CONTINUE
50 CONTINUE
55 CONTINUE
60 CONTINUE
70 K = N
IF (A(MD,N) .EQ. 0.D0) GO TO 80
RETURN
80 IER = K
IP(N) = 0
RETURN
C----------------------- END OF SUBROUTINE DECB ------------------------
END
SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP)
INTEGER N, NDIM, ML, MU, IP(N)
DOUBLE PRECISION A(NDIM,N), B(N)
C-----------------------------------------------------------------------
C SOLUTION OF LINEAR SYSTEM, A*X = B .
C INPUT..
C N ORDER OF MATRIX A.
C NDIM DECLARED DIMENSION OF ARRAY A .
C A TRIANGULARIZED MATRIX OBTAINED FROM DECB.
C ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
C B RIGHT HAND SIDE VECTOR.
C IP PIVOT VECTOR OBTAINED FROM DECB.
C DO NOT USE IF DECB HAS SET IER .NE. 0.
C OUTPUT..
C B SOLUTION VECTOR, X .
C-----------------------------------------------------------------------
INTEGER I, IMD, K, KB, KMD, LM, M, MD, MD1, MDL, MDM, NM1
DOUBLE PRECISION T
MD = ML + MU + 1
MD1 = MD + 1
MDM = MD - 1
NM1 = N - 1
IF (ML .EQ. 0) GO TO 25
IF (N .EQ. 1) GO TO 50
DO 20 K = 1,NM1
M = IP(K)
T = B(M)
B(M) = B(K)
B(K) = T
MDL = MIN(ML,N-K) + MD
DO 10 I = MD1,MDL
IMD = I + K - MD
10 B(IMD) = B(IMD) + A(I,K)*T
20 CONTINUE
25 CONTINUE
DO 40 KB = 1,NM1
K = N + 1 - KB
B(K) = B(K)/A(MD,K)
T = -B(K)
KMD = MD - K
LM = MAX(1,KMD+1)
DO 30 I = LM,MDM
IMD = I - KMD
30 B(IMD) = B(IMD) + A(I,K)*T
40 CONTINUE
50 B(1) = B(1)/A(MD,1)
RETURN
C----------------------- END OF SUBROUTINE SOLB ------------------------
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