LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dgbsvxx.f
Go to the documentation of this file.
1 *> \brief <b> DGBSVXX 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 DGBSVXX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbsvxx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbsvxx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbsvxx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGBSVXX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB,
22 * LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX,
23 * RCOND, RPVGRW, BERR, N_ERR_BNDS,
24 * ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
25 * WORK, IWORK, INFO )
26 *
27 * .. Scalar Arguments ..
28 * CHARACTER EQUED, FACT, TRANS
29 * INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, NRHS, NPARAMS,
30 * $ N_ERR_BNDS, KL, KU
31 * DOUBLE PRECISION RCOND, RPVGRW
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IPIV( * ), IWORK( * )
35 * DOUBLE PRECISION AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
36 * $ X( LDX , * ),WORK( * )
37 * DOUBLE PRECISION R( * ), C( * ), PARAMS( * ), BERR( * ),
38 * $ ERR_BNDS_NORM( NRHS, * ),
39 * $ ERR_BNDS_COMP( NRHS, * )
40 * ..
41 *
42 *
43 *> \par Purpose:
44 * =============
45 *>
46 *> \verbatim
47 *>
48 *> DGBSVXX uses the LU factorization to compute the solution to a
49 *> double precision system of linear equations A * X = B, where A is an
50 *> N-by-N matrix and X and B are N-by-NRHS matrices.
51 *>
52 *> If requested, both normwise and maximum componentwise error bounds
53 *> are returned. DGBSVXX will return a solution with a tiny
54 *> guaranteed error (O(eps) where eps is the working machine
55 *> precision) unless the matrix is very ill-conditioned, in which
56 *> case a warning is returned. Relevant condition numbers also are
57 *> calculated and returned.
58 *>
59 *> DGBSVXX accepts user-provided factorizations and equilibration
60 *> factors; see the definitions of the FACT and EQUED options.
61 *> Solving with refinement and using a factorization from a previous
62 *> DGBSVXX call will also produce a solution with either O(eps)
63 *> errors or warnings, but we cannot make that claim for general
64 *> user-provided factorizations and equilibration factors if they
65 *> differ from what DGBSVXX would itself produce.
66 *> \endverbatim
67 *
68 *> \par Description:
69 * =================
70 *>
71 *> \verbatim
72 *>
73 *> The following steps are performed:
74 *>
75 *> 1. If FACT = 'E', double precision scaling factors are computed to equilibrate
76 *> the system:
77 *>
78 *> TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B
79 *> TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B
80 *> TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B
81 *>
82 *> Whether or not the system will be equilibrated depends on the
83 *> scaling of the matrix A, but if equilibration is used, A is
84 *> overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N')
85 *> or diag(C)*B (if TRANS = 'T' or 'C').
86 *>
87 *> 2. If FACT = 'N' or 'E', the LU decomposition is used to factor
88 *> the matrix A (after equilibration if FACT = 'E') as
89 *>
90 *> A = P * L * U,
91 *>
92 *> where P is a permutation matrix, L is a unit lower triangular
93 *> matrix, and U is upper triangular.
94 *>
95 *> 3. If some U(i,i)=0, so that U is exactly singular, then the
96 *> routine returns with INFO = i. Otherwise, the factored form of A
97 *> is used to estimate the condition number of the matrix A (see
98 *> argument RCOND). If the reciprocal of the condition number is less
99 *> than machine precision, the routine still goes on to solve for X
100 *> and compute error bounds as described below.
101 *>
102 *> 4. The system of equations is solved for X using the factored form
103 *> of A.
104 *>
105 *> 5. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
106 *> the routine will use iterative refinement to try to get a small
107 *> error and error bounds. Refinement calculates the residual to at
108 *> least twice the working precision.
109 *>
110 *> 6. If equilibration was used, the matrix X is premultiplied by
111 *> diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so
112 *> that it solves the original system before equilibration.
113 *> \endverbatim
114 *
115 * Arguments:
116 * ==========
117 *
118 *> \verbatim
119 *> Some optional parameters are bundled in the PARAMS array. These
120 *> settings determine how refinement is performed, but often the
121 *> defaults are acceptable. If the defaults are acceptable, users
122 *> can pass NPARAMS = 0 which prevents the source code from accessing
123 *> the PARAMS argument.
124 *> \endverbatim
125 *>
126 *> \param[in] FACT
127 *> \verbatim
128 *> FACT is CHARACTER*1
129 *> Specifies whether or not the factored form of the matrix A is
130 *> supplied on entry, and if not, whether the matrix A should be
131 *> equilibrated before it is factored.
132 *> = 'F': On entry, AF and IPIV contain the factored form of A.
133 *> If EQUED is not 'N', the matrix A has been
134 *> equilibrated with scaling factors given by R and C.
135 *> A, AF, and IPIV are not modified.
136 *> = 'N': The matrix A will be copied to AF and factored.
137 *> = 'E': The matrix A will be equilibrated if necessary, then
138 *> copied to AF and factored.
139 *> \endverbatim
140 *>
141 *> \param[in] TRANS
142 *> \verbatim
143 *> TRANS is CHARACTER*1
144 *> Specifies the form of the system of equations:
145 *> = 'N': A * X = B (No transpose)
146 *> = 'T': A**T * X = B (Transpose)
147 *> = 'C': A**H * X = B (Conjugate Transpose = Transpose)
148 *> \endverbatim
149 *>
150 *> \param[in] N
151 *> \verbatim
152 *> N is INTEGER
153 *> The number of linear equations, i.e., the order of the
154 *> matrix A. N >= 0.
155 *> \endverbatim
156 *>
157 *> \param[in] KL
158 *> \verbatim
159 *> KL is INTEGER
160 *> The number of subdiagonals within the band of A. KL >= 0.
161 *> \endverbatim
162 *>
163 *> \param[in] KU
164 *> \verbatim
165 *> KU is INTEGER
166 *> The number of superdiagonals within the band of A. KU >= 0.
167 *> \endverbatim
168 *>
169 *> \param[in] NRHS
170 *> \verbatim
171 *> NRHS is INTEGER
172 *> The number of right hand sides, i.e., the number of columns
173 *> of the matrices B and X. NRHS >= 0.
174 *> \endverbatim
175 *>
176 *> \param[in,out] AB
177 *> \verbatim
178 *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
179 *> On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
180 *> The j-th column of A is stored in the j-th column of the
181 *> array AB as follows:
182 *> AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
183 *>
184 *> If FACT = 'F' and EQUED is not 'N', then AB must have been
185 *> equilibrated by the scaling factors in R and/or C. AB is not
186 *> modified if FACT = 'F' or 'N', or if FACT = 'E' and
187 *> EQUED = 'N' on exit.
188 *>
189 *> On exit, if EQUED .ne. 'N', A is scaled as follows:
190 *> EQUED = 'R': A := diag(R) * A
191 *> EQUED = 'C': A := A * diag(C)
192 *> EQUED = 'B': A := diag(R) * A * diag(C).
193 *> \endverbatim
194 *>
195 *> \param[in] LDAB
196 *> \verbatim
197 *> LDAB is INTEGER
198 *> The leading dimension of the array AB. LDAB >= KL+KU+1.
199 *> \endverbatim
200 *>
201 *> \param[in,out] AFB
202 *> \verbatim
203 *> AFB is DOUBLE PRECISION array, dimension (LDAFB,N)
204 *> If FACT = 'F', then AFB is an input argument and on entry
205 *> contains details of the LU factorization of the band matrix
206 *> A, as computed by DGBTRF. U is stored as an upper triangular
207 *> band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1,
208 *> and the multipliers used during the factorization are stored
209 *> in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is
210 *> the factored form of the equilibrated matrix A.
211 *>
212 *> If FACT = 'N', then AF is an output argument and on exit
213 *> returns the factors L and U from the factorization A = P*L*U
214 *> of the original matrix A.
215 *>
216 *> If FACT = 'E', then AF is an output argument and on exit
217 *> returns the factors L and U from the factorization A = P*L*U
218 *> of the equilibrated matrix A (see the description of A for
219 *> the form of the equilibrated matrix).
220 *> \endverbatim
221 *>
222 *> \param[in] LDAFB
223 *> \verbatim
224 *> LDAFB is INTEGER
225 *> The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1.
226 *> \endverbatim
227 *>
228 *> \param[in,out] IPIV
229 *> \verbatim
230 *> IPIV is INTEGER array, dimension (N)
231 *> If FACT = 'F', then IPIV is an input argument and on entry
232 *> contains the pivot indices from the factorization A = P*L*U
233 *> as computed by DGETRF; row i of the matrix was interchanged
234 *> with row IPIV(i).
235 *>
236 *> If FACT = 'N', then IPIV is an output argument and on exit
237 *> contains the pivot indices from the factorization A = P*L*U
238 *> of the original matrix A.
239 *>
240 *> If FACT = 'E', then IPIV is an output argument and on exit
241 *> contains the pivot indices from the factorization A = P*L*U
242 *> of the equilibrated matrix A.
243 *> \endverbatim
244 *>
245 *> \param[in,out] EQUED
246 *> \verbatim
247 *> EQUED is CHARACTER*1
248 *> Specifies the form of equilibration that was done.
249 *> = 'N': No equilibration (always true if FACT = 'N').
250 *> = 'R': Row equilibration, i.e., A has been premultiplied by
251 *> diag(R).
252 *> = 'C': Column equilibration, i.e., A has been postmultiplied
253 *> by diag(C).
254 *> = 'B': Both row and column equilibration, i.e., A has been
255 *> replaced by diag(R) * A * diag(C).
256 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
257 *> output argument.
258 *> \endverbatim
259 *>
260 *> \param[in,out] R
261 *> \verbatim
262 *> R is DOUBLE PRECISION array, dimension (N)
263 *> The row scale factors for A. If EQUED = 'R' or 'B', A is
264 *> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
265 *> is not accessed. R is an input argument if FACT = 'F';
266 *> otherwise, R is an output argument. If FACT = 'F' and
267 *> EQUED = 'R' or 'B', each element of R must be positive.
268 *> If R is output, each element of R is a power of the radix.
269 *> If R is input, each element of R should be a power of the radix
270 *> to ensure a reliable solution and error estimates. Scaling by
271 *> powers of the radix does not cause rounding errors unless the
272 *> result underflows or overflows. Rounding errors during scaling
273 *> lead to refining with a matrix that is not equivalent to the
274 *> input matrix, producing error estimates that may not be
275 *> reliable.
276 *> \endverbatim
277 *>
278 *> \param[in,out] C
279 *> \verbatim
280 *> C is DOUBLE PRECISION array, dimension (N)
281 *> The column scale factors for A. If EQUED = 'C' or 'B', A is
282 *> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
283 *> is not accessed. C is an input argument if FACT = 'F';
284 *> otherwise, C is an output argument. If FACT = 'F' and
285 *> EQUED = 'C' or 'B', each element of C must be positive.
286 *> If C is output, each element of C is a power of the radix.
287 *> If C is input, each element of C should be a power of the radix
288 *> to ensure a reliable solution and error estimates. Scaling by
289 *> powers of the radix does not cause rounding errors unless the
290 *> result underflows or overflows. Rounding errors during scaling
291 *> lead to refining with a matrix that is not equivalent to the
292 *> input matrix, producing error estimates that may not be
293 *> reliable.
294 *> \endverbatim
295 *>
296 *> \param[in,out] B
297 *> \verbatim
298 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
299 *> On entry, the N-by-NRHS right hand side matrix B.
300 *> On exit,
301 *> if EQUED = 'N', B is not modified;
302 *> if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
303 *> diag(R)*B;
304 *> if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
305 *> overwritten by diag(C)*B.
306 *> \endverbatim
307 *>
308 *> \param[in] LDB
309 *> \verbatim
310 *> LDB is INTEGER
311 *> The leading dimension of the array B. LDB >= max(1,N).
312 *> \endverbatim
313 *>
314 *> \param[out] X
315 *> \verbatim
316 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
317 *> If INFO = 0, the N-by-NRHS solution matrix X to the original
318 *> system of equations. Note that A and B are modified on exit
319 *> if EQUED .ne. 'N', and the solution to the equilibrated system is
320 *> inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
321 *> inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
322 *> \endverbatim
323 *>
324 *> \param[in] LDX
325 *> \verbatim
326 *> LDX is INTEGER
327 *> The leading dimension of the array X. LDX >= max(1,N).
328 *> \endverbatim
329 *>
330 *> \param[out] RCOND
331 *> \verbatim
332 *> RCOND is DOUBLE PRECISION
333 *> Reciprocal scaled condition number. This is an estimate of the
334 *> reciprocal Skeel condition number of the matrix A after
335 *> equilibration (if done). If this is less than the machine
336 *> precision (in particular, if it is zero), the matrix is singular
337 *> to working precision. Note that the error may still be small even
338 *> if this number is very small and the matrix appears ill-
339 *> conditioned.
340 *> \endverbatim
341 *>
342 *> \param[out] RPVGRW
343 *> \verbatim
344 *> RPVGRW is DOUBLE PRECISION
345 *> Reciprocal pivot growth. On exit, this contains the reciprocal
346 *> pivot growth factor norm(A)/norm(U). The "max absolute element"
347 *> norm is used. If this is much less than 1, then the stability of
348 *> the LU factorization of the (equilibrated) matrix A could be poor.
349 *> This also means that the solution X, estimated condition numbers,
350 *> and error bounds could be unreliable. If factorization fails with
351 *> 0<INFO<=N, then this contains the reciprocal pivot growth factor
352 *> for the leading INFO columns of A. In DGESVX, this quantity is
353 *> returned in WORK(1).
354 *> \endverbatim
355 *>
356 *> \param[out] BERR
357 *> \verbatim
358 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
359 *> Componentwise relative backward error. This is the
360 *> componentwise relative backward error of each solution vector X(j)
361 *> (i.e., the smallest relative change in any element of A or B that
362 *> makes X(j) an exact solution).
363 *> \endverbatim
364 *>
365 *> \param[in] N_ERR_BNDS
366 *> \verbatim
367 *> N_ERR_BNDS is INTEGER
368 *> Number of error bounds to return for each right hand side
369 *> and each type (normwise or componentwise). See ERR_BNDS_NORM and
370 *> ERR_BNDS_COMP below.
371 *> \endverbatim
372 *>
373 *> \param[out] ERR_BNDS_NORM
374 *> \verbatim
375 *> ERR_BNDS_NORM is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
376 *> For each right-hand side, this array contains information about
377 *> various error bounds and condition numbers corresponding to the
378 *> normwise relative error, which is defined as follows:
379 *>
380 *> Normwise relative error in the ith solution vector:
381 *> max_j (abs(XTRUE(j,i) - X(j,i)))
382 *> ------------------------------
383 *> max_j abs(X(j,i))
384 *>
385 *> The array is indexed by the type of error information as described
386 *> below. There currently are up to three pieces of information
387 *> returned.
388 *>
389 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
390 *> right-hand side.
391 *>
392 *> The second index in ERR_BNDS_NORM(:,err) contains the following
393 *> three fields:
394 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
395 *> reciprocal condition number is less than the threshold
396 *> sqrt(n) * dlamch('Epsilon').
397 *>
398 *> err = 2 "Guaranteed" error bound: The estimated forward error,
399 *> almost certainly within a factor of 10 of the true error
400 *> so long as the next entry is greater than the threshold
401 *> sqrt(n) * dlamch('Epsilon'). This error bound should only
402 *> be trusted if the previous boolean is true.
403 *>
404 *> err = 3 Reciprocal condition number: Estimated normwise
405 *> reciprocal condition number. Compared with the threshold
406 *> sqrt(n) * dlamch('Epsilon') to determine if the error
407 *> estimate is "guaranteed". These reciprocal condition
408 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
409 *> appropriately scaled matrix Z.
410 *> Let Z = S*A, where S scales each row by a power of the
411 *> radix so all absolute row sums of Z are approximately 1.
412 *>
413 *> See Lapack Working Note 165 for further details and extra
414 *> cautions.
415 *> \endverbatim
416 *>
417 *> \param[out] ERR_BNDS_COMP
418 *> \verbatim
419 *> ERR_BNDS_COMP is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
420 *> For each right-hand side, this array contains information about
421 *> various error bounds and condition numbers corresponding to the
422 *> componentwise relative error, which is defined as follows:
423 *>
424 *> Componentwise relative error in the ith solution vector:
425 *> abs(XTRUE(j,i) - X(j,i))
426 *> max_j ----------------------
427 *> abs(X(j,i))
428 *>
429 *> The array is indexed by the right-hand side i (on which the
430 *> componentwise relative error depends), and the type of error
431 *> information as described below. There currently are up to three
432 *> pieces of information returned for each right-hand side. If
433 *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
434 *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
435 *> the first (:,N_ERR_BNDS) entries are returned.
436 *>
437 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
438 *> right-hand side.
439 *>
440 *> The second index in ERR_BNDS_COMP(:,err) contains the following
441 *> three fields:
442 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
443 *> reciprocal condition number is less than the threshold
444 *> sqrt(n) * dlamch('Epsilon').
445 *>
446 *> err = 2 "Guaranteed" error bound: The estimated forward error,
447 *> almost certainly within a factor of 10 of the true error
448 *> so long as the next entry is greater than the threshold
449 *> sqrt(n) * dlamch('Epsilon'). This error bound should only
450 *> be trusted if the previous boolean is true.
451 *>
452 *> err = 3 Reciprocal condition number: Estimated componentwise
453 *> reciprocal condition number. Compared with the threshold
454 *> sqrt(n) * dlamch('Epsilon') to determine if the error
455 *> estimate is "guaranteed". These reciprocal condition
456 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
457 *> appropriately scaled matrix Z.
458 *> Let Z = S*(A*diag(x)), where x is the solution for the
459 *> current right-hand side and S scales each row of
460 *> A*diag(x) by a power of the radix so all absolute row
461 *> sums of Z are approximately 1.
462 *>
463 *> See Lapack Working Note 165 for further details and extra
464 *> cautions.
465 *> \endverbatim
466 *>
467 *> \param[in] NPARAMS
468 *> \verbatim
469 *> NPARAMS is INTEGER
470 *> Specifies the number of parameters set in PARAMS. If .LE. 0, the
471 *> PARAMS array is never referenced and default values are used.
472 *> \endverbatim
473 *>
474 *> \param[in,out] PARAMS
475 *> \verbatim
476 *> PARAMS is DOUBLE PRECISION array, dimension (NPARAMS)
477 *> Specifies algorithm parameters. If an entry is .LT. 0.0, then
478 *> that entry will be filled with default value used for that
479 *> parameter. Only positions up to NPARAMS are accessed; defaults
480 *> are used for higher-numbered parameters.
481 *>
482 *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
483 *> refinement or not.
484 *> Default: 1.0D+0
485 *> = 0.0 : No refinement is performed, and no error bounds are
486 *> computed.
487 *> = 1.0 : Use the extra-precise refinement algorithm.
488 *> (other values are reserved for future use)
489 *>
490 *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
491 *> computations allowed for refinement.
492 *> Default: 10
493 *> Aggressive: Set to 100 to permit convergence using approximate
494 *> factorizations or factorizations other than LU. If
495 *> the factorization uses a technique other than
496 *> Gaussian elimination, the guarantees in
497 *> err_bnds_norm and err_bnds_comp may no longer be
498 *> trustworthy.
499 *>
500 *> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
501 *> will attempt to find a solution with small componentwise
502 *> relative error in the double-precision algorithm. Positive
503 *> is true, 0.0 is false.
504 *> Default: 1.0 (attempt componentwise convergence)
505 *> \endverbatim
506 *>
507 *> \param[out] WORK
508 *> \verbatim
509 *> WORK is DOUBLE PRECISION array, dimension (4*N)
510 *> \endverbatim
511 *>
512 *> \param[out] IWORK
513 *> \verbatim
514 *> IWORK is INTEGER array, dimension (N)
515 *> \endverbatim
516 *>
517 *> \param[out] INFO
518 *> \verbatim
519 *> INFO is INTEGER
520 *> = 0: Successful exit. The solution to every right-hand side is
521 *> guaranteed.
522 *> < 0: If INFO = -i, the i-th argument had an illegal value
523 *> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
524 *> has been completed, but the factor U is exactly singular, so
525 *> the solution and error bounds could not be computed. RCOND = 0
526 *> is returned.
527 *> = N+J: The solution corresponding to the Jth right-hand side is
528 *> not guaranteed. The solutions corresponding to other right-
529 *> hand sides K with K > J may not be guaranteed as well, but
530 *> only the first such right-hand side is reported. If a small
531 *> componentwise error is not requested (PARAMS(3) = 0.0) then
532 *> the Jth right-hand side is the first with a normwise error
533 *> bound that is not guaranteed (the smallest J such
534 *> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
535 *> the Jth right-hand side is the first with either a normwise or
536 *> componentwise error bound that is not guaranteed (the smallest
537 *> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
538 *> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
539 *> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
540 *> about all of the right-hand sides check ERR_BNDS_NORM or
541 *> ERR_BNDS_COMP.
542 *> \endverbatim
543 *
544 * Authors:
545 * ========
546 *
547 *> \author Univ. of Tennessee
548 *> \author Univ. of California Berkeley
549 *> \author Univ. of Colorado Denver
550 *> \author NAG Ltd.
551 *
552 *> \date April 2012
553 *
554 *> \ingroup doubleGBsolve
555 *
556 * =====================================================================
557  SUBROUTINE dgbsvxx( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB,
558  $ ldafb, ipiv, equed, r, c, b, ldb, x, ldx,
559  $ rcond, rpvgrw, berr, n_err_bnds,
560  $ err_bnds_norm, err_bnds_comp, nparams, params,
561  $ work, iwork, info )
562 *
563 * -- LAPACK driver routine (version 3.4.1) --
564 * -- LAPACK is a software package provided by Univ. of Tennessee, --
565 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
566 * April 2012
567 *
568 * .. Scalar Arguments ..
569  CHARACTER equed, fact, trans
570  INTEGER info, ldab, ldafb, ldb, ldx, n, nrhs, nparams,
571  $ n_err_bnds, kl, ku
572  DOUBLE PRECISION rcond, rpvgrw
573 * ..
574 * .. Array Arguments ..
575  INTEGER ipiv( * ), iwork( * )
576  DOUBLE PRECISION ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
577  $ x( ldx , * ),work( * )
578  DOUBLE PRECISION r( * ), c( * ), params( * ), berr( * ),
579  $ err_bnds_norm( nrhs, * ),
580  $ err_bnds_comp( nrhs, * )
581 * ..
582 *
583 * ==================================================================
584 *
585 * .. Parameters ..
586  DOUBLE PRECISION zero, one
587  parameter( zero = 0.0d+0, one = 1.0d+0 )
588  INTEGER final_nrm_err_i, final_cmp_err_i, berr_i
589  INTEGER rcond_i, nrm_rcond_i, nrm_err_i, cmp_rcond_i
590  INTEGER cmp_err_i, piv_growth_i
591  parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
592  $ berr_i = 3 )
593  parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
594  parameter( cmp_rcond_i = 7, cmp_err_i = 8,
595  $ piv_growth_i = 9 )
596 * ..
597 * .. Local Scalars ..
598  LOGICAL colequ, equil, nofact, notran, rowequ
599  INTEGER infequ, i, j
600  DOUBLE PRECISION amax, bignum, colcnd, rcmax, rcmin,
601  $ rowcnd, smlnum
602 * ..
603 * .. External Functions ..
604  EXTERNAL lsame, dlamch, dla_gbrpvgrw
605  LOGICAL lsame
606  DOUBLE PRECISION dlamch, dla_gbrpvgrw
607 * ..
608 * .. External Subroutines ..
609  EXTERNAL dgbequb, dgbtrf, dgbtrs, dlacpy, dlaqgb,
611 * ..
612 * .. Intrinsic Functions ..
613  INTRINSIC max, min
614 * ..
615 * .. Executable Statements ..
616 *
617  info = 0
618  nofact = lsame( fact, 'N' )
619  equil = lsame( fact, 'E' )
620  notran = lsame( trans, 'N' )
621  smlnum = dlamch( 'Safe minimum' )
622  bignum = one / smlnum
623  IF( nofact .OR. equil ) THEN
624  equed = 'N'
625  rowequ = .false.
626  colequ = .false.
627  ELSE
628  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
629  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
630  END IF
631 *
632 * Default is failure. If an input parameter is wrong or
633 * factorization fails, make everything look horrible. Only the
634 * pivot growth is set here, the rest is initialized in DGBRFSX.
635 *
636  rpvgrw = zero
637 *
638 * Test the input parameters. PARAMS is not tested until DGBRFSX.
639 *
640  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
641  $ lsame( fact, 'F' ) ) THEN
642  info = -1
643  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
644  $ lsame( trans, 'C' ) ) THEN
645  info = -2
646  ELSE IF( n.LT.0 ) THEN
647  info = -3
648  ELSE IF( kl.LT.0 ) THEN
649  info = -4
650  ELSE IF( ku.LT.0 ) THEN
651  info = -5
652  ELSE IF( nrhs.LT.0 ) THEN
653  info = -6
654  ELSE IF( ldab.LT.kl+ku+1 ) THEN
655  info = -8
656  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
657  info = -10
658  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
659  $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
660  info = -12
661  ELSE
662  IF( rowequ ) THEN
663  rcmin = bignum
664  rcmax = zero
665  DO 10 j = 1, n
666  rcmin = min( rcmin, r( j ) )
667  rcmax = max( rcmax, r( j ) )
668  10 CONTINUE
669  IF( rcmin.LE.zero ) THEN
670  info = -13
671  ELSE IF( n.GT.0 ) THEN
672  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
673  ELSE
674  rowcnd = one
675  END IF
676  END IF
677  IF( colequ .AND. info.EQ.0 ) THEN
678  rcmin = bignum
679  rcmax = zero
680  DO 20 j = 1, n
681  rcmin = min( rcmin, c( j ) )
682  rcmax = max( rcmax, c( j ) )
683  20 CONTINUE
684  IF( rcmin.LE.zero ) THEN
685  info = -14
686  ELSE IF( n.GT.0 ) THEN
687  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
688  ELSE
689  colcnd = one
690  END IF
691  END IF
692  IF( info.EQ.0 ) THEN
693  IF( ldb.LT.max( 1, n ) ) THEN
694  info = -15
695  ELSE IF( ldx.LT.max( 1, n ) ) THEN
696  info = -16
697  END IF
698  END IF
699  END IF
700 *
701  IF( info.NE.0 ) THEN
702  CALL xerbla( 'DGBSVXX', -info )
703  RETURN
704  END IF
705 *
706  IF( equil ) THEN
707 *
708 * Compute row and column scalings to equilibrate the matrix A.
709 *
710  CALL dgbequb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,
711  $ amax, infequ )
712  IF( infequ.EQ.0 ) THEN
713 *
714 * Equilibrate the matrix.
715 *
716  CALL dlaqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,
717  $ amax, equed )
718  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
719  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
720  END IF
721 *
722 * If the scaling factors are not applied, set them to 1.0.
723 *
724  IF ( .NOT.rowequ ) THEN
725  DO j = 1, n
726  r( j ) = 1.0d+0
727  END DO
728  END IF
729  IF ( .NOT.colequ ) THEN
730  DO j = 1, n
731  c( j ) = 1.0d+0
732  END DO
733  END IF
734  END IF
735 *
736 * Scale the right hand side.
737 *
738  IF( notran ) THEN
739  IF( rowequ ) CALL dlascl2(n, nrhs, r, b, ldb)
740  ELSE
741  IF( colequ ) CALL dlascl2(n, nrhs, c, b, ldb)
742  END IF
743 *
744  IF( nofact .OR. equil ) THEN
745 *
746 * Compute the LU factorization of A.
747 *
748  DO 40, j = 1, n
749  DO 30, i = kl+1, 2*kl+ku+1
750  afb( i, j ) = ab( i-kl, j )
751  30 CONTINUE
752  40 CONTINUE
753  CALL dgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
754 *
755 * Return if INFO is non-zero.
756 *
757  IF( info.GT.0 ) THEN
758 *
759 * Pivot in column INFO is exactly 0
760 * Compute the reciprocal pivot growth factor of the
761 * leading rank-deficient INFO columns of A.
762 *
763  rpvgrw = dla_gbrpvgrw( n, kl, ku, info, ab, ldab, afb,
764  $ ldafb )
765  RETURN
766  END IF
767  END IF
768 *
769 * Compute the reciprocal pivot growth factor RPVGRW.
770 *
771  rpvgrw = dla_gbrpvgrw( n, kl, ku, n, ab, ldab, afb, ldafb )
772 *
773 * Compute the solution matrix X.
774 *
775  CALL dlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
776  CALL dgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,
777  $ info )
778 *
779 * Use iterative refinement to improve the computed solution and
780 * compute error bounds and backward error estimates for it.
781 *
782  CALL dgbrfsx( trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb,
783  $ ipiv, r, c, b, ldb, x, ldx, rcond, berr,
784  $ n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params,
785  $ work, iwork, info )
786 *
787 * Scale solutions.
788 *
789  IF ( colequ .AND. notran ) THEN
790  CALL dlascl2( n, nrhs, c, x, ldx )
791  ELSE IF ( rowequ .AND. .NOT.notran ) THEN
792  CALL dlascl2( n, nrhs, r, x, ldx )
793  END IF
794 *
795  RETURN
796 *
797 * End of DGBSVXX
798 *
799  END