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