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