001:       REAL FUNCTION CLA_GBRCOND_C( TRANS, N, KL, KU, AB, LDAB, AFB,
002:      $                             LDAFB, IPIV, C, CAPPLY, 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          TRANS
017:       LOGICAL            CAPPLY
018:       INTEGER            N, KL, KU, KD, KE, LDAB, LDAFB, INFO
019: *     ..
020: *     .. Array Arguments ..
021:       INTEGER            IPIV( * )
022:       COMPLEX            AB( LDAB, * ), AFB( LDAFB, * ), WORK( * )
023:       REAL               C( * ), RWORK( * )
024: *     ..
025: *
026: *  Purpose
027: *  =======
028: *
029: *     CLA_GBRCOND_C Computes the infinity norm condition number of
030: *     op(A) * inv(diag(C)) where C is a REAL vector.
031: *
032: *  Arguments
033: *  =========
034: *
035: *     TRANS   (input) CHARACTER*1
036: *     Specifies the form of the system of equations:
037: *       = 'N':  A * X = B     (No transpose)
038: *       = 'T':  A**T * X = B  (Transpose)
039: *       = 'C':  A**H * X = B  (Conjugate Transpose = Transpose)
040: *
041: *     N       (input) INTEGER
042: *     The number of linear equations, i.e., the order of the
043: *     matrix A.  N >= 0.
044: *
045: *     KL      (input) INTEGER
046: *     The number of subdiagonals within the band of A.  KL >= 0.
047: *
048: *     KU      (input) INTEGER
049: *     The number of superdiagonals within the band of A.  KU >= 0.
050: *
051: *     AB      (input) COMPLEX array, dimension (LDAB,N)
052: *     On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
053: *     The j-th column of A is stored in the j-th column of the
054: *     array AB as follows:
055: *     AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
056: *
057: *     LDAB    (input) INTEGER
058: *     The leading dimension of the array AB.  LDAB >= KL+KU+1.
059: *
060: *     AFB     (input) COMPLEX array, dimension (LDAFB,N)
061: *     Details of the LU factorization of the band matrix A, as
062: *     computed by CGBTRF.  U is stored as an upper triangular
063: *     band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
064: *     and the multipliers used during the factorization are stored
065: *     in rows KL+KU+2 to 2*KL+KU+1.
066: *
067: *     LDAFB   (input) INTEGER
068: *     The leading dimension of the array AFB.  LDAFB >= 2*KL+KU+1.
069: *
070: *     IPIV    (input) INTEGER array, dimension (N)
071: *     The pivot indices from the factorization A = P*L*U
072: *     as computed by CGBTRF; row i of the matrix was interchanged
073: *     with row IPIV(i).
074: *
075: *     C       (input) REAL array, dimension (N)
076: *     The vector C in the formula op(A) * inv(diag(C)).
077: *
078: *     CAPPLY  (input) LOGICAL
079: *     If .TRUE. then access the vector C in the formula above.
080: *
081: *     INFO    (output) INTEGER
082: *       = 0:  Successful exit.
083: *     i > 0:  The ith argument is invalid.
084: *
085: *     WORK    (input) COMPLEX array, dimension (2*N).
086: *     Workspace.
087: *
088: *     RWORK   (input) REAL array, dimension (N).
089: *     Workspace.
090: *
091: *  =====================================================================
092: *
093: *     .. Local Scalars ..
094:       LOGICAL            NOTRANS
095:       INTEGER            KASE, I, J
096:       REAL               AINVNM, ANORM, TMP
097:       COMPLEX            ZDUM
098: *     ..
099: *     .. Local Arrays ..
100:       INTEGER            ISAVE( 3 )
101: *     ..
102: *     .. External Functions ..
103:       LOGICAL            LSAME
104:       EXTERNAL           LSAME
105: *     ..
106: *     .. External Subroutines ..
107:       EXTERNAL           CLACN2, CGBTRS, XERBLA
108: *     ..
109: *     .. Intrinsic Functions ..
110:       INTRINSIC          ABS, MAX
111: *     ..
112: *     .. Statement Functions ..
113:       REAL               CABS1
114: *     ..
115: *     .. Statement Function Definitions ..
116:       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
117: *     ..
118: *     .. Executable Statements ..
119:       CLA_GBRCOND_C = 0.0E+0
120: *
121:       INFO = 0
122:       NOTRANS = LSAME( TRANS, 'N' )
123:       IF ( .NOT. NOTRANS .AND. .NOT. LSAME( TRANS, 'T' ) .AND. .NOT.
124:      $     LSAME( TRANS, 'C' ) ) THEN
125:          INFO = -1
126:       ELSE IF( N.LT.0 ) THEN
127:          INFO = -2
128:       ELSE IF( KL.LT.0 .OR. KL.GT.N-1 ) THEN
129:          INFO = -3
130:       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
131:          INFO = -4
132:       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
133:          INFO = -6
134:       ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
135:          INFO = -8
136:       END IF
137:       IF( INFO.NE.0 ) THEN
138:          CALL XERBLA( 'CLA_GBRCOND_C', -INFO )
139:          RETURN
140:       END IF
141: *
142: *     Compute norm of op(A)*op2(C).
143: *
144:       ANORM = 0.0E+0
145:       KD = KU + 1
146:       KE = KL + 1
147:       IF ( NOTRANS ) THEN
148:          DO I = 1, N
149:             TMP = 0.0E+0
150:             IF ( CAPPLY ) THEN
151:                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
152:                   TMP = TMP + CABS1( AB( KD+I-J, J ) ) / C( J )
153:                END DO
154:             ELSE
155:                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
156:                   TMP = TMP + CABS1( AB( KD+I-J, J ) )
157:                END DO
158:             END IF
159:             RWORK( I ) = TMP
160:             ANORM = MAX( ANORM, TMP )
161:          END DO
162:       ELSE
163:          DO I = 1, N
164:             TMP = 0.0E+0
165:             IF ( CAPPLY ) THEN
166:                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
167:                   TMP = TMP + CABS1( AB( KE-I+J, I ) ) / C( J )
168:                END DO
169:             ELSE
170:                DO J = MAX( I-KL, 1 ), MIN( I+KU, N )
171:                   TMP = TMP + CABS1( AB( KE-I+J, I ) )
172:                END DO
173:             END IF
174:             RWORK( I ) = TMP
175:             ANORM = MAX( ANORM, TMP )
176:          END DO
177:       END IF
178: *
179: *     Quick return if possible.
180: *
181:       IF( N.EQ.0 ) THEN
182:          CLA_GBRCOND_C = 1.0E+0
183:          RETURN
184:       ELSE IF( ANORM .EQ. 0.0E+0 ) THEN
185:          RETURN
186:       END IF
187: *
188: *     Estimate the norm of inv(op(A)).
189: *
190:       AINVNM = 0.0E+0
191: *
192:       KASE = 0
193:    10 CONTINUE
194:       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
195:       IF( KASE.NE.0 ) THEN
196:          IF( KASE.EQ.2 ) THEN
197: *
198: *           Multiply by R.
199: *
200:             DO I = 1, N
201:                WORK( I ) = WORK( I ) * RWORK( I )
202:             END DO
203: *
204:             IF ( NOTRANS ) THEN
205:                CALL CGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
206:      $              IPIV, WORK, N, INFO )
207:             ELSE
208:                CALL CGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
209:      $              LDAFB, IPIV, WORK, N, INFO )
210:             ENDIF
211: *
212: *           Multiply by inv(C).
213: *
214:             IF ( CAPPLY ) THEN
215:                DO I = 1, N
216:                   WORK( I ) = WORK( I ) * C( I )
217:                END DO
218:             END IF
219:          ELSE
220: *
221: *           Multiply by inv(C').
222: *
223:             IF ( CAPPLY ) THEN
224:                DO I = 1, N
225:                   WORK( I ) = WORK( I ) * C( I )
226:                END DO
227:             END IF
228: *
229:             IF ( NOTRANS ) THEN
230:                CALL CGBTRS( 'Conjugate transpose', N, KL, KU, 1, AFB,
231:      $              LDAFB, IPIV,  WORK, N, INFO )
232:             ELSE
233:                CALL CGBTRS( 'No transpose', N, KL, KU, 1, AFB, LDAFB,
234:      $              IPIV, WORK, N, INFO )
235:             END IF
236: *
237: *           Multiply by R.
238: *
239:             DO I = 1, N
240:                WORK( I ) = WORK( I ) * RWORK( I )
241:             END DO
242:          END IF
243:          GO TO 10
244:       END IF
245: *
246: *     Compute the estimate of the reciprocal condition number.
247: *
248:       IF( AINVNM .NE. 0.0E+0 )
249:      $   CLA_GBRCOND_C = 1.0E+0 / AINVNM
250: *
251:       RETURN
252: *
253:       END
254: