LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
cgbsvxx.f
Go to the documentation of this file.
1 *> \brief <b> CGBSVXX 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbsvxx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbsvxx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbsvxx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGBSVXX( 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, RWORK, INFO )
26 *
27 * .. Scalar Arguments ..
28 * CHARACTER EQUED, FACT, TRANS
29 * INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, NRHS, NPARAMS,
30 * \$ N_ERR_BNDS
31 * REAL RCOND, RPVGRW
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IPIV( * )
35 * COMPLEX AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
36 * \$ X( LDX , * ),WORK( * )
37 * REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
38 * \$ ERR_BNDS_NORM( NRHS, * ),
39 * \$ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
40 * ..
41 *
42 *
43 *> \par Purpose:
44 * =============
45 *>
46 *> \verbatim
47 *>
48 *> CGBSVXX uses the LU factorization to compute the solution to a
49 *> complex 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. CGBSVXX 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 *> CGBSVXX 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 *> CGBSVXX 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 CGBSVXX 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', real 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 COMPLEX 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 COMPLEX 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 CGBTRF. 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 SGETRF; 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 REAL 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 REAL 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 COMPLEX 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 COMPLEX 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 REAL
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 REAL
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 SGESVX, this quantity is
353 *> returned in WORK(1).
354 *> \endverbatim
355 *>
356 *> \param[out] BERR
357 *> \verbatim
358 *> BERR is REAL 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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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) * slamch('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) * slamch('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) * slamch('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 REAL 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.0
485 *> = 0.0 : No refinement is performed, and no error bounds are
486 *> computed.
487 *> = 1.0 : Use the double-precision refinement algorithm,
488 *> possibly with doubled-single computations if the
489 *> compilation environment does not support DOUBLE
490 *> PRECISION.
491 *> (other values are reserved for future use)
492 *>
493 *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
494 *> computations allowed for refinement.
495 *> Default: 10
496 *> Aggressive: Set to 100 to permit convergence using approximate
497 *> factorizations or factorizations other than LU. If
498 *> the factorization uses a technique other than
499 *> Gaussian elimination, the guarantees in
500 *> err_bnds_norm and err_bnds_comp may no longer be
501 *> trustworthy.
502 *>
503 *> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
504 *> will attempt to find a solution with small componentwise
505 *> relative error in the double-precision algorithm. Positive
506 *> is true, 0.0 is false.
507 *> Default: 1.0 (attempt componentwise convergence)
508 *> \endverbatim
509 *>
510 *> \param[out] WORK
511 *> \verbatim
512 *> WORK is COMPLEX array, dimension (2*N)
513 *> \endverbatim
514 *>
515 *> \param[out] RWORK
516 *> \verbatim
517 *> RWORK is REAL array, dimension (2*N)
518 *> \endverbatim
519 *>
520 *> \param[out] INFO
521 *> \verbatim
522 *> INFO is INTEGER
523 *> = 0: Successful exit. The solution to every right-hand side is
524 *> guaranteed.
525 *> < 0: If INFO = -i, the i-th argument had an illegal value
526 *> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
527 *> has been completed, but the factor U is exactly singular, so
528 *> the solution and error bounds could not be computed. RCOND = 0
529 *> is returned.
530 *> = N+J: The solution corresponding to the Jth right-hand side is
531 *> not guaranteed. The solutions corresponding to other right-
532 *> hand sides K with K > J may not be guaranteed as well, but
533 *> only the first such right-hand side is reported. If a small
534 *> componentwise error is not requested (PARAMS(3) = 0.0) then
535 *> the Jth right-hand side is the first with a normwise error
536 *> bound that is not guaranteed (the smallest J such
537 *> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
538 *> the Jth right-hand side is the first with either a normwise or
539 *> componentwise error bound that is not guaranteed (the smallest
540 *> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
541 *> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
542 *> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
543 *> about all of the right-hand sides check ERR_BNDS_NORM or
544 *> ERR_BNDS_COMP.
545 *> \endverbatim
546 *
547 * Authors:
548 * ========
549 *
550 *> \author Univ. of Tennessee
551 *> \author Univ. of California Berkeley
552 *> \author Univ. of Colorado Denver
553 *> \author NAG Ltd.
554 *
555 *> \date April 2012
556 *
557 *> \ingroup complexGBsolve
558 *
559 * =====================================================================
560  SUBROUTINE cgbsvxx( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB,
561  \$ ldafb, ipiv, equed, r, c, b, ldb, x, ldx,
562  \$ rcond, rpvgrw, berr, n_err_bnds,
563  \$ err_bnds_norm, err_bnds_comp, nparams, params,
564  \$ work, rwork, info )
565 *
566 * -- LAPACK driver routine (version 3.4.1) --
567 * -- LAPACK is a software package provided by Univ. of Tennessee, --
568 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
569 * April 2012
570 *
571 * .. Scalar Arguments ..
572  CHARACTER EQUED, FACT, TRANS
573  INTEGER INFO, LDAB, LDAFB, LDB, LDX, N, NRHS, NPARAMS,
574  \$ n_err_bnds
575  REAL RCOND, RPVGRW
576 * ..
577 * .. Array Arguments ..
578  INTEGER IPIV( * )
579  COMPLEX AB( ldab, * ), AFB( ldafb, * ), B( ldb, * ),
580  \$ x( ldx , * ),work( * )
581  REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
582  \$ err_bnds_norm( nrhs, * ),
583  \$ err_bnds_comp( nrhs, * ), rwork( * )
584 * ..
585 *
586 * ==================================================================
587 *
588 * .. Parameters ..
589  REAL ZERO, ONE
590  parameter ( zero = 0.0e+0, one = 1.0e+0 )
591  INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
592  INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
593  INTEGER CMP_ERR_I, PIV_GROWTH_I
594  parameter ( final_nrm_err_i = 1, final_cmp_err_i = 2,
595  \$ berr_i = 3 )
596  parameter ( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
597  parameter ( cmp_rcond_i = 7, cmp_err_i = 8,
598  \$ piv_growth_i = 9 )
599 * ..
600 * .. Local Scalars ..
601  LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
602  INTEGER INFEQU, I, J, KL, KU
603  REAL AMAX, BIGNUM, COLCND, RCMAX, RCMIN,
604  \$ rowcnd, smlnum
605 * ..
606 * .. External Functions ..
607  EXTERNAL lsame, slamch, cla_gbrpvgrw
608  LOGICAL LSAME
609  REAL SLAMCH, CLA_GBRPVGRW
610 * ..
611 * .. External Subroutines ..
612  EXTERNAL cgbequb, cgbtrf, cgbtrs, clacpy, claqgb,
614 * ..
615 * .. Intrinsic Functions ..
616  INTRINSIC max, min
617 * ..
618 * .. Executable Statements ..
619 *
620  info = 0
621  nofact = lsame( fact, 'N' )
622  equil = lsame( fact, 'E' )
623  notran = lsame( trans, 'N' )
624  smlnum = slamch( 'Safe minimum' )
625  bignum = one / smlnum
626  IF( nofact .OR. equil ) THEN
627  equed = 'N'
628  rowequ = .false.
629  colequ = .false.
630  ELSE
631  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
632  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
633  END IF
634 *
635 * Default is failure. If an input parameter is wrong or
636 * factorization fails, make everything look horrible. Only the
637 * pivot growth is set here, the rest is initialized in CGBRFSX.
638 *
639  rpvgrw = zero
640 *
641 * Test the input parameters. PARAMS is not tested until SGERFSX.
642 *
643  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
644  \$ lsame( fact, 'F' ) ) THEN
645  info = -1
646  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
647  \$ lsame( trans, 'C' ) ) THEN
648  info = -2
649  ELSE IF( n.LT.0 ) THEN
650  info = -3
651  ELSE IF( kl.LT.0 ) THEN
652  info = -4
653  ELSE IF( ku.LT.0 ) THEN
654  info = -5
655  ELSE IF( nrhs.LT.0 ) THEN
656  info = -6
657  ELSE IF( ldab.LT.kl+ku+1 ) THEN
658  info = -8
659  ELSE IF( ldafb.LT.2*kl+ku+1 ) THEN
660  info = -10
661  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
662  \$ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
663  info = -12
664  ELSE
665  IF( rowequ ) THEN
666  rcmin = bignum
667  rcmax = zero
668  DO 10 j = 1, n
669  rcmin = min( rcmin, r( j ) )
670  rcmax = max( rcmax, r( j ) )
671  10 CONTINUE
672  IF( rcmin.LE.zero ) THEN
673  info = -13
674  ELSE IF( n.GT.0 ) THEN
675  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
676  ELSE
677  rowcnd = one
678  END IF
679  END IF
680  IF( colequ .AND. info.EQ.0 ) THEN
681  rcmin = bignum
682  rcmax = zero
683  DO 20 j = 1, n
684  rcmin = min( rcmin, c( j ) )
685  rcmax = max( rcmax, c( j ) )
686  20 CONTINUE
687  IF( rcmin.LE.zero ) THEN
688  info = -14
689  ELSE IF( n.GT.0 ) THEN
690  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
691  ELSE
692  colcnd = one
693  END IF
694  END IF
695  IF( info.EQ.0 ) THEN
696  IF( ldb.LT.max( 1, n ) ) THEN
697  info = -15
698  ELSE IF( ldx.LT.max( 1, n ) ) THEN
699  info = -16
700  END IF
701  END IF
702  END IF
703 *
704  IF( info.NE.0 ) THEN
705  CALL xerbla( 'CGBSVXX', -info )
706  RETURN
707  END IF
708 *
709  IF( equil ) THEN
710 *
711 * Compute row and column scalings to equilibrate the matrix A.
712 *
713  CALL cgbequb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,
714  \$ amax, infequ )
715  IF( infequ.EQ.0 ) THEN
716 *
717 * Equilibrate the matrix.
718 *
719  CALL claqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,
720  \$ amax, equed )
721  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
722  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
723  END IF
724 *
725 * If the scaling factors are not applied, set them to 1.0.
726 *
727  IF ( .NOT.rowequ ) THEN
728  DO j = 1, n
729  r( j ) = 1.0
730  END DO
731  END IF
732  IF ( .NOT.colequ ) THEN
733  DO j = 1, n
734  c( j ) = 1.0
735  END DO
736  END IF
737  END IF
738 *
739 * Scale the right-hand side.
740 *
741  IF( notran ) THEN
742  IF( rowequ ) CALL clascl2( n, nrhs, r, b, ldb )
743  ELSE
744  IF( colequ ) CALL clascl2( n, nrhs, c, b, ldb )
745  END IF
746 *
747  IF( nofact .OR. equil ) THEN
748 *
749 * Compute the LU factorization of A.
750 *
751  DO 40, j = 1, n
752  DO 30, i = kl+1, 2*kl+ku+1
753  afb( i, j ) = ab( i-kl, j )
754  30 CONTINUE
755  40 CONTINUE
756  CALL cgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
757 *
758 * Return if INFO is non-zero.
759 *
760  IF( info.GT.0 ) THEN
761 *
762 * Pivot in column INFO is exactly 0
763 * Compute the reciprocal pivot growth factor of the
764 * leading rank-deficient INFO columns of A.
765 *
766  rpvgrw = cla_gbrpvgrw( n, kl, ku, info, ab, ldab, afb,
767  \$ ldafb )
768  RETURN
769  END IF
770  END IF
771 *
772 * Compute the reciprocal pivot growth factor RPVGRW.
773 *
774  rpvgrw = cla_gbrpvgrw( n, kl, ku, n, ab, ldab, afb, ldafb )
775 *
776 * Compute the solution matrix X.
777 *
778  CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
779  CALL cgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,
780  \$ info )
781 *
782 * Use iterative refinement to improve the computed solution and
783 * compute error bounds and backward error estimates for it.
784 *
785  CALL cgbrfsx( trans, equed, n, kl, ku, nrhs, ab, ldab, afb, ldafb,
786  \$ ipiv, r, c, b, ldb, x, ldx, rcond, berr,
787  \$ n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params,
788  \$ work, rwork, info )
789
790 *
791 * Scale solutions.
792 *
793  IF ( colequ .AND. notran ) THEN
794  CALL clascl2( n, nrhs, c, x, ldx )
795  ELSE IF ( rowequ .AND. .NOT.notran ) THEN
796  CALL clascl2( n, nrhs, r, x, ldx )
797  END IF
798 *
799  RETURN
800 *
801 * End of CGBSVXX
802 *
803  END
subroutine cgbtrf(M, N, KL, KU, AB, LDAB, IPIV, INFO)
CGBTRF
Definition: cgbtrf.f:146
subroutine cgbequb(M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)
CGBEQUB
Definition: cgbequb.f:163
subroutine cgbrfsx(TRANS, EQUED, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO)
CGBRFSX
Definition: cgbrfsx.f:442
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine claqgb(M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, EQUED)
CLAQGB scales a general band matrix, using row and column scaling factors computed by sgbequ...
Definition: claqgb.f:162
subroutine clascl2(M, N, D, X, LDX)
CLASCL2 performs diagonal scaling on a vector.
Definition: clascl2.f:93
subroutine cgbsvxx(FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO)
CGBSVXX computes the solution to system of linear equations A * X = B for GB matrices ...
Definition: cgbsvxx.f:565
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
real function cla_gbrpvgrw(N, KL, KU, NCOLS, AB, LDAB, AFB, LDAFB)
CLA_GBRPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a general banded matrix...
Definition: cla_gbrpvgrw.f:119
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine cgbtrs(TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO)
CGBTRS
Definition: cgbtrs.f:140