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