Annotation of rpl/lapack/lapack/dptcon.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: INTEGER INFO, N
! 10: DOUBLE PRECISION ANORM, RCOND
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION D( * ), E( * ), WORK( * )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * DPTCON computes the reciprocal of the condition number (in the
! 20: * 1-norm) of a real symmetric positive definite tridiagonal matrix
! 21: * using the factorization A = L*D*L**T or A = U**T*D*U computed by
! 22: * DPTTRF.
! 23: *
! 24: * Norm(inv(A)) is computed by a direct method, and the reciprocal of
! 25: * the condition number is computed as
! 26: * RCOND = 1 / (ANORM * norm(inv(A))).
! 27: *
! 28: * Arguments
! 29: * =========
! 30: *
! 31: * N (input) INTEGER
! 32: * The order of the matrix A. N >= 0.
! 33: *
! 34: * D (input) DOUBLE PRECISION array, dimension (N)
! 35: * The n diagonal elements of the diagonal matrix D from the
! 36: * factorization of A, as computed by DPTTRF.
! 37: *
! 38: * E (input) DOUBLE PRECISION array, dimension (N-1)
! 39: * The (n-1) off-diagonal elements of the unit bidiagonal factor
! 40: * U or L from the factorization of A, as computed by DPTTRF.
! 41: *
! 42: * ANORM (input) DOUBLE PRECISION
! 43: * The 1-norm of the original matrix A.
! 44: *
! 45: * RCOND (output) DOUBLE PRECISION
! 46: * The reciprocal of the condition number of the matrix A,
! 47: * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
! 48: * 1-norm of inv(A) computed in this routine.
! 49: *
! 50: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
! 51: *
! 52: * INFO (output) INTEGER
! 53: * = 0: successful exit
! 54: * < 0: if INFO = -i, the i-th argument had an illegal value
! 55: *
! 56: * Further Details
! 57: * ===============
! 58: *
! 59: * The method used is described in Nicholas J. Higham, "Efficient
! 60: * Algorithms for Computing the Condition Number of a Tridiagonal
! 61: * Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
! 62: *
! 63: * =====================================================================
! 64: *
! 65: * .. Parameters ..
! 66: DOUBLE PRECISION ONE, ZERO
! 67: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 68: * ..
! 69: * .. Local Scalars ..
! 70: INTEGER I, IX
! 71: DOUBLE PRECISION AINVNM
! 72: * ..
! 73: * .. External Functions ..
! 74: INTEGER IDAMAX
! 75: EXTERNAL IDAMAX
! 76: * ..
! 77: * .. External Subroutines ..
! 78: EXTERNAL XERBLA
! 79: * ..
! 80: * .. Intrinsic Functions ..
! 81: INTRINSIC ABS
! 82: * ..
! 83: * .. Executable Statements ..
! 84: *
! 85: * Test the input arguments.
! 86: *
! 87: INFO = 0
! 88: IF( N.LT.0 ) THEN
! 89: INFO = -1
! 90: ELSE IF( ANORM.LT.ZERO ) THEN
! 91: INFO = -4
! 92: END IF
! 93: IF( INFO.NE.0 ) THEN
! 94: CALL XERBLA( 'DPTCON', -INFO )
! 95: RETURN
! 96: END IF
! 97: *
! 98: * Quick return if possible
! 99: *
! 100: RCOND = ZERO
! 101: IF( N.EQ.0 ) THEN
! 102: RCOND = ONE
! 103: RETURN
! 104: ELSE IF( ANORM.EQ.ZERO ) THEN
! 105: RETURN
! 106: END IF
! 107: *
! 108: * Check that D(1:N) is positive.
! 109: *
! 110: DO 10 I = 1, N
! 111: IF( D( I ).LE.ZERO )
! 112: $ RETURN
! 113: 10 CONTINUE
! 114: *
! 115: * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
! 116: *
! 117: * m(i,j) = abs(A(i,j)), i = j,
! 118: * m(i,j) = -abs(A(i,j)), i .ne. j,
! 119: *
! 120: * and e = [ 1, 1, ..., 1 ]'. Note M(A) = M(L)*D*M(L)'.
! 121: *
! 122: * Solve M(L) * x = e.
! 123: *
! 124: WORK( 1 ) = ONE
! 125: DO 20 I = 2, N
! 126: WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
! 127: 20 CONTINUE
! 128: *
! 129: * Solve D * M(L)' * x = b.
! 130: *
! 131: WORK( N ) = WORK( N ) / D( N )
! 132: DO 30 I = N - 1, 1, -1
! 133: WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
! 134: 30 CONTINUE
! 135: *
! 136: * Compute AINVNM = max(x(i)), 1<=i<=n.
! 137: *
! 138: IX = IDAMAX( N, WORK, 1 )
! 139: AINVNM = ABS( WORK( IX ) )
! 140: *
! 141: * Compute the reciprocal condition number.
! 142: *
! 143: IF( AINVNM.NE.ZERO )
! 144: $ RCOND = ( ONE / AINVNM ) / ANORM
! 145: *
! 146: RETURN
! 147: *
! 148: * End of DPTCON
! 149: *
! 150: END
CVSweb interface <joel.bertrand@systella.fr>