001:       SUBROUTINE ZGBRFS( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB,
002:      $                   IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK,
003:      $                   INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
011: *
012: *     .. Scalar Arguments ..
013:       CHARACTER          TRANS
014:       INTEGER            INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
015: *     ..
016: *     .. Array Arguments ..
017:       INTEGER            IPIV( * )
018:       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * )
019:       COMPLEX*16         AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
020:      $                   WORK( * ), X( LDX, * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  ZGBRFS improves the computed solution to a system of linear
027: *  equations when the coefficient matrix is banded, and provides
028: *  error bounds and backward error estimates for the solution.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  TRANS   (input) CHARACTER*1
034: *          Specifies the form of the system of equations:
035: *          = 'N':  A * X = B     (No transpose)
036: *          = 'T':  A**T * X = B  (Transpose)
037: *          = 'C':  A**H * X = B  (Conjugate transpose)
038: *
039: *  N       (input) INTEGER
040: *          The order of the matrix A.  N >= 0.
041: *
042: *  KL      (input) INTEGER
043: *          The number of subdiagonals within the band of A.  KL >= 0.
044: *
045: *  KU      (input) INTEGER
046: *          The number of superdiagonals within the band of A.  KU >= 0.
047: *
048: *  NRHS    (input) INTEGER
049: *          The number of right hand sides, i.e., the number of columns
050: *          of the matrices B and X.  NRHS >= 0.
051: *
052: *  AB      (input) COMPLEX*16 array, dimension (LDAB,N)
053: *          The original band matrix A, stored 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 band
064: *          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
065: *          the multipliers used during the factorization are stored in
066: *          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 ZGBTRF; for 1<=i<=N, row i of the
073: *          matrix was interchanged with row IPIV(i).
074: *
075: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
076: *          The right hand side matrix B.
077: *
078: *  LDB     (input) INTEGER
079: *          The leading dimension of the array B.  LDB >= max(1,N).
080: *
081: *  X       (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
082: *          On entry, the solution matrix X, as computed by ZGBTRS.
083: *          On exit, the improved solution matrix X.
084: *
085: *  LDX     (input) INTEGER
086: *          The leading dimension of the array X.  LDX >= max(1,N).
087: *
088: *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
089: *          The estimated forward error bound for each solution vector
090: *          X(j) (the j-th column of the solution matrix X).
091: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
092: *          is an estimated upper bound for the magnitude of the largest
093: *          element in (X(j) - XTRUE) divided by the magnitude of the
094: *          largest element in X(j).  The estimate is as reliable as
095: *          the estimate for RCOND, and is almost always a slight
096: *          overestimate of the true error.
097: *
098: *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
099: *          The componentwise relative backward error of each solution
100: *          vector X(j) (i.e., the smallest relative change in
101: *          any element of A or B that makes X(j) an exact solution).
102: *
103: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
104: *
105: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
106: *
107: *  INFO    (output) INTEGER
108: *          = 0:  successful exit
109: *          < 0:  if INFO = -i, the i-th argument had an illegal value
110: *
111: *  Internal Parameters
112: *  ===================
113: *
114: *  ITMAX is the maximum number of steps of iterative refinement.
115: *
116: *  =====================================================================
117: *
118: *     .. Parameters ..
119:       INTEGER            ITMAX
120:       PARAMETER          ( ITMAX = 5 )
121:       DOUBLE PRECISION   ZERO
122:       PARAMETER          ( ZERO = 0.0D+0 )
123:       COMPLEX*16         CONE
124:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
125:       DOUBLE PRECISION   TWO
126:       PARAMETER          ( TWO = 2.0D+0 )
127:       DOUBLE PRECISION   THREE
128:       PARAMETER          ( THREE = 3.0D+0 )
129: *     ..
130: *     .. Local Scalars ..
131:       LOGICAL            NOTRAN
132:       CHARACTER          TRANSN, TRANST
133:       INTEGER            COUNT, I, J, K, KASE, KK, NZ
134:       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
135:       COMPLEX*16         ZDUM
136: *     ..
137: *     .. Local Arrays ..
138:       INTEGER            ISAVE( 3 )
139: *     ..
140: *     .. External Subroutines ..
141:       EXTERNAL           XERBLA, ZAXPY, ZCOPY, ZGBMV, ZGBTRS, ZLACN2
142: *     ..
143: *     .. Intrinsic Functions ..
144:       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
145: *     ..
146: *     .. External Functions ..
147:       LOGICAL            LSAME
148:       DOUBLE PRECISION   DLAMCH
149:       EXTERNAL           LSAME, DLAMCH
150: *     ..
151: *     .. Statement Functions ..
152:       DOUBLE PRECISION   CABS1
153: *     ..
154: *     .. Statement Function definitions ..
155:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
156: *     ..
157: *     .. Executable Statements ..
158: *
159: *     Test the input parameters.
160: *
161:       INFO = 0
162:       NOTRAN = LSAME( TRANS, 'N' )
163:       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
164:      $    LSAME( TRANS, 'C' ) ) THEN
165:          INFO = -1
166:       ELSE IF( N.LT.0 ) THEN
167:          INFO = -2
168:       ELSE IF( KL.LT.0 ) THEN
169:          INFO = -3
170:       ELSE IF( KU.LT.0 ) THEN
171:          INFO = -4
172:       ELSE IF( NRHS.LT.0 ) THEN
173:          INFO = -5
174:       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
175:          INFO = -7
176:       ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
177:          INFO = -9
178:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
179:          INFO = -12
180:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
181:          INFO = -14
182:       END IF
183:       IF( INFO.NE.0 ) THEN
184:          CALL XERBLA( 'ZGBRFS', -INFO )
185:          RETURN
186:       END IF
187: *
188: *     Quick return if possible
189: *
190:       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
191:          DO 10 J = 1, NRHS
192:             FERR( J ) = ZERO
193:             BERR( J ) = ZERO
194:    10    CONTINUE
195:          RETURN
196:       END IF
197: *
198:       IF( NOTRAN ) THEN
199:          TRANSN = 'N'
200:          TRANST = 'C'
201:       ELSE
202:          TRANSN = 'C'
203:          TRANST = 'N'
204:       END IF
205: *
206: *     NZ = maximum number of nonzero elements in each row of A, plus 1
207: *
208:       NZ = MIN( KL+KU+2, N+1 )
209:       EPS = DLAMCH( 'Epsilon' )
210:       SAFMIN = DLAMCH( 'Safe minimum' )
211:       SAFE1 = NZ*SAFMIN
212:       SAFE2 = SAFE1 / EPS
213: *
214: *     Do for each right hand side
215: *
216:       DO 140 J = 1, NRHS
217: *
218:          COUNT = 1
219:          LSTRES = THREE
220:    20    CONTINUE
221: *
222: *        Loop until stopping criterion is satisfied.
223: *
224: *        Compute residual R = B - op(A) * X,
225: *        where op(A) = A, A**T, or A**H, depending on TRANS.
226: *
227:          CALL ZCOPY( N, B( 1, J ), 1, WORK, 1 )
228:          CALL ZGBMV( TRANS, N, N, KL, KU, -CONE, AB, LDAB, X( 1, J ), 1,
229:      $               CONE, WORK, 1 )
230: *
231: *        Compute componentwise relative backward error from formula
232: *
233: *        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
234: *
235: *        where abs(Z) is the componentwise absolute value of the matrix
236: *        or vector Z.  If the i-th component of the denominator is less
237: *        than SAFE2, then SAFE1 is added to the i-th components of the
238: *        numerator and denominator before dividing.
239: *
240:          DO 30 I = 1, N
241:             RWORK( I ) = CABS1( B( I, J ) )
242:    30    CONTINUE
243: *
244: *        Compute abs(op(A))*abs(X) + abs(B).
245: *
246:          IF( NOTRAN ) THEN
247:             DO 50 K = 1, N
248:                KK = KU + 1 - K
249:                XK = CABS1( X( K, J ) )
250:                DO 40 I = MAX( 1, K-KU ), MIN( N, K+KL )
251:                   RWORK( I ) = RWORK( I ) + CABS1( AB( KK+I, K ) )*XK
252:    40          CONTINUE
253:    50       CONTINUE
254:          ELSE
255:             DO 70 K = 1, N
256:                S = ZERO
257:                KK = KU + 1 - K
258:                DO 60 I = MAX( 1, K-KU ), MIN( N, K+KL )
259:                   S = S + CABS1( AB( KK+I, K ) )*CABS1( X( I, J ) )
260:    60          CONTINUE
261:                RWORK( K ) = RWORK( K ) + S
262:    70       CONTINUE
263:          END IF
264:          S = ZERO
265:          DO 80 I = 1, N
266:             IF( RWORK( I ).GT.SAFE2 ) THEN
267:                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
268:             ELSE
269:                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
270:      $             ( RWORK( I )+SAFE1 ) )
271:             END IF
272:    80    CONTINUE
273:          BERR( J ) = S
274: *
275: *        Test stopping criterion. Continue iterating if
276: *           1) The residual BERR(J) is larger than machine epsilon, and
277: *           2) BERR(J) decreased by at least a factor of 2 during the
278: *              last iteration, and
279: *           3) At most ITMAX iterations tried.
280: *
281:          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
282:      $       COUNT.LE.ITMAX ) THEN
283: *
284: *           Update solution and try again.
285: *
286:             CALL ZGBTRS( TRANS, N, KL, KU, 1, AFB, LDAFB, IPIV, WORK, N,
287:      $                   INFO )
288:             CALL ZAXPY( N, CONE, WORK, 1, X( 1, J ), 1 )
289:             LSTRES = BERR( J )
290:             COUNT = COUNT + 1
291:             GO TO 20
292:          END IF
293: *
294: *        Bound error from formula
295: *
296: *        norm(X - XTRUE) / norm(X) .le. FERR =
297: *        norm( abs(inv(op(A)))*
298: *           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
299: *
300: *        where
301: *          norm(Z) is the magnitude of the largest component of Z
302: *          inv(op(A)) is the inverse of op(A)
303: *          abs(Z) is the componentwise absolute value of the matrix or
304: *             vector Z
305: *          NZ is the maximum number of nonzeros in any row of A, plus 1
306: *          EPS is machine epsilon
307: *
308: *        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
309: *        is incremented by SAFE1 if the i-th component of
310: *        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
311: *
312: *        Use ZLACN2 to estimate the infinity-norm of the matrix
313: *           inv(op(A)) * diag(W),
314: *        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
315: *
316:          DO 90 I = 1, N
317:             IF( RWORK( I ).GT.SAFE2 ) THEN
318:                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
319:             ELSE
320:                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
321:      $                      SAFE1
322:             END IF
323:    90    CONTINUE
324: *
325:          KASE = 0
326:   100    CONTINUE
327:          CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
328:          IF( KASE.NE.0 ) THEN
329:             IF( KASE.EQ.1 ) THEN
330: *
331: *              Multiply by diag(W)*inv(op(A)**H).
332: *
333:                CALL ZGBTRS( TRANST, N, KL, KU, 1, AFB, LDAFB, IPIV,
334:      $                      WORK, N, INFO )
335:                DO 110 I = 1, N
336:                   WORK( I ) = RWORK( I )*WORK( I )
337:   110          CONTINUE
338:             ELSE
339: *
340: *              Multiply by inv(op(A))*diag(W).
341: *
342:                DO 120 I = 1, N
343:                   WORK( I ) = RWORK( I )*WORK( I )
344:   120          CONTINUE
345:                CALL ZGBTRS( TRANSN, N, KL, KU, 1, AFB, LDAFB, IPIV,
346:      $                      WORK, N, INFO )
347:             END IF
348:             GO TO 100
349:          END IF
350: *
351: *        Normalize error.
352: *
353:          LSTRES = ZERO
354:          DO 130 I = 1, N
355:             LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
356:   130    CONTINUE
357:          IF( LSTRES.NE.ZERO )
358:      $      FERR( J ) = FERR( J ) / LSTRES
359: *
360:   140 CONTINUE
361: *
362:       RETURN
363: *
364: *     End of ZGBRFS
365: *
366:       END
367: