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