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