网格计算

开发平台：

Windows_Unix

lapack_agmg_win.f90：源码内容

SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
use MSIMSL
!
! -- LAPACK routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
! ..
! .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * )
! ..
!
! Purpose
! =======
!
! DGETF2 computes an LU factorization of a general m-by-n matrix A
! using partial pivoting with row interchanges.
!
! The factorization has the form
! A = P * L * U
! where P is a permutation matrix, L is lower triangular with unit
! diagonal elements (lower trapezoidal if m > n), and U is upper
! triangular (upper trapezoidal if m < n).
!
! This is the right-looking Level 2 BLAS version of the algorithm.
!
! Arguments
! =========
!
! M (input) INTEGER
! The number of rows of the matrix A. M >= 0.
!
! N (input) INTEGER
! The number of columns of the matrix A. N >= 0.
!
! A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! On entry, the m by n matrix to be factored.
! On exit, the factors L and U from the factorization
! A = P*L*U; the unit diagonal elements of L are not stored.
!
! LDA (input) INTEGER
! The leading dimension of the array A. LDA >= max(1,M).
!
! IPIV (output) INTEGER array, dimension (min(M,N))
! The pivot indices; for 1 <= i <= min(M,N), row i of the
! matrix was interchanged with row IPIV(i).
!
! INFO (output) INTEGER
! = 0: successful exit
! < 0: if INFO = -k, the k-th argument had an illegal value
! > 0: if INFO = k, U(k,k) is exactly zero. The factorization
! has been completed, but the factor U is exactly
! singular, and division by zero will occur if it is used
! to solve a system of equations.
!
! =====================================================================
!
! .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! ..
! .. Local Scalars ..
DOUBLE PRECISION SFMIN
INTEGER I, J, JP
! ..
! .. External Functions ..
DOUBLE PRECISION DLAMCH
!msimslINTEGER IDAMAX
EXTERNAL DLAMCH!msimsl, IDAMAX
! ..
! .. External Subroutines ..
!msimslEXTERNAL DGER, DSCAL, DSWAP, XERBLA
EXTERNAL XERBLA
! ..
! .. Intrinsic Functions ..
INTRINSIC MAX, MIN
! ..
! .. Executable Statements ..
!
! Test the input parameters.
!
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGETF2', -INFO )
RETURN
END IF
!
! Quick return if possible
!
IF( M.EQ.0 .OR. N.EQ.0 ) &
& RETURN
!
! Compute machine safe minimum
!
SFMIN = DLAMCH('S')
!
DO 10 J = 1, MIN( M, N )
!
! Find pivot and test for singularity.
!
JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 )
IPIV( J ) = JP
IF( A( JP, J ).NE.ZERO ) THEN
!
! Apply the interchange to columns 1:N.
!
IF( JP.NE.J ) &
& CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
!
! Compute elements J+1:M of J-th column.
!
IF( J.LT.M ) THEN
IF( ABS(A( J, J )) .GE. SFMIN ) THEN
CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
ELSE
DO 20 I = 1, M-J
A( J+I, J ) = A( J+I, J ) / A( J, J )
20 CONTINUE
END IF
END IF
!
ELSE IF( INFO.EQ.0 ) THEN
!
INFO = J
END IF
!
IF( J.LT.MIN( M, N ) ) THEN
!
! Update trailing submatrix.
!
CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA,&
& A( J+1, J+1 ), LDA )
END IF
10 END DO
RETURN
!
! End of DGETF2
!
END
SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
use MSIMSL
!
! -- LAPACK routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
! ..
! .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * )
! ..
!
! Purpose
! =======
!
! DGETRF computes an LU factorization of a general M-by-N matrix A
! using partial pivoting with row interchanges.
!
! The factorization has the form
! A = P * L * U
! where P is a permutation matrix, L is lower triangular with unit
! diagonal elements (lower trapezoidal if m > n), and U is upper
! triangular (upper trapezoidal if m < n).
!
! This is the right-looking Level 3 BLAS version of the algorithm.
!
! Arguments
! =========
!
! M (input) INTEGER
! The number of rows of the matrix A. M >= 0.
!
! N (input) INTEGER
! The number of columns of the matrix A. N >= 0.
!
! A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! On entry, the M-by-N matrix to be factored.
! On exit, the factors L and U from the factorization
! A = P*L*U; the unit diagonal elements of L are not stored.
!
! LDA (input) INTEGER
! The leading dimension of the array A. LDA >= max(1,M).
!
! IPIV (output) INTEGER array, dimension (min(M,N))
! The pivot indices; for 1 <= i <= min(M,N), row i of the
! matrix was interchanged with row IPIV(i).
!
! INFO (output) INTEGER
! = 0: successful exit
! < 0: if INFO = -i, the i-th argument had an illegal value
! > 0: if INFO = i, U(i,i) is exactly zero. The factorization
! has been completed, but the factor U is exactly
! singular, and division by zero will occur if it is used
! to solve a system of equations.
!
! =====================================================================
!
! .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
! ..
! .. Local Scalars ..
INTEGER I, IINFO, J, JB, NB
! ..
! .. External Subroutines ..
!msimslEXTERNAL DGEMM, DTRSM,
EXTERNAL DGETF2, DLASWP, XERBLA
! ..
! .. External Functions ..
INTEGER ILAENV
EXTERNAL ILAENV
! ..
! .. Intrinsic Functions ..
INTRINSIC MAX, MIN
! ..
! .. Executable Statements ..
!
! Test the input parameters.
!
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGETRF', -INFO )
RETURN
END IF
!
! Quick return if possible
!
IF( M.EQ.0 .OR. N.EQ.0 ) &
& RETURN
!
! Determine the block size for this environment.
!
NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 )
IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
!
! Use unblocked code.
!
CALL DGETF2( M, N, A, LDA, IPIV, INFO )
ELSE
!
! Use blocked code.
!
DO 20 J = 1, MIN( M, N ), NB
JB = MIN( MIN( M, N )-J+1, NB )
!
! Factor diagonal and subdiagonal blocks and test for exact
! singularity.
!
CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
!
! Adjust INFO and the pivot indices.
!
IF( INFO.EQ.0 .AND. IINFO.GT.0 ) &
& INFO = IINFO + J - 1
DO 10 I = J, MIN( M, J+JB-1 )
IPIV( I ) = J - 1 + IPIV( I )
10 CONTINUE
!
! Apply interchanges to columns 1:J-1.
!
CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
!
IF( J+JB.LE.N ) THEN
!
! Apply interchanges to columns J+JB:N.
!
CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, &
& IPIV, 1 )
!
! Compute block row of U.
!
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB, &
& N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), &
& LDA )
IF( J+JB.LE.M ) THEN
!
! Update trailing submatrix.
!
CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1, &
& N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, &
& A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), &
& LDA )
END IF
END IF
20 CONTINUE
END IF
RETURN
!
! End of DGETRF
!
END
SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
use MSIMSL
!
! -- LAPACK routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. Scalar Arguments ..
CHARACTER TRANS
INTEGER INFO, LDA, LDB, N, NRHS
! ..
! .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
! ..
!
! Purpose
! =======
!
! DGETRS solves a system of linear equations
! A * X = B or A' * X = B
! with a general N-by-N matrix A using the LU factorization computed
! by DGETRF.
!
! Arguments
! =========
!
! TRANS (input) CHARACTER*1
! Specifies the form of the system of equations:
! = 'N': A * X = B (No transpose)
! = 'T': A'* X = B (Transpose)
! = 'C': A'* X = B (Conjugate transpose = Transpose)
!
! N (input) INTEGER
! The order of the matrix A. N >= 0.
!
! NRHS (input) INTEGER
! The number of right hand sides, i.e., the number of columns
! of the matrix B. NRHS >= 0.
!
! A (input) DOUBLE PRECISION array, dimension (LDA,N)
! The factors L and U from the factorization A = P*L*U
! as computed by DGETRF.
!
! LDA (input) INTEGER
! The leading dimension of the array A. LDA >= max(1,N).
!
! IPIV (input) INTEGER array, dimension (N)
! The pivot indices from DGETRF; for 1<=i<=N, row i of the
! matrix was interchanged with row IPIV(i).
!
! B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
! On entry, the right hand side matrix B.
! On exit, the solution matrix X.
!
! LDB (input) INTEGER
! The leading dimension of the array B. LDB >= max(1,N).
!
! INFO (output) INTEGER
! = 0: successful exit
! < 0: if INFO = -i, the i-th argument had an illegal value
!
! =====================================================================
!
! .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
! ..
! .. Local Scalars ..
LOGICAL NOTRAN
! ..
! .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
! ..
! .. External Subroutines ..
!msimslEXTERNAL , DTRSM,
EXTERNAL DLASWP, XERBLA
! ..
! .. Intrinsic Functions ..
INTRINSIC MAX
! ..
! .. Executable Statements ..
!
! Test the input parameters.
!
INFO = 0
NOTRAN = LSAME( TRANS, 'N' )
IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. &
& LSAME( TRANS, 'C' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( NRHS.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -8
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGETRS', -INFO )
RETURN
END IF
!
! Quick return if possible
!
IF( N.EQ.0 .OR. NRHS.EQ.0 ) &
& RETURN
!
IF( NOTRAN ) THEN
!
! Solve A * X = B.
!
! Apply row interchanges to the right hand sides.
!
CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, 1 )
!
! Solve L*X = B, overwriting B with X.
!
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS, &
& ONE, A, LDA, B, LDB )
!
! Solve U*X = B, overwriting B with X.
!
CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N, &
& NRHS, ONE, A, LDA, B, LDB )
ELSE
!
! Solve A' * X = B.
!
! Solve U'*X = B, overwriting B with X.
!
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS, &
& ONE, A, LDA, B, LDB )
!
! Solve L'*X = B, overwriting B with X.
!
CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE,&
& A, LDA, B, LDB )
!
! Apply row interchanges to the solution vectors.
!
CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, -1 )
END IF
!
RETURN
!
! End of DGETRS
!
END
SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX )
!
! -- LAPACK auxiliary routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. Scalar Arguments ..
INTEGER INCX, K1, K2, LDA, N
! ..
! .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * )
! ..
!
! Purpose
! =======
!
! DLASWP performs a series of row interchanges on the matrix A.
! One row interchange is initiated for each of rows K1 through K2 of A.
!
! Arguments
! =========
!
! N (input) INTEGER
! The number of columns of the matrix A.
!
! A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! On entry, the matrix of column dimension N to which the row
! interchanges will be applied.
! On exit, the permuted matrix.
!
! LDA (input) INTEGER
! The leading dimension of the array A.
!
! K1 (input) INTEGER
! The first element of IPIV for which a row interchange will
! be done.
!
! K2 (input) INTEGER
! The last element of IPIV for which a row interchange will
! be done.
!
! IPIV (input) INTEGER array, dimension (K2*abs(INCX))
! The vector of pivot indices. Only the elements in positions
! K1 through K2 of IPIV are accessed.
! IPIV(K) = L implies rows K and L are to be interchanged.
!
! INCX (input) INTEGER
! The increment between successive values of IPIV. If IPIV
! is negative, the pivots are applied in reverse order.
!
! Further Details
! ===============
!
! Modified by
! R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA
!
! =====================================================================
!
! .. Local Scalars ..
INTEGER I, I1, I2, INC, IP, IX, IX0, J, K, N32
DOUBLE PRECISION TEMP
! ..
! .. Executable Statements ..
!
! Interchange row I with row IPIV(I) for each of rows K1 through K2.
!
IF( INCX.GT.0 ) THEN
IX0 = K1
I1 = K1
I2 = K2
INC = 1
ELSE IF( INCX.LT.0 ) THEN
IX0 = 1 + ( 1-K2 )*INCX
I1 = K2
I2 = K1
INC = -1
ELSE
RETURN
END IF
!
N32 = ( N / 32 )*32
IF( N32.NE.0 ) THEN
DO 30 J = 1, N32, 32
IX = IX0
DO 20 I = I1, I2, INC
IP = IPIV( IX )
IF( IP.NE.I ) THEN
DO 10 K = J, J + 31
TEMP = A( I, K )
A( I, K ) = A( IP, K )
A( IP, K ) = TEMP
10 CONTINUE
END IF
IX = IX + INCX
20 CONTINUE
30 CONTINUE
END IF
IF( N32.NE.N ) THEN
N32 = N32 + 1
IX = IX0
DO 50 I = I1, I2, INC
IP = IPIV( IX )
IF( IP.NE.I ) THEN
DO 40 K = N32, N
TEMP = A( I, K )
A( I, K ) = A( IP, K )
A( IP, K ) = TEMP
40 CONTINUE
END IF
IX = IX + INCX
50 CONTINUE
END IF
!
RETURN
!
! End of DLASWP
!
END
SUBROUTINE DTPTRS( UPLO, TRANS, DIAG, N, NRHS, AP, B, LDB, INFO )
use MSIMSL
!
! -- LAPACK routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. Scalar Arguments ..
CHARACTER DIAG, TRANS, UPLO
INTEGER INFO, LDB, N, NRHS
! ..
! .. Array Arguments ..
DOUBLE PRECISION AP( * ), B( LDB, * )
! ..
!
! Purpose
! =======
!
! DTPTRS solves a triangular system of the form
!
! A * X = B or A**T * X = B,
!
! where A is a triangular matrix of order N stored in packed format,
! and B is an N-by-NRHS matrix. A check is made to verify that A is
! nonsingular.
!
! Arguments
! =========
!
! UPLO (input) CHARACTER*1
! = 'U': A is upper triangular;
! = 'L': A is lower triangular.
!
! TRANS (input) CHARACTER*1
! Specifies the form of the system of equations:
! = 'N': A * X = B (No transpose)
! = 'T': A**T * X = B (Transpose)
! = 'C': A**H * X = B (Conjugate transpose = Transpose)
!
! DIAG (input) CHARACTER*1
! = 'N': A is non-unit triangular;
! = 'U': A is unit triangular.
!
! N (input) INTEGER
! The order of the matrix A. N >= 0.
!
! NRHS (input) INTEGER
! The number of right hand sides, i.e., the number of columns
! of the matrix B. NRHS >= 0.
!
! AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
! The upper or lower triangular matrix A, packed columnwise in
! a linear array. The j-th column of A is stored in the array
! AP as follows:
! if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
! if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!
! B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
! On entry, the right hand side matrix B.
! On exit, if INFO = 0, the solution matrix X.
!
! LDB (input) INTEGER
! The leading dimension of the array B. LDB >= max(1,N).
!
! INFO (output) INTEGER
! = 0: successful exit
! < 0: if INFO = -i, the i-th argument had an illegal value
! > 0: if INFO = i, the i-th diagonal element of A is zero,
! indicating that the matrix is singular and the
! solutions X have not been computed.
!
! =====================================================================
!
! .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D+0 )
! ..
! .. Local Scalars ..
LOGICAL NOUNIT, UPPER
INTEGER J, JC
! ..
! .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
! ..
! .. External Subroutines ..
EXTERNAL DTPSV, XERBLA
! ..
! .. Intrinsic Functions ..
INTRINSIC MAX
! ..
! .. Executable Statements ..
!
! Test the input parameters.
!
INFO = 0
UPPER = LSAME( UPLO, 'U' )
NOUNIT = LSAME( DIAG, 'N' )
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT. &
& LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
INFO = -2
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
INFO = -3
ELSE IF( N.LT.0 ) THEN
INFO = -4
ELSE IF( NRHS.LT.0 ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -8
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DTPTRS', -INFO )
RETURN
END IF
!
! Quick return if possible
!
IF( N.EQ.0 ) &
& RETURN
!
! Check for singularity.
!
IF( NOUNIT ) THEN
IF( UPPER ) THEN
JC = 1
DO 10 INFO = 1, N
IF( AP( JC+INFO-1 ).EQ.ZERO ) &
& RETURN
JC = JC + INFO
10 CONTINUE
ELSE
JC = 1
DO 20 INFO = 1, N
IF( AP( JC ).EQ.ZERO ) &
& RETURN
JC = JC + N - INFO + 1
20 CONTINUE
END IF
END IF
INFO = 0
!
! Solve A * x = b or A' * x = b.
!
DO 30 J = 1, NRHS
CALL DTPSV( UPLO, TRANS, DIAG, N, AP, B( 1, J ), 1 )
30 END DO
!
RETURN
!
! End of DTPTRS
!
END
DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
!
! -- LAPACK auxiliary routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. 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
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
FIRST = .FALSE.
RETURN
!
! End of DLAMCH
!
END
!
!***********************************************************************
!
SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
!
! -- LAPACK auxiliary routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. 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
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
FIRST = .FALSE.
RETURN
!
! End of DLAMC1
!
END
!
!***********************************************************************
!
SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
!
! -- LAPACK auxiliary routine (version 3.1) --
! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
! November 2006
!
! .. 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
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
FIRST = .FALSE.
!**
! 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