LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zgbsvx.f
Go to the documentation of this file.
1 *> \brief <b> ZGBSVX computes the solution to system of linear equations A * X = B for GB matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGBSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbsvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbsvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbsvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGBSVX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB,
22 * LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX,
23 * RCOND, FERR, BERR, WORK, RWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER EQUED, FACT, TRANS
27 * INTEGER INFO, KL, KU, LDAB, LDAFB, LDB, LDX, N, NRHS
28 * DOUBLE PRECISION RCOND
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IPIV( * )
32 * DOUBLE PRECISION BERR( * ), C( * ), FERR( * ), R( * ),
33 * $ RWORK( * )
34 * COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
35 * $ WORK( * ), X( LDX, * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> ZGBSVX uses the LU factorization to compute the solution to a complex
45 *> system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
46 *> where A is a band matrix of order N with KL subdiagonals and KU
47 *> superdiagonals, and X and B are N-by-NRHS matrices.
48 *>
49 *> Error bounds on the solution and a condition estimate are also
50 *> provided.
51 *> \endverbatim
52 *
53 *> \par Description:
54 * =================
55 *>
56 *> \verbatim
57 *>
58 *> The following steps are performed by this subroutine:
59 *>
60 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate
61 *> the system:
62 *> TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
63 *> TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
64 *> TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
65 *> Whether or not the system will be equilibrated depends on the
66 *> scaling of the matrix A, but if equilibration is used, A is
67 *> overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
68 *> or diag(C)*B (if TRANS = 'T' or 'C').
69 *>
70 *> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the
71 *> matrix A (after equilibration if FACT = 'E') as
72 *> A = L * U,
73 *> where L is a product of permutation and unit lower triangular
74 *> matrices with KL subdiagonals, and U is upper triangular with
75 *> KL+KU superdiagonals.
76 *>
77 *> 3. If some U(i,i)=0, so that U is exactly singular, then the routine
78 *> returns with INFO = i. Otherwise, the factored form of A is used
79 *> to estimate the condition number of the matrix A. If the
80 *> reciprocal of the condition number is less than machine precision,
81 *> INFO = N+1 is returned as a warning, but the routine still goes on
82 *> to solve for X and compute error bounds as described below.
83 *>
84 *> 4. The system of equations is solved for X using the factored form
85 *> of A.
86 *>
87 *> 5. Iterative refinement is applied to improve the computed solution
88 *> matrix and calculate error bounds and backward error estimates
89 *> for it.
90 *>
91 *> 6. If equilibration was used, the matrix X is premultiplied by
92 *> diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
93 *> that it solves the original system before equilibration.
94 *> \endverbatim
95 *
96 * Arguments:
97 * ==========
98 *
99 *> \param[in] FACT
100 *> \verbatim
101 *> FACT is CHARACTER*1
102 *> Specifies whether or not the factored form of the matrix A is
103 *> supplied on entry, and if not, whether the matrix A should be
104 *> equilibrated before it is factored.
105 *> = 'F': On entry, AFB and IPIV contain the factored form of
106 *> A. If EQUED is not 'N', the matrix A has been
107 *> equilibrated with scaling factors given by R and C.
108 *> AB, AFB, and IPIV are not modified.
109 *> = 'N': The matrix A will be copied to AFB and factored.
110 *> = 'E': The matrix A will be equilibrated if necessary, then
111 *> copied to AFB and factored.
112 *> \endverbatim
113 *>
114 *> \param[in] TRANS
115 *> \verbatim
116 *> TRANS is CHARACTER*1
117 *> Specifies the form of the system of equations.
118 *> = 'N': A * X = B (No transpose)
119 *> = 'T': A**T * X = B (Transpose)
120 *> = 'C': A**H * X = B (Conjugate transpose)
121 *> \endverbatim
122 *>
123 *> \param[in] N
124 *> \verbatim
125 *> N is INTEGER
126 *> The number of linear equations, i.e., the order of the
127 *> matrix A. N >= 0.
128 *> \endverbatim
129 *>
130 *> \param[in] KL
131 *> \verbatim
132 *> KL is INTEGER
133 *> The number of subdiagonals within the band of A. KL >= 0.
134 *> \endverbatim
135 *>
136 *> \param[in] KU
137 *> \verbatim
138 *> KU is INTEGER
139 *> The number of superdiagonals within the band of A. KU >= 0.
140 *> \endverbatim
141 *>
142 *> \param[in] NRHS
143 *> \verbatim
144 *> NRHS is INTEGER
145 *> The number of right hand sides, i.e., the number of columns
146 *> of the matrices B and X. NRHS >= 0.
147 *> \endverbatim
148 *>
149 *> \param[in,out] AB
150 *> \verbatim
151 *> AB is COMPLEX*16 array, dimension (LDAB,N)
152 *> On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
153 *> The j-th column of A is stored in the j-th column of the
154 *> array AB as follows:
155 *> AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
156 *>
157 *> If FACT = 'F' and EQUED is not 'N', then A must have been
158 *> equilibrated by the scaling factors in R and/or C. AB is not
159 *> modified if FACT = 'F' or 'N', or if FACT = 'E' and
160 *> EQUED = 'N' on exit.
161 *>
162 *> On exit, if EQUED .ne. 'N', A is scaled as follows:
163 *> EQUED = 'R': A := diag(R) * A
164 *> EQUED = 'C': A := A * diag(C)
165 *> EQUED = 'B': A := diag(R) * A * diag(C).
166 *> \endverbatim
167 *>
168 *> \param[in] LDAB
169 *> \verbatim
170 *> LDAB is INTEGER
171 *> The leading dimension of the array AB. LDAB >= KL+KU+1.
172 *> \endverbatim
173 *>
174 *> \param[in,out] AFB
175 *> \verbatim
176 *> AFB is COMPLEX*16 array, dimension (LDAFB,N)
177 *> If FACT = 'F', then AFB is an input argument and on entry
178 *> contains details of the LU factorization of the band matrix
179 *> A, as computed by ZGBTRF. U is stored as an upper triangular
180 *> band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
181 *> and the multipliers used during the factorization are stored
182 *> in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is
183 *> the factored form of the equilibrated matrix A.
184 *>
185 *> If FACT = 'N', then AFB is an output argument and on exit
186 *> returns details of the LU factorization of A.
187 *>
188 *> If FACT = 'E', then AFB is an output argument and on exit
189 *> returns details of the LU factorization of the equilibrated
190 *> matrix A (see the description of AB for the form of the
191 *> equilibrated matrix).
192 *> \endverbatim
193 *>
194 *> \param[in] LDAFB
195 *> \verbatim
196 *> LDAFB is INTEGER
197 *> The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
198 *> \endverbatim
199 *>
200 *> \param[in,out] IPIV
201 *> \verbatim
202 *> IPIV is INTEGER array, dimension (N)
203 *> If FACT = 'F', then IPIV is an input argument and on entry
204 *> contains the pivot indices from the factorization A = L*U
205 *> as computed by ZGBTRF; row i of the matrix was interchanged
206 *> with row IPIV(i).
207 *>
208 *> If FACT = 'N', then IPIV is an output argument and on exit
209 *> contains the pivot indices from the factorization A = L*U
210 *> of the original matrix A.
211 *>
212 *> If FACT = 'E', then IPIV is an output argument and on exit
213 *> contains the pivot indices from the factorization A = L*U
214 *> of the equilibrated matrix A.
215 *> \endverbatim
216 *>
217 *> \param[in,out] EQUED
218 *> \verbatim
219 *> EQUED is CHARACTER*1
220 *> Specifies the form of equilibration that was done.
221 *> = 'N': No equilibration (always true if FACT = 'N').
222 *> = 'R': Row equilibration, i.e., A has been premultiplied by
223 *> diag(R).
224 *> = 'C': Column equilibration, i.e., A has been postmultiplied
225 *> by diag(C).
226 *> = 'B': Both row and column equilibration, i.e., A has been
227 *> replaced by diag(R) * A * diag(C).
228 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
229 *> output argument.
230 *> \endverbatim
231 *>
232 *> \param[in,out] R
233 *> \verbatim
234 *> R is DOUBLE PRECISION array, dimension (N)
235 *> The row scale factors for A. If EQUED = 'R' or 'B', A is
236 *> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
237 *> is not accessed. R is an input argument if FACT = 'F';
238 *> otherwise, R is an output argument. If FACT = 'F' and
239 *> EQUED = 'R' or 'B', each element of R must be positive.
240 *> \endverbatim
241 *>
242 *> \param[in,out] C
243 *> \verbatim
244 *> C is DOUBLE PRECISION array, dimension (N)
245 *> The column scale factors for A. If EQUED = 'C' or 'B', A is
246 *> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
247 *> is not accessed. C is an input argument if FACT = 'F';
248 *> otherwise, C is an output argument. If FACT = 'F' and
249 *> EQUED = 'C' or 'B', each element of C must be positive.
250 *> \endverbatim
251 *>
252 *> \param[in,out] B
253 *> \verbatim
254 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
255 *> On entry, the right hand side matrix B.
256 *> On exit,
257 *> if EQUED = 'N', B is not modified;
258 *> if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
259 *> diag(R)*B;
260 *> if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
261 *> overwritten by diag(C)*B.
262 *> \endverbatim
263 *>
264 *> \param[in] LDB
265 *> \verbatim
266 *> LDB is INTEGER
267 *> The leading dimension of the array B. LDB >= max(1,N).
268 *> \endverbatim
269 *>
270 *> \param[out] X
271 *> \verbatim
272 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
273 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
274 *> to the original system of equations. Note that A and B are
275 *> modified on exit if EQUED .ne. 'N', and the solution to the
276 *> equilibrated system is inv(diag(C))*X if TRANS = 'N' and
277 *> EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
278 *> and EQUED = 'R' or 'B'.
279 *> \endverbatim
280 *>
281 *> \param[in] LDX
282 *> \verbatim
283 *> LDX is INTEGER
284 *> The leading dimension of the array X. LDX >= max(1,N).
285 *> \endverbatim
286 *>
287 *> \param[out] RCOND
288 *> \verbatim
289 *> RCOND is DOUBLE PRECISION
290 *> The estimate of the reciprocal condition number of the matrix
291 *> A after equilibration (if done). If RCOND is less than the
292 *> machine precision (in particular, if RCOND = 0), the matrix
293 *> is singular to working precision. This condition is
294 *> indicated by a return code of INFO > 0.
295 *> \endverbatim
296 *>
297 *> \param[out] FERR
298 *> \verbatim
299 *> FERR is DOUBLE PRECISION array, dimension (NRHS)
300 *> The estimated forward error bound for each solution vector
301 *> X(j) (the j-th column of the solution matrix X).
302 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
303 *> is an estimated upper bound for the magnitude of the largest
304 *> element in (X(j) - XTRUE) divided by the magnitude of the
305 *> largest element in X(j). The estimate is as reliable as
306 *> the estimate for RCOND, and is almost always a slight
307 *> overestimate of the true error.
308 *> \endverbatim
309 *>
310 *> \param[out] BERR
311 *> \verbatim
312 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
313 *> The componentwise relative backward error of each solution
314 *> vector X(j) (i.e., the smallest relative change in
315 *> any element of A or B that makes X(j) an exact solution).
316 *> \endverbatim
317 *>
318 *> \param[out] WORK
319 *> \verbatim
320 *> WORK is COMPLEX*16 array, dimension (2*N)
321 *> \endverbatim
322 *>
323 *> \param[out] RWORK
324 *> \verbatim
325 *> RWORK is DOUBLE PRECISION array, dimension (N)
326 *> On exit, RWORK(1) contains the reciprocal pivot growth
327 *> factor norm(A)/norm(U). The "max absolute element" norm is
328 *> used. If RWORK(1) is much less than 1, then the stability
329 *> of the LU factorization of the (equilibrated) matrix A
330 *> could be poor. This also means that the solution X, condition
331 *> estimator RCOND, and forward error bound FERR could be
332 *> unreliable. If factorization fails with 0<INFO<=N, then
333 *> RWORK(1) contains the reciprocal pivot growth factor for the
334 *> leading INFO columns of A.
335 *> \endverbatim
336 *>
337 *> \param[out] INFO
338 *> \verbatim
339 *> INFO is INTEGER
340 *> = 0: successful exit
341 *> < 0: if INFO = -i, the i-th argument had an illegal value
342 *> > 0: if INFO = i, and i is
343 *> <= N: U(i,i) is exactly zero. The factorization
344 *> has been completed, but the factor U is exactly
345 *> singular, so the solution and error bounds
346 *> could not be computed. RCOND = 0 is returned.
347 *> = N+1: U is nonsingular, but RCOND is less than machine
348 *> precision, meaning that the matrix is singular
349 *> to working precision. Nevertheless, the
350 *> solution and error bounds are computed because
351 *> there are a number of situations where the
352 *> computed solution can be more accurate than the
353 *> value of RCOND would suggest.
354 *> \endverbatim
355 *
356 * Authors:
357 * ========
358 *
359 *> \author Univ. of Tennessee
360 *> \author Univ. of California Berkeley
361 *> \author Univ. of Colorado Denver
362 *> \author NAG Ltd.
363 *
364 *> \date April 2012
365 *
366 *> \ingroup complex16GBsolve
367 *
368 * =====================================================================
369  SUBROUTINE zgbsvx( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB,
370  $ ldafb, ipiv, equed, r, c, b, ldb, x, ldx,
371  $ rcond, ferr, berr, work, rwork, info )
372 *
373 * -- LAPACK driver routine (version 3.4.1) --
374 * -- LAPACK is a software package provided by Univ. of Tennessee, --
375 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
376 * April 2012
377 *
378 * .. Scalar Arguments ..
379  CHARACTER equed, fact, trans
380  INTEGER info, kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
381  DOUBLE PRECISION rcond
382 * ..
383 * .. Array Arguments ..
384  INTEGER ipiv( * )
385  DOUBLE PRECISION berr( * ), c( * ), ferr( * ), r( * ),
386  $ rwork( * )
387  COMPLEX*16 ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
388  $ work( * ), x( ldx, * )
389 * ..
390 *
391 * =====================================================================
392 * Moved setting of INFO = N+1 so INFO does not subsequently get
393 * overwritten. Sven, 17 Mar 05.
394 * =====================================================================
395 *
396 * .. Parameters ..
397  DOUBLE PRECISION zero, one
398  parameter( zero = 0.0d+0, one = 1.0d+0 )
399 * ..
400 * .. Local Scalars ..
401  LOGICAL colequ, equil, nofact, notran, rowequ
402  CHARACTER norm
403  INTEGER i, infequ, j, j1, j2
404  DOUBLE PRECISION amax, anorm, bignum, colcnd, rcmax, rcmin,
405  $ rowcnd, rpvgrw, smlnum
406 * ..
407 * .. External Functions ..
408  LOGICAL lsame
409  DOUBLE PRECISION dlamch, zlangb, zlantb
410  EXTERNAL lsame, dlamch, zlangb, zlantb
411 * ..
412 * .. External Subroutines ..
413  EXTERNAL xerbla, zcopy, zgbcon, zgbequ, zgbrfs, zgbtrf,
414  $ zgbtrs, zlacpy, zlaqgb
415 * ..
416 * .. Intrinsic Functions ..
417  INTRINSIC abs, max, min
418 * ..
419 * .. Executable Statements ..
420 *
421  info = 0
422  nofact = lsame( fact, 'N' )
423  equil = lsame( fact, 'E' )
424  notran = lsame( trans, 'N' )
425  IF( nofact .OR. equil ) THEN
426  equed = 'N'
427  rowequ = .false.
428  colequ = .false.
429  ELSE
430  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
431  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
432  smlnum = dlamch( 'Safe minimum' )
433  bignum = one / smlnum
434  END IF
435 *
436 * Test the input parameters.
437 *
438  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
439  $ THEN
440  info = -1
441  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
442  $ lsame( trans, 'C' ) ) THEN
443  info = -2
444  ELSE IF( n.LT.0 ) THEN
445  info = -3
446  ELSE IF( kl.LT.0 ) THEN
447  info = -4
448  ELSE IF( ku.LT.0 ) THEN
449  info = -5
450  ELSE IF( nrhs.LT.0 ) THEN
451  info = -6
452  ELSE IF( ldab.LT.kl+ku+1 ) THEN
453  info = -8
454  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
455  info = -10
456  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
457  $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
458  info = -12
459  ELSE
460  IF( rowequ ) THEN
461  rcmin = bignum
462  rcmax = zero
463  DO 10 j = 1, n
464  rcmin = min( rcmin, r( j ) )
465  rcmax = max( rcmax, r( j ) )
466  10 continue
467  IF( rcmin.LE.zero ) THEN
468  info = -13
469  ELSE IF( n.GT.0 ) THEN
470  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
471  ELSE
472  rowcnd = one
473  END IF
474  END IF
475  IF( colequ .AND. info.EQ.0 ) THEN
476  rcmin = bignum
477  rcmax = zero
478  DO 20 j = 1, n
479  rcmin = min( rcmin, c( j ) )
480  rcmax = max( rcmax, c( j ) )
481  20 continue
482  IF( rcmin.LE.zero ) THEN
483  info = -14
484  ELSE IF( n.GT.0 ) THEN
485  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
486  ELSE
487  colcnd = one
488  END IF
489  END IF
490  IF( info.EQ.0 ) THEN
491  IF( ldb.LT.max( 1, n ) ) THEN
492  info = -16
493  ELSE IF( ldx.LT.max( 1, n ) ) THEN
494  info = -18
495  END IF
496  END IF
497  END IF
498 *
499  IF( info.NE.0 ) THEN
500  CALL xerbla( 'ZGBSVX', -info )
501  return
502  END IF
503 *
504  IF( equil ) THEN
505 *
506 * Compute row and column scalings to equilibrate the matrix A.
507 *
508  CALL zgbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,
509  $ amax, infequ )
510  IF( infequ.EQ.0 ) THEN
511 *
512 * Equilibrate the matrix.
513 *
514  CALL zlaqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,
515  $ amax, equed )
516  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
517  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
518  END IF
519  END IF
520 *
521 * Scale the right hand side.
522 *
523  IF( notran ) THEN
524  IF( rowequ ) THEN
525  DO 40 j = 1, nrhs
526  DO 30 i = 1, n
527  b( i, j ) = r( i )*b( i, j )
528  30 continue
529  40 continue
530  END IF
531  ELSE IF( colequ ) THEN
532  DO 60 j = 1, nrhs
533  DO 50 i = 1, n
534  b( i, j ) = c( i )*b( i, j )
535  50 continue
536  60 continue
537  END IF
538 *
539  IF( nofact .OR. equil ) THEN
540 *
541 * Compute the LU factorization of the band matrix A.
542 *
543  DO 70 j = 1, n
544  j1 = max( j-ku, 1 )
545  j2 = min( j+kl, n )
546  CALL zcopy( j2-j1+1, ab( ku+1-j+j1, j ), 1,
547  $ afb( kl+ku+1-j+j1, j ), 1 )
548  70 continue
549 *
550  CALL zgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
551 *
552 * Return if INFO is non-zero.
553 *
554  IF( info.GT.0 ) THEN
555 *
556 * Compute the reciprocal pivot growth factor of the
557 * leading rank-deficient INFO columns of A.
558 *
559  anorm = zero
560  DO 90 j = 1, info
561  DO 80 i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
562  anorm = max( anorm, abs( ab( i, j ) ) )
563  80 continue
564  90 continue
565  rpvgrw = zlantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),
566  $ afb( max( 1, kl+ku+2-info ), 1 ), ldafb,
567  $ rwork )
568  IF( rpvgrw.EQ.zero ) THEN
569  rpvgrw = one
570  ELSE
571  rpvgrw = anorm / rpvgrw
572  END IF
573  rwork( 1 ) = rpvgrw
574  rcond = zero
575  return
576  END IF
577  END IF
578 *
579 * Compute the norm of the matrix A and the
580 * reciprocal pivot growth factor RPVGRW.
581 *
582  IF( notran ) THEN
583  norm = '1'
584  ELSE
585  norm = 'I'
586  END IF
587  anorm = zlangb( norm, n, kl, ku, ab, ldab, rwork )
588  rpvgrw = zlantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, rwork )
589  IF( rpvgrw.EQ.zero ) THEN
590  rpvgrw = one
591  ELSE
592  rpvgrw = zlangb( 'M', n, kl, ku, ab, ldab, rwork ) / rpvgrw
593  END IF
594 *
595 * Compute the reciprocal of the condition number of A.
596 *
597  CALL zgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,
598  $ work, rwork, info )
599 *
600 * Compute the solution matrix X.
601 *
602  CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
603  CALL zgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,
604  $ info )
605 *
606 * Use iterative refinement to improve the computed solution and
607 * compute error bounds and backward error estimates for it.
608 *
609  CALL zgbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,
610  $ b, ldb, x, ldx, ferr, berr, work, rwork, info )
611 *
612 * Transform the solution matrix X to a solution of the original
613 * system.
614 *
615  IF( notran ) THEN
616  IF( colequ ) THEN
617  DO 110 j = 1, nrhs
618  DO 100 i = 1, n
619  x( i, j ) = c( i )*x( i, j )
620  100 continue
621  110 continue
622  DO 120 j = 1, nrhs
623  ferr( j ) = ferr( j ) / colcnd
624  120 continue
625  END IF
626  ELSE IF( rowequ ) THEN
627  DO 140 j = 1, nrhs
628  DO 130 i = 1, n
629  x( i, j ) = r( i )*x( i, j )
630  130 continue
631  140 continue
632  DO 150 j = 1, nrhs
633  ferr( j ) = ferr( j ) / rowcnd
634  150 continue
635  END IF
636 *
637 * Set INFO = N+1 if the matrix is singular to working precision.
638 *
639  IF( rcond.LT.dlamch( 'Epsilon' ) )
640  $ info = n + 1
641 *
642  rwork( 1 ) = rpvgrw
643  return
644 *
645 * End of ZGBSVX
646 *
647  END