001:       DOUBLE PRECISION FUNCTION ZLA_GERCOND_X( TRANS, N, A, LDA, AF,
002:      $                                         LDAF, IPIV, X, INFO,
003:      $                                         WORK, 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          TRANS
017:       INTEGER            N, LDA, LDAF, INFO
018: *     ..
019: *     .. Array Arguments ..
020:       INTEGER            IPIV( * )
021:       COMPLEX*16         A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
022:       DOUBLE PRECISION   RWORK( * )
023: *     ..
024: *
025: *  Purpose
026: *  =======
027: *
028: *     ZLA_GERCOND_X computes the infinity norm condition number of
029: *     op(A) * diag(X) where X is a COMPLEX*16 vector.
030: *
031: *  Arguments
032: *  =========
033: *
034: *     TRANS   (input) CHARACTER*1
035: *     Specifies the form of the system of equations:
036: *       = 'N':  A * X = B     (No transpose)
037: *       = 'T':  A**T * X = B  (Transpose)
038: *       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
039: *
040: *     N       (input) INTEGER
041: *     The number of linear equations, i.e., the order of the
042: *     matrix A.  N >= 0.
043: *
044: *     A       (input) COMPLEX*16 array, dimension (LDA,N)
045: *     On entry, the N-by-N matrix A.
046: *
047: *     LDA     (input) INTEGER
048: *     The leading dimension of the array A.  LDA >= max(1,N).
049: *
050: *     AF      (input) COMPLEX*16 array, dimension (LDAF,N)
051: *     The factors L and U from the factorization
052: *     A = P*L*U as computed by ZGETRF.
053: *
054: *     LDAF    (input) INTEGER
055: *     The leading dimension of the array AF.  LDAF >= max(1,N).
056: *
057: *     IPIV    (input) INTEGER array, dimension (N)
058: *     The pivot indices from the factorization A = P*L*U
059: *     as computed by ZGETRF; row i of the matrix was interchanged
060: *     with row IPIV(i).
061: *
062: *     X       (input) COMPLEX*16 array, dimension (N)
063: *     The vector X in the formula op(A) * diag(X).
064: *
065: *     INFO    (output) INTEGER
066: *       = 0:  Successful exit.
067: *     i > 0:  The ith argument is invalid.
068: *
069: *     WORK    (input) COMPLEX*16 array, dimension (2*N).
070: *     Workspace.
071: *
072: *     RWORK   (input) DOUBLE PRECISION array, dimension (N).
073: *     Workspace.
074: *
075: *  =====================================================================
076: *
077: *     .. Local Scalars ..
078:       LOGICAL            NOTRANS
079:       INTEGER            KASE
080:       DOUBLE PRECISION   AINVNM, ANORM, TMP
081:       INTEGER            I, J
082:       COMPLEX*16         ZDUM
083: *     ..
084: *     .. Local Arrays ..
085:       INTEGER            ISAVE( 3 )
086: *     ..
087: *     .. External Functions ..
088:       LOGICAL            LSAME
089:       EXTERNAL           LSAME
090: *     ..
091: *     .. External Subroutines ..
092:       EXTERNAL           ZLACN2, ZGETRS, XERBLA
093: *     ..
094: *     .. Intrinsic Functions ..
095:       INTRINSIC          ABS, MAX, REAL, DIMAG
096: *     ..
097: *     .. Statement Functions ..
098:       DOUBLE PRECISION   CABS1
099: *     ..
100: *     .. Statement Function Definitions ..
101:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
102: *     ..
103: *     .. Executable Statements ..
104: *
105:       ZLA_GERCOND_X = 0.0D+0
106: *
107:       INFO = 0
108:       NOTRANS = LSAME( TRANS, 'N' )
109:       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
110:      $     LSAME( TRANS, 'C' ) ) THEN
111:          INFO = -1
112:       ELSE IF( N.LT.0 ) THEN
113:          INFO = -2
114:       END IF
115:       IF( INFO.NE.0 ) THEN
116:          CALL XERBLA( 'ZLA_GERCOND_X', -INFO )
117:          RETURN
118:       END IF
119: *
120: *     Compute norm of op(A)*op2(C).
121: *
122:       ANORM = 0.0D+0
123:       IF ( NOTRANS ) THEN
124:          DO I = 1, N
125:             TMP = 0.0D+0
126:             DO J = 1, N
127:                TMP = TMP + CABS1( A( I, J ) * X( J ) )
128:             END DO
129:             RWORK( I ) = TMP
130:             ANORM = MAX( ANORM, TMP )
131:          END DO
132:       ELSE
133:          DO I = 1, N
134:             TMP = 0.0D+0
135:             DO J = 1, N
136:                TMP = TMP + CABS1( A( J, I ) * X( J ) )
137:             END DO
138:             RWORK( I ) = TMP
139:             ANORM = MAX( ANORM, TMP )
140:          END DO
141:       END IF
142: *
143: *     Quick return if possible.
144: *
145:       IF( N.EQ.0 ) THEN
146:          ZLA_GERCOND_X = 1.0D+0
147:          RETURN
148:       ELSE IF( ANORM .EQ. 0.0D+0 ) THEN
149:          RETURN
150:       END IF
151: *
152: *     Estimate the norm of inv(op(A)).
153: *
154:       AINVNM = 0.0D+0
155: *
156:       KASE = 0
157:    10 CONTINUE
158:       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
159:       IF( KASE.NE.0 ) THEN
160:          IF( KASE.EQ.2 ) THEN
161: *           Multiply by R.
162:             DO I = 1, N
163:                WORK( I ) = WORK( I ) * RWORK( I )
164:             END DO
165: *
166:             IF ( NOTRANS ) THEN
167:                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
168:      $            WORK, N, INFO )
169:             ELSE
170:                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
171:      $            WORK, N, INFO )
172:             ENDIF
173: *
174: *           Multiply by inv(X).
175: *
176:             DO I = 1, N
177:                WORK( I ) = WORK( I ) / X( I )
178:             END DO
179:          ELSE
180: *
181: *           Multiply by inv(X').
182: *
183:             DO I = 1, N
184:                WORK( I ) = WORK( I ) / X( I )
185:             END DO
186: *
187:             IF ( NOTRANS ) THEN
188:                CALL ZGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
189:      $            WORK, N, INFO )
190:             ELSE
191:                CALL ZGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
192:      $            WORK, N, INFO )
193:             END IF
194: *
195: *           Multiply by R.
196: *
197:             DO I = 1, N
198:                WORK( I ) = WORK( I ) * RWORK( I )
199:             END DO
200:          END IF
201:          GO TO 10
202:       END IF
203: *
204: *     Compute the estimate of the reciprocal condition number.
205: *
206:       IF( AINVNM .NE. 0.0D+0 )
207:      $   ZLA_GERCOND_X = 1.0D+0 / AINVNM
208: *
209:       RETURN
210: *
211:       END
212: