001:       SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, N
009:       DOUBLE PRECISION   ANORM, RCOND
010: *     ..
011: *     .. Array Arguments ..
012:       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  DPTCON computes the reciprocal of the condition number (in the
019: *  1-norm) of a real symmetric positive definite tridiagonal matrix
020: *  using the factorization A = L*D*L**T or A = U**T*D*U computed by
021: *  DPTTRF.
022: *
023: *  Norm(inv(A)) is computed by a direct method, and the reciprocal of
024: *  the condition number is computed as
025: *               RCOND = 1 / (ANORM * norm(inv(A))).
026: *
027: *  Arguments
028: *  =========
029: *
030: *  N       (input) INTEGER
031: *          The order of the matrix A.  N >= 0.
032: *
033: *  D       (input) DOUBLE PRECISION array, dimension (N)
034: *          The n diagonal elements of the diagonal matrix D from the
035: *          factorization of A, as computed by DPTTRF.
036: *
037: *  E       (input) DOUBLE PRECISION array, dimension (N-1)
038: *          The (n-1) off-diagonal elements of the unit bidiagonal factor
039: *          U or L from the factorization of A,  as computed by DPTTRF.
040: *
041: *  ANORM   (input) DOUBLE PRECISION
042: *          The 1-norm of the original matrix A.
043: *
044: *  RCOND   (output) DOUBLE PRECISION
045: *          The reciprocal of the condition number of the matrix A,
046: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
047: *          1-norm of inv(A) computed in this routine.
048: *
049: *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
050: *
051: *  INFO    (output) INTEGER
052: *          = 0:  successful exit
053: *          < 0:  if INFO = -i, the i-th argument had an illegal value
054: *
055: *  Further Details
056: *  ===============
057: *
058: *  The method used is described in Nicholas J. Higham, "Efficient
059: *  Algorithms for Computing the Condition Number of a Tridiagonal
060: *  Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
061: *
062: *  =====================================================================
063: *
064: *     .. Parameters ..
065:       DOUBLE PRECISION   ONE, ZERO
066:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
067: *     ..
068: *     .. Local Scalars ..
069:       INTEGER            I, IX
070:       DOUBLE PRECISION   AINVNM
071: *     ..
072: *     .. External Functions ..
073:       INTEGER            IDAMAX
074:       EXTERNAL           IDAMAX
075: *     ..
076: *     .. External Subroutines ..
077:       EXTERNAL           XERBLA
078: *     ..
079: *     .. Intrinsic Functions ..
080:       INTRINSIC          ABS
081: *     ..
082: *     .. Executable Statements ..
083: *
084: *     Test the input arguments.
085: *
086:       INFO = 0
087:       IF( N.LT.0 ) THEN
088:          INFO = -1
089:       ELSE IF( ANORM.LT.ZERO ) THEN
090:          INFO = -4
091:       END IF
092:       IF( INFO.NE.0 ) THEN
093:          CALL XERBLA( 'DPTCON', -INFO )
094:          RETURN
095:       END IF
096: *
097: *     Quick return if possible
098: *
099:       RCOND = ZERO
100:       IF( N.EQ.0 ) THEN
101:          RCOND = ONE
102:          RETURN
103:       ELSE IF( ANORM.EQ.ZERO ) THEN
104:          RETURN
105:       END IF
106: *
107: *     Check that D(1:N) is positive.
108: *
109:       DO 10 I = 1, N
110:          IF( D( I ).LE.ZERO )
111:      $      RETURN
112:    10 CONTINUE
113: *
114: *     Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
115: *
116: *        m(i,j) =  abs(A(i,j)), i = j,
117: *        m(i,j) = -abs(A(i,j)), i .ne. j,
118: *
119: *     and e = [ 1, 1, ..., 1 ]'.  Note M(A) = M(L)*D*M(L)'.
120: *
121: *     Solve M(L) * x = e.
122: *
123:       WORK( 1 ) = ONE
124:       DO 20 I = 2, N
125:          WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
126:    20 CONTINUE
127: *
128: *     Solve D * M(L)' * x = b.
129: *
130:       WORK( N ) = WORK( N ) / D( N )
131:       DO 30 I = N - 1, 1, -1
132:          WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
133:    30 CONTINUE
134: *
135: *     Compute AINVNM = max(x(i)), 1<=i<=n.
136: *
137:       IX = IDAMAX( N, WORK, 1 )
138:       AINVNM = ABS( WORK( IX ) )
139: *
140: *     Compute the reciprocal condition number.
141: *
142:       IF( AINVNM.NE.ZERO )
143:      $   RCOND = ( ONE / AINVNM ) / ANORM
144: *
145:       RETURN
146: *
147: *     End of DPTCON
148: *
149:       END
150: