001:       REAL FUNCTION CLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV, C,
002:      $                             CAPPLY, 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 Aguments ..
015:       CHARACTER          TRANS
016:       LOGICAL            CAPPLY
017:       INTEGER            N, LDA, LDAF, INFO
018: *     ..
019: *     .. Array Arguments ..
020:       INTEGER            IPIV( * )
021:       COMPLEX            A( LDA, * ), AF( LDAF, * ), WORK( * )
022:       REAL               C( * ), RWORK( * )
023: *     ..
024: *
025: *  Purpose
026: *  =======
027: * 
028: *     CLA_GERCOND_C computes the infinity norm condition number of
029: *     op(A) * inv(diag(C)) where C is a REAL 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 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 array, dimension (LDAF,N)
051: *     The factors L and U from the factorization
052: *     A = P*L*U as computed by CGETRF.
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 CGETRF; row i of the matrix was interchanged
060: *     with row IPIV(i).
061: *
062: *     C       (input) REAL array, dimension (N)
063: *     The vector C in the formula op(A) * inv(diag(C)).
064: *
065: *     CAPPLY  (input) LOGICAL
066: *     If .TRUE. then access the vector C in the formula above.
067: *
068: *     INFO    (output) INTEGER
069: *       = 0:  Successful exit.
070: *     i > 0:  The ith argument is invalid.
071: *
072: *     WORK    (input) COMPLEX array, dimension (2*N).
073: *     Workspace.
074: *
075: *     RWORK   (input) REAL array, dimension (N).
076: *     Workspace.
077: *
078: *  =====================================================================
079: *
080: *     .. Local Scalars ..
081:       LOGICAL            NOTRANS
082:       INTEGER            KASE, I, J
083:       REAL               AINVNM, ANORM, TMP
084:       COMPLEX            ZDUM
085: *     ..
086: *     .. Local Arrays ..
087:       INTEGER            ISAVE( 3 )
088: *     ..
089: *     .. External Functions ..
090:       LOGICAL            LSAME
091:       EXTERNAL           LSAME
092: *     ..
093: *     .. External Subroutines ..
094:       EXTERNAL           CLACN2, CGETRS, XERBLA
095: *     ..
096: *     .. Intrinsic Functions ..
097:       INTRINSIC          ABS, MAX, REAL, AIMAG
098: *     ..
099: *     .. Statement Functions ..
100:       REAL               CABS1
101: *     ..
102: *     .. Statement Function Definitions ..
103:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
104: *     ..
105: *     .. Executable Statements ..
106:       CLA_GERCOND_C = 0.0E+0
107: *
108:       INFO = 0
109:       NOTRANS = LSAME( TRANS, 'N' )
110:       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
111:      $     LSAME( TRANS, 'C' ) ) THEN
112:       ELSE IF( N.LT.0 ) THEN
113:          INFO = -2
114:       END IF
115:       IF( INFO.NE.0 ) THEN
116:          CALL XERBLA( 'CLA_GERCOND_C', -INFO )
117:          RETURN
118:       END IF
119: *
120: *     Compute norm of op(A)*op2(C).
121: *
122:       ANORM = 0.0E+0
123:       IF ( NOTRANS ) THEN
124:          DO I = 1, N
125:             TMP = 0.0E+0
126:             IF ( CAPPLY ) THEN
127:                DO J = 1, N
128:                   TMP = TMP + CABS1( A( I, J ) ) / C( J )
129:                END DO
130:             ELSE
131:                DO J = 1, N
132:                   TMP = TMP + CABS1( A( I, J ) )
133:                END DO
134:             END IF
135:             RWORK( I ) = TMP
136:             ANORM = MAX( ANORM, TMP )
137:          END DO
138:       ELSE
139:          DO I = 1, N
140:             TMP = 0.0E+0
141:             IF ( CAPPLY ) THEN
142:                DO J = 1, N
143:                   TMP = TMP + CABS1( A( J, I ) ) / C( J )
144:                END DO
145:             ELSE
146:                DO J = 1, N
147:                   TMP = TMP + CABS1( A( J, I ) )
148:                END DO
149:             END IF
150:             RWORK( I ) = TMP
151:             ANORM = MAX( ANORM, TMP )
152:          END DO
153:       END IF
154: *
155: *     Quick return if possible.
156: *
157:       IF( N.EQ.0 ) THEN
158:          CLA_GERCOND_C = 1.0E+0
159:          RETURN
160:       ELSE IF( ANORM .EQ. 0.0E+0 ) THEN
161:          RETURN
162:       END IF
163: *
164: *     Estimate the norm of inv(op(A)).
165: *
166:       AINVNM = 0.0E+0
167: *
168:       KASE = 0
169:    10 CONTINUE
170:       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
171:       IF( KASE.NE.0 ) THEN
172:          IF( KASE.EQ.2 ) THEN
173: *
174: *           Multiply by R.
175: *
176:             DO I = 1, N
177:                WORK( I ) = WORK( I ) * RWORK( I )
178:             END DO
179: *
180:             IF (NOTRANS) THEN
181:                CALL CGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
182:      $            WORK, N, INFO )
183:             ELSE
184:                CALL CGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
185:      $            WORK, N, INFO )
186:             ENDIF
187: *
188: *           Multiply by inv(C).
189: *
190:             IF ( CAPPLY ) THEN
191:                DO I = 1, N
192:                   WORK( I ) = WORK( I ) * C( I )
193:                END DO
194:             END IF
195:          ELSE
196: *
197: *           Multiply by inv(C').
198: *
199:             IF ( CAPPLY ) THEN
200:                DO I = 1, N
201:                   WORK( I ) = WORK( I ) * C( I )
202:                END DO
203:             END IF
204: *
205:             IF ( NOTRANS ) THEN
206:                CALL CGETRS( 'Conjugate transpose', N, 1, AF, LDAF, IPIV,
207:      $            WORK, N, INFO )
208:             ELSE
209:                CALL CGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
210:      $            WORK, N, INFO )
211:             END IF
212: *
213: *           Multiply by R.
214: *
215:             DO I = 1, N
216:                WORK( I ) = WORK( I ) * RWORK( I )
217:             END DO
218:          END IF
219:          GO TO 10
220:       END IF
221: *
222: *     Compute the estimate of the reciprocal condition number.
223: *
224:       IF( AINVNM .NE. 0.0E+0 )
225:      $   CLA_GERCOND_C = 1.0E+0 / AINVNM
226: *
227:       RETURN
228: *
229:       END
230: