001:       SUBROUTINE ZPTCON( N, D, E, ANORM, RCOND, RWORK, 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( * ), RWORK( * )
013:       COMPLEX*16         E( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZPTCON computes the reciprocal of the condition number (in the
020: *  1-norm) of a complex Hermitian positive definite tridiagonal matrix
021: *  using the factorization A = L*D*L**H or A = U**H*D*U computed by
022: *  ZPTTRF.
023: *
024: *  Norm(inv(A)) is computed by a direct method, and the reciprocal of
025: *  the condition number is computed as
026: *                   RCOND = 1 / (ANORM * norm(inv(A))).
027: *
028: *  Arguments
029: *  =========
030: *
031: *  N       (input) INTEGER
032: *          The order of the matrix A.  N >= 0.
033: *
034: *  D       (input) DOUBLE PRECISION array, dimension (N)
035: *          The n diagonal elements of the diagonal matrix D from the
036: *          factorization of A, as computed by ZPTTRF.
037: *
038: *  E       (input) COMPLEX*16 array, dimension (N-1)
039: *          The (n-1) off-diagonal elements of the unit bidiagonal factor
040: *          U or L from the factorization of A, as computed by ZPTTRF.
041: *
042: *  ANORM   (input) DOUBLE PRECISION
043: *          The 1-norm of the original matrix A.
044: *
045: *  RCOND   (output) DOUBLE PRECISION
046: *          The reciprocal of the condition number of the matrix A,
047: *          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
048: *          1-norm of inv(A) computed in this routine.
049: *
050: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
051: *
052: *  INFO    (output) INTEGER
053: *          = 0:  successful exit
054: *          < 0:  if INFO = -i, the i-th argument had an illegal value
055: *
056: *  Further Details
057: *  ===============
058: *
059: *  The method used is described in Nicholas J. Higham, "Efficient
060: *  Algorithms for Computing the Condition Number of a Tridiagonal
061: *  Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
062: *
063: *  =====================================================================
064: *
065: *     .. Parameters ..
066:       DOUBLE PRECISION   ONE, ZERO
067:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
068: *     ..
069: *     .. Local Scalars ..
070:       INTEGER            I, IX
071:       DOUBLE PRECISION   AINVNM
072: *     ..
073: *     .. External Functions ..
074:       INTEGER            IDAMAX
075:       EXTERNAL           IDAMAX
076: *     ..
077: *     .. External Subroutines ..
078:       EXTERNAL           XERBLA
079: *     ..
080: *     .. Intrinsic Functions ..
081:       INTRINSIC          ABS
082: *     ..
083: *     .. Executable Statements ..
084: *
085: *     Test the input arguments.
086: *
087:       INFO = 0
088:       IF( N.LT.0 ) THEN
089:          INFO = -1
090:       ELSE IF( ANORM.LT.ZERO ) THEN
091:          INFO = -4
092:       END IF
093:       IF( INFO.NE.0 ) THEN
094:          CALL XERBLA( 'ZPTCON', -INFO )
095:          RETURN
096:       END IF
097: *
098: *     Quick return if possible
099: *
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:       RWORK( 1 ) = ONE
125:       DO 20 I = 2, N
126:          RWORK( I ) = ONE + RWORK( I-1 )*ABS( E( I-1 ) )
127:    20 CONTINUE
128: *
129: *     Solve D * M(L)' * x = b.
130: *
131:       RWORK( N ) = RWORK( N ) / D( N )
132:       DO 30 I = N - 1, 1, -1
133:          RWORK( I ) = RWORK( I ) / D( I ) + RWORK( I+1 )*ABS( E( I ) )
134:    30 CONTINUE
135: *
136: *     Compute AINVNM = max(x(i)), 1<=i<=n.
137: *
138:       IX = IDAMAX( N, RWORK, 1 )
139:       AINVNM = ABS( RWORK( 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 ZPTCON
149: *
150:       END
151: