LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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
subroutine zgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
ZGBTRF
Definition: zgbtrf.f:146
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
ZGBTRS
Definition: zgbtrs.f:140
subroutine zgbequ(M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)
ZGBEQU
Definition: zgbequ.f:156
subroutine zgbrfs(TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZGBRFS
Definition: zgbrfs.f:208
subroutine zgbcon(NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, WORK, RWORK, INFO)
ZGBCON
Definition: zgbcon.f:149
subroutine zlaqgb(M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, EQUED)
ZLAQGB scales a general band matrix, using row and column scaling factors computed by sgbequ...
Definition: zlaqgb.f:162
subroutine zgbsvx(FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO)
ZGBSVX computes the solution to system of linear equations A * X = B for GB matrices ...
Definition: zgbsvx.f:372