001:       DOUBLE PRECISION FUNCTION ZLA_PORCOND_X( UPLO, N, A, LDA, AF,
002:      $                                         LDAF, X, INFO, WORK,
003:      $                                         RWORK )
004: *
005: *     -- LAPACK routine (version 3.2.1)                                 --
006: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
007: *     -- Jason Riedy of Univ. of California Berkeley.                 --
008: *     -- April 2009                                                   --
009: *
010: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
011: *     -- Univ. of California Berkeley and NAG Ltd.                    --
012: *
013:       IMPLICIT NONE
014: *     ..
015: *     .. Scalar Arguments ..
016:       CHARACTER          UPLO
017:       INTEGER            N, LDA, LDAF, INFO
018: *     ..
019: *     .. Array Arguments ..
020:       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
021:       DOUBLE PRECISION   RWORK( * )
022: *     ..
023: *
024: *  Purpose
025: *  =======
026: *
027: *     ZLA_PORCOND_X Computes the infinity norm condition number of
028: *     op(A) * diag(X) where X is a COMPLEX*16 vector.
029: *
030: *  Arguments
031: *  =========
032: *
033: *     UPLO    (input) CHARACTER*1
034: *       = 'U':  Upper triangle of A is stored;
035: *       = 'L':  Lower triangle of A is stored.
036: *
037: *     N       (input) INTEGER
038: *     The number of linear equations, i.e., the order of the
039: *     matrix A.  N >= 0.
040: *
041: *     A       (input) COMPLEX*16 array, dimension (LDA,N)
042: *     On entry, the N-by-N matrix A.
043: *
044: *     LDA     (input) INTEGER
045: *     The leading dimension of the array A.  LDA >= max(1,N).
046: *
047: *     AF      (input) COMPLEX*16 array, dimension (LDAF,N)
048: *     The triangular factor U or L from the Cholesky factorization
049: *     A = U**T*U or A = L*L**T, as computed by ZPOTRF.
050: *
051: *     LDAF    (input) INTEGER
052: *     The leading dimension of the array AF.  LDAF >= max(1,N).
053: *
054: *     X       (input) COMPLEX*16 array, dimension (N)
055: *     The vector X in the formula op(A) * diag(X).
056: *
057: *     INFO    (output) INTEGER
058: *       = 0:  Successful exit.
059: *     i > 0:  The ith argument is invalid.
060: *
061: *     WORK    (input) COMPLEX*16 array, dimension (2*N).
062: *     Workspace.
063: *
064: *     RWORK   (input) DOUBLE PRECISION array, dimension (N).
065: *     Workspace.
066: *
067: *  =====================================================================
068: *
069: *     .. Local Scalars ..
070:       INTEGER            KASE, I, J
071:       DOUBLE PRECISION   AINVNM, ANORM, TMP
072:       LOGICAL            UP
073:       COMPLEX*16         ZDUM
074: *     ..
075: *     .. Local Arrays ..
076:       INTEGER            ISAVE( 3 )
077: *     ..
078: *     .. External Functions ..
079:       LOGICAL            LSAME
080:       EXTERNAL           LSAME
081: *     ..
082: *     .. External Subroutines ..
083:       EXTERNAL           ZLACN2, ZPOTRS, XERBLA
084: *     ..
085: *     .. Intrinsic Functions ..
086:       INTRINSIC          ABS, MAX, REAL, DIMAG
087: *     ..
088: *     .. Statement Functions ..
089:       DOUBLE PRECISION CABS1
090: *     ..
091: *     .. Statement Function Definitions ..
092:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
093: *     ..
094: *     .. Executable Statements ..
095: *
096:       ZLA_PORCOND_X = 0.0D+0
097: *
098:       INFO = 0
099:       IF( N.LT.0 ) THEN
100:          INFO = -2
101:       END IF
102:       IF( INFO.NE.0 ) THEN
103:          CALL XERBLA( 'ZLA_PORCOND_X', -INFO )
104:          RETURN
105:       END IF
106:       UP = .FALSE.
107:       IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
108: *
109: *     Compute norm of op(A)*op2(C).
110: *
111:       ANORM = 0.0D+0
112:       IF ( UP ) THEN
113:          DO I = 1, N
114:             TMP = 0.0D+0
115:             DO J = 1, I
116:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
117:             END DO
118:             DO J = I+1, N
119:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
120:             END DO
121:             RWORK( I ) = TMP
122:             ANORM = MAX( ANORM, TMP )
123:          END DO
124:       ELSE
125:          DO I = 1, N
126:             TMP = 0.0D+0
127:             DO J = 1, I
128:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
129:             END DO
130:             DO J = I+1, N
131:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
132:             END DO
133:             RWORK( I ) = TMP
134:             ANORM = MAX( ANORM, TMP )
135:          END DO
136:       END IF
137: *
138: *     Quick return if possible.
139: *
140:       IF( N.EQ.0 ) THEN
141:          ZLA_PORCOND_X = 1.0D+0
142:          RETURN
143:       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
144:          RETURN
145:       END IF
146: *
147: *     Estimate the norm of inv(op(A)).
148: *
149:       AINVNM = 0.0D+0
150: *
151:       KASE = 0
152:    10 CONTINUE
153:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
154:       IF( KASE.NE.0 ) THEN
155:          IF( KASE.EQ.2 ) THEN
156: *
157: *           Multiply by R.
158: *
159:             DO I = 1, N
160:                WORK( I ) = WORK( I ) * RWORK( I )
161:             END DO
162: *
163:             IF ( UP ) THEN
164:                CALL ZPOTRS( 'U', N, 1, AF, LDAF,
165:      $            WORK, N, INFO )
166:             ELSE
167:                CALL ZPOTRS( 'L', N, 1, AF, LDAF,
168:      $            WORK, N, INFO )
169:             ENDIF
170: *
171: *           Multiply by inv(X).
172: *
173:             DO I = 1, N
174:                WORK( I ) = WORK( I ) / X( I )
175:             END DO
176:          ELSE
177: *
178: *           Multiply by inv(X').
179: *
180:             DO I = 1, N
181:                WORK( I ) = WORK( I ) / X( I )
182:             END DO
183: *
184:             IF ( UP ) THEN
185:                CALL ZPOTRS( 'U', N, 1, AF, LDAF,
186:      $            WORK, N, INFO )
187:             ELSE
188:                CALL ZPOTRS( 'L', N, 1, AF, LDAF,
189:      $            WORK, N, INFO )
190:             END IF
191: *
192: *           Multiply by R.
193: *
194:             DO I = 1, N
195:                WORK( I ) = WORK( I ) * RWORK( I )
196:             END DO
197:          END IF
198:          GO TO 10
199:       END IF
200: *
201: *     Compute the estimate of the reciprocal condition number.
202: *
203:       IF( AINVNM .NE. 0.0D+0 )
204:      $   ZLA_PORCOND_X = 1.0D+0 / AINVNM
205: *
206:       RETURN
207: *
208:       END
209: