LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
cgesvxx.f
Go to the documentation of this file.
1 *> \brief <b> CGESVXX computes the solution to system of linear equations A * X = B for GE 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 CGESVXX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvxx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvxx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvxx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGESVXX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
22 * EQUED, R, C, B, LDB, X, LDX, RCOND, RPVGRW,
23 * BERR, N_ERR_BNDS, ERR_BNDS_NORM,
24 * ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK,
25 * INFO )
26 *
27 * .. Scalar Arguments ..
28 * CHARACTER EQUED, FACT, TRANS
29 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
30 * $ N_ERR_BNDS
31 * REAL RCOND, RPVGRW
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IPIV( * )
35 * COMPLEX A( LDA, * ), AF( LDAF, * ), 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 *> CGESVXX 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. CGESVXX 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 *> CGESVXX 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 *> CGESVXX 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 CGESVXX 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)
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] NRHS
158 *> \verbatim
159 *> NRHS is INTEGER
160 *> The number of right hand sides, i.e., the number of columns
161 *> of the matrices B and X. NRHS >= 0.
162 *> \endverbatim
163 *>
164 *> \param[in,out] A
165 *> \verbatim
166 *> A is COMPLEX array, dimension (LDA,N)
167 *> On entry, the N-by-N matrix A. If FACT = 'F' and EQUED is
168 *> not 'N', then A must have been equilibrated by the scaling
169 *> factors in R and/or C. A is not modified if FACT = 'F' or
170 *> 'N', or if FACT = 'E' and EQUED = 'N' on exit.
171 *>
172 *> On exit, if EQUED .ne. 'N', A is scaled as follows:
173 *> EQUED = 'R': A := diag(R) * A
174 *> EQUED = 'C': A := A * diag(C)
175 *> EQUED = 'B': A := diag(R) * A * diag(C).
176 *> \endverbatim
177 *>
178 *> \param[in] LDA
179 *> \verbatim
180 *> LDA is INTEGER
181 *> The leading dimension of the array A. LDA >= max(1,N).
182 *> \endverbatim
183 *>
184 *> \param[in,out] AF
185 *> \verbatim
186 *> AF is COMPLEX array, dimension (LDAF,N)
187 *> If FACT = 'F', then AF is an input argument and on entry
188 *> contains the factors L and U from the factorization
189 *> A = P*L*U as computed by CGETRF. If EQUED .ne. 'N', then
190 *> AF is the factored form of the equilibrated matrix A.
191 *>
192 *> If FACT = 'N', then AF is an output argument and on exit
193 *> returns the factors L and U from the factorization A = P*L*U
194 *> of the original matrix A.
195 *>
196 *> If FACT = 'E', then AF is an output argument and on exit
197 *> returns the factors L and U from the factorization A = P*L*U
198 *> of the equilibrated matrix A (see the description of A for
199 *> the form of the equilibrated matrix).
200 *> \endverbatim
201 *>
202 *> \param[in] LDAF
203 *> \verbatim
204 *> LDAF is INTEGER
205 *> The leading dimension of the array AF. LDAF >= max(1,N).
206 *> \endverbatim
207 *>
208 *> \param[in,out] IPIV
209 *> \verbatim
210 *> IPIV is INTEGER array, dimension (N)
211 *> If FACT = 'F', then IPIV is an input argument and on entry
212 *> contains the pivot indices from the factorization A = P*L*U
213 *> as computed by CGETRF; row i of the matrix was interchanged
214 *> with row IPIV(i).
215 *>
216 *> If FACT = 'N', then IPIV is an output argument and on exit
217 *> contains the pivot indices from the factorization A = P*L*U
218 *> of the original matrix A.
219 *>
220 *> If FACT = 'E', then IPIV is an output argument and on exit
221 *> contains the pivot indices from the factorization A = P*L*U
222 *> of the equilibrated matrix A.
223 *> \endverbatim
224 *>
225 *> \param[in,out] EQUED
226 *> \verbatim
227 *> EQUED is CHARACTER*1
228 *> Specifies the form of equilibration that was done.
229 *> = 'N': No equilibration (always true if FACT = 'N').
230 *> = 'R': Row equilibration, i.e., A has been premultiplied by
231 *> diag(R).
232 *> = 'C': Column equilibration, i.e., A has been postmultiplied
233 *> by diag(C).
234 *> = 'B': Both row and column equilibration, i.e., A has been
235 *> replaced by diag(R) * A * diag(C).
236 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
237 *> output argument.
238 *> \endverbatim
239 *>
240 *> \param[in,out] R
241 *> \verbatim
242 *> R is REAL array, dimension (N)
243 *> The row scale factors for A. If EQUED = 'R' or 'B', A is
244 *> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
245 *> is not accessed. R is an input argument if FACT = 'F';
246 *> otherwise, R is an output argument. If FACT = 'F' and
247 *> EQUED = 'R' or 'B', each element of R must be positive.
248 *> If R is output, each element of R is a power of the radix.
249 *> If R is input, each element of R should be a power of the radix
250 *> to ensure a reliable solution and error estimates. Scaling by
251 *> powers of the radix does not cause rounding errors unless the
252 *> result underflows or overflows. Rounding errors during scaling
253 *> lead to refining with a matrix that is not equivalent to the
254 *> input matrix, producing error estimates that may not be
255 *> reliable.
256 *> \endverbatim
257 *>
258 *> \param[in,out] C
259 *> \verbatim
260 *> C is REAL array, dimension (N)
261 *> The column scale factors for A. If EQUED = 'C' or 'B', A is
262 *> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
263 *> is not accessed. C is an input argument if FACT = 'F';
264 *> otherwise, C is an output argument. If FACT = 'F' and
265 *> EQUED = 'C' or 'B', each element of C must be positive.
266 *> If C is output, each element of C is a power of the radix.
267 *> If C is input, each element of C should be a power of the radix
268 *> to ensure a reliable solution and error estimates. Scaling by
269 *> powers of the radix does not cause rounding errors unless the
270 *> result underflows or overflows. Rounding errors during scaling
271 *> lead to refining with a matrix that is not equivalent to the
272 *> input matrix, producing error estimates that may not be
273 *> reliable.
274 *> \endverbatim
275 *>
276 *> \param[in,out] B
277 *> \verbatim
278 *> B is COMPLEX array, dimension (LDB,NRHS)
279 *> On entry, the N-by-NRHS right hand side matrix B.
280 *> On exit,
281 *> if EQUED = 'N', B is not modified;
282 *> if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
283 *> diag(R)*B;
284 *> if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
285 *> overwritten by diag(C)*B.
286 *> \endverbatim
287 *>
288 *> \param[in] LDB
289 *> \verbatim
290 *> LDB is INTEGER
291 *> The leading dimension of the array B. LDB >= max(1,N).
292 *> \endverbatim
293 *>
294 *> \param[out] X
295 *> \verbatim
296 *> X is COMPLEX array, dimension (LDX,NRHS)
297 *> If INFO = 0, the N-by-NRHS solution matrix X to the original
298 *> system of equations. Note that A and B are modified on exit
299 *> if EQUED .ne. 'N', and the solution to the equilibrated system is
300 *> inv(diag(C))*X if TRANS = 'N' and EQUED = 'C' or 'B', or
301 *> inv(diag(R))*X if TRANS = 'T' or 'C' and EQUED = 'R' or 'B'.
302 *> \endverbatim
303 *>
304 *> \param[in] LDX
305 *> \verbatim
306 *> LDX is INTEGER
307 *> The leading dimension of the array X. LDX >= max(1,N).
308 *> \endverbatim
309 *>
310 *> \param[out] RCOND
311 *> \verbatim
312 *> RCOND is REAL
313 *> Reciprocal scaled condition number. This is an estimate of the
314 *> reciprocal Skeel condition number of the matrix A after
315 *> equilibration (if done). If this is less than the machine
316 *> precision (in particular, if it is zero), the matrix is singular
317 *> to working precision. Note that the error may still be small even
318 *> if this number is very small and the matrix appears ill-
319 *> conditioned.
320 *> \endverbatim
321 *>
322 *> \param[out] RPVGRW
323 *> \verbatim
324 *> RPVGRW is REAL
325 *> Reciprocal pivot growth. On exit, this contains the reciprocal
326 *> pivot growth factor norm(A)/norm(U). The "max absolute element"
327 *> norm is used. If this is much less than 1, then the stability of
328 *> the LU factorization of the (equilibrated) matrix A could be poor.
329 *> This also means that the solution X, estimated condition numbers,
330 *> and error bounds could be unreliable. If factorization fails with
331 *> 0<INFO<=N, then this contains the reciprocal pivot growth factor
332 *> for the leading INFO columns of A. In CGESVX, this quantity is
333 *> returned in WORK(1).
334 *> \endverbatim
335 *>
336 *> \param[out] BERR
337 *> \verbatim
338 *> BERR is REAL array, dimension (NRHS)
339 *> Componentwise relative backward error. This is the
340 *> componentwise relative backward error of each solution vector X(j)
341 *> (i.e., the smallest relative change in any element of A or B that
342 *> makes X(j) an exact solution).
343 *> \endverbatim
344 *>
345 *> \param[in] N_ERR_BNDS
346 *> \verbatim
347 *> N_ERR_BNDS is INTEGER
348 *> Number of error bounds to return for each right hand side
349 *> and each type (normwise or componentwise). See ERR_BNDS_NORM and
350 *> ERR_BNDS_COMP below.
351 *> \endverbatim
352 *>
353 *> \param[out] ERR_BNDS_NORM
354 *> \verbatim
355 *> ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
356 *> For each right-hand side, this array contains information about
357 *> various error bounds and condition numbers corresponding to the
358 *> normwise relative error, which is defined as follows:
359 *>
360 *> Normwise relative error in the ith solution vector:
361 *> max_j (abs(XTRUE(j,i) - X(j,i)))
362 *> ------------------------------
363 *> max_j abs(X(j,i))
364 *>
365 *> The array is indexed by the type of error information as described
366 *> below. There currently are up to three pieces of information
367 *> returned.
368 *>
369 *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
370 *> right-hand side.
371 *>
372 *> The second index in ERR_BNDS_NORM(:,err) contains the following
373 *> three fields:
374 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
375 *> reciprocal condition number is less than the threshold
376 *> sqrt(n) * slamch('Epsilon').
377 *>
378 *> err = 2 "Guaranteed" error bound: The estimated forward error,
379 *> almost certainly within a factor of 10 of the true error
380 *> so long as the next entry is greater than the threshold
381 *> sqrt(n) * slamch('Epsilon'). This error bound should only
382 *> be trusted if the previous boolean is true.
383 *>
384 *> err = 3 Reciprocal condition number: Estimated normwise
385 *> reciprocal condition number. Compared with the threshold
386 *> sqrt(n) * slamch('Epsilon') to determine if the error
387 *> estimate is "guaranteed". These reciprocal condition
388 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
389 *> appropriately scaled matrix Z.
390 *> Let Z = S*A, where S scales each row by a power of the
391 *> radix so all absolute row sums of Z are approximately 1.
392 *>
393 *> See Lapack Working Note 165 for further details and extra
394 *> cautions.
395 *> \endverbatim
396 *>
397 *> \param[out] ERR_BNDS_COMP
398 *> \verbatim
399 *> ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
400 *> For each right-hand side, this array contains information about
401 *> various error bounds and condition numbers corresponding to the
402 *> componentwise relative error, which is defined as follows:
403 *>
404 *> Componentwise relative error in the ith solution vector:
405 *> abs(XTRUE(j,i) - X(j,i))
406 *> max_j ----------------------
407 *> abs(X(j,i))
408 *>
409 *> The array is indexed by the right-hand side i (on which the
410 *> componentwise relative error depends), and the type of error
411 *> information as described below. There currently are up to three
412 *> pieces of information returned for each right-hand side. If
413 *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
414 *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS .LT. 3, then at most
415 *> the first (:,N_ERR_BNDS) entries are returned.
416 *>
417 *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
418 *> right-hand side.
419 *>
420 *> The second index in ERR_BNDS_COMP(:,err) contains the following
421 *> three fields:
422 *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
423 *> reciprocal condition number is less than the threshold
424 *> sqrt(n) * slamch('Epsilon').
425 *>
426 *> err = 2 "Guaranteed" error bound: The estimated forward error,
427 *> almost certainly within a factor of 10 of the true error
428 *> so long as the next entry is greater than the threshold
429 *> sqrt(n) * slamch('Epsilon'). This error bound should only
430 *> be trusted if the previous boolean is true.
431 *>
432 *> err = 3 Reciprocal condition number: Estimated componentwise
433 *> reciprocal condition number. Compared with the threshold
434 *> sqrt(n) * slamch('Epsilon') to determine if the error
435 *> estimate is "guaranteed". These reciprocal condition
436 *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
437 *> appropriately scaled matrix Z.
438 *> Let Z = S*(A*diag(x)), where x is the solution for the
439 *> current right-hand side and S scales each row of
440 *> A*diag(x) by a power of the radix so all absolute row
441 *> sums of Z are approximately 1.
442 *>
443 *> See Lapack Working Note 165 for further details and extra
444 *> cautions.
445 *> \endverbatim
446 *>
447 *> \param[in] NPARAMS
448 *> \verbatim
449 *> NPARAMS is INTEGER
450 *> Specifies the number of parameters set in PARAMS. If .LE. 0, the
451 *> PARAMS array is never referenced and default values are used.
452 *> \endverbatim
453 *>
454 *> \param[in,out] PARAMS
455 *> \verbatim
456 *> PARAMS is REAL array, dimension NPARAMS
457 *> Specifies algorithm parameters. If an entry is .LT. 0.0, then
458 *> that entry will be filled with default value used for that
459 *> parameter. Only positions up to NPARAMS are accessed; defaults
460 *> are used for higher-numbered parameters.
461 *>
462 *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
463 *> refinement or not.
464 *> Default: 1.0
465 *> = 0.0 : No refinement is performed, and no error bounds are
466 *> computed.
467 *> = 1.0 : Use the double-precision refinement algorithm,
468 *> possibly with doubled-single computations if the
469 *> compilation environment does not support DOUBLE
470 *> PRECISION.
471 *> (other values are reserved for future use)
472 *>
473 *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
474 *> computations allowed for refinement.
475 *> Default: 10
476 *> Aggressive: Set to 100 to permit convergence using approximate
477 *> factorizations or factorizations other than LU. If
478 *> the factorization uses a technique other than
479 *> Gaussian elimination, the guarantees in
480 *> err_bnds_norm and err_bnds_comp may no longer be
481 *> trustworthy.
482 *>
483 *> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
484 *> will attempt to find a solution with small componentwise
485 *> relative error in the double-precision algorithm. Positive
486 *> is true, 0.0 is false.
487 *> Default: 1.0 (attempt componentwise convergence)
488 *> \endverbatim
489 *>
490 *> \param[out] WORK
491 *> \verbatim
492 *> WORK is COMPLEX array, dimension (2*N)
493 *> \endverbatim
494 *>
495 *> \param[out] RWORK
496 *> \verbatim
497 *> RWORK is REAL array, dimension (2*N)
498 *> \endverbatim
499 *>
500 *> \param[out] INFO
501 *> \verbatim
502 *> INFO is INTEGER
503 *> = 0: Successful exit. The solution to every right-hand side is
504 *> guaranteed.
505 *> < 0: If INFO = -i, the i-th argument had an illegal value
506 *> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
507 *> has been completed, but the factor U is exactly singular, so
508 *> the solution and error bounds could not be computed. RCOND = 0
509 *> is returned.
510 *> = N+J: The solution corresponding to the Jth right-hand side is
511 *> not guaranteed. The solutions corresponding to other right-
512 *> hand sides K with K > J may not be guaranteed as well, but
513 *> only the first such right-hand side is reported. If a small
514 *> componentwise error is not requested (PARAMS(3) = 0.0) then
515 *> the Jth right-hand side is the first with a normwise error
516 *> bound that is not guaranteed (the smallest J such
517 *> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
518 *> the Jth right-hand side is the first with either a normwise or
519 *> componentwise error bound that is not guaranteed (the smallest
520 *> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
521 *> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
522 *> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
523 *> about all of the right-hand sides check ERR_BNDS_NORM or
524 *> ERR_BNDS_COMP.
525 *> \endverbatim
526 *
527 * Authors:
528 * ========
529 *
530 *> \author Univ. of Tennessee
531 *> \author Univ. of California Berkeley
532 *> \author Univ. of Colorado Denver
533 *> \author NAG Ltd.
534 *
535 *> \date April 2012
536 *
537 *> \ingroup complexGEsolve
538 *
539 * =====================================================================
540  SUBROUTINE cgesvxx( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV,
541  $ equed, r, c, b, ldb, x, ldx, rcond, rpvgrw,
542  $ berr, n_err_bnds, err_bnds_norm,
543  $ err_bnds_comp, nparams, params, work, rwork,
544  $ info )
545 *
546 * -- LAPACK driver routine (version 3.4.1) --
547 * -- LAPACK is a software package provided by Univ. of Tennessee, --
548 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
549 * April 2012
550 *
551 * .. Scalar Arguments ..
552  CHARACTER EQUED, FACT, TRANS
553  INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
554  $ n_err_bnds
555  REAL RCOND, RPVGRW
556 * ..
557 * .. Array Arguments ..
558  INTEGER IPIV( * )
559  COMPLEX A( lda, * ), AF( ldaf, * ), B( ldb, * ),
560  $ x( ldx , * ),work( * )
561  REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
562  $ err_bnds_norm( nrhs, * ),
563  $ err_bnds_comp( nrhs, * ), rwork( * )
564 * ..
565 *
566 * ==================================================================
567 *
568 * .. Parameters ..
569  REAL ZERO, ONE
570  parameter ( zero = 0.0e+0, one = 1.0e+0 )
571  INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
572  INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
573  INTEGER CMP_ERR_I, PIV_GROWTH_I
574  parameter ( final_nrm_err_i = 1, final_cmp_err_i = 2,
575  $ berr_i = 3 )
576  parameter ( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
577  parameter ( cmp_rcond_i = 7, cmp_err_i = 8,
578  $ piv_growth_i = 9 )
579 * ..
580 * .. Local Scalars ..
581  LOGICAL COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
582  INTEGER INFEQU, J
583  REAL AMAX, BIGNUM, COLCND, RCMAX, RCMIN,
584  $ rowcnd, smlnum
585 * ..
586 * .. External Functions ..
587  EXTERNAL lsame, slamch, cla_gerpvgrw
588  LOGICAL LSAME
589  REAL SLAMCH, CLA_GERPVGRW
590 * ..
591 * .. External Subroutines ..
592  EXTERNAL cgeequb, cgetrf, cgetrs, clacpy, claqge,
594 * ..
595 * .. Intrinsic Functions ..
596  INTRINSIC max, min
597 * ..
598 * .. Executable Statements ..
599 *
600  info = 0
601  nofact = lsame( fact, 'N' )
602  equil = lsame( fact, 'E' )
603  notran = lsame( trans, 'N' )
604  smlnum = slamch( 'Safe minimum' )
605  bignum = one / smlnum
606  IF( nofact .OR. equil ) THEN
607  equed = 'N'
608  rowequ = .false.
609  colequ = .false.
610  ELSE
611  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
612  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
613  END IF
614 *
615 * Default is failure. If an input parameter is wrong or
616 * factorization fails, make everything look horrible. Only the
617 * pivot growth is set here, the rest is initialized in CGERFSX.
618 *
619  rpvgrw = zero
620 *
621 * Test the input parameters. PARAMS is not tested until CGERFSX.
622 *
623  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
624  $ lsame( fact, 'F' ) ) THEN
625  info = -1
626  ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
627  $ lsame( trans, 'C' ) ) THEN
628  info = -2
629  ELSE IF( n.LT.0 ) THEN
630  info = -3
631  ELSE IF( nrhs.LT.0 ) THEN
632  info = -4
633  ELSE IF( lda.LT.max( 1, n ) ) THEN
634  info = -6
635  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
636  info = -8
637  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
638  $ ( rowequ .OR. colequ .OR. lsame( equed, 'N' ) ) ) THEN
639  info = -10
640  ELSE
641  IF( rowequ ) THEN
642  rcmin = bignum
643  rcmax = zero
644  DO 10 j = 1, n
645  rcmin = min( rcmin, r( j ) )
646  rcmax = max( rcmax, r( j ) )
647  10 CONTINUE
648  IF( rcmin.LE.zero ) THEN
649  info = -11
650  ELSE IF( n.GT.0 ) THEN
651  rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
652  ELSE
653  rowcnd = one
654  END IF
655  END IF
656  IF( colequ .AND. info.EQ.0 ) THEN
657  rcmin = bignum
658  rcmax = zero
659  DO 20 j = 1, n
660  rcmin = min( rcmin, c( j ) )
661  rcmax = max( rcmax, c( j ) )
662  20 CONTINUE
663  IF( rcmin.LE.zero ) THEN
664  info = -12
665  ELSE IF( n.GT.0 ) THEN
666  colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
667  ELSE
668  colcnd = one
669  END IF
670  END IF
671  IF( info.EQ.0 ) THEN
672  IF( ldb.LT.max( 1, n ) ) THEN
673  info = -14
674  ELSE IF( ldx.LT.max( 1, n ) ) THEN
675  info = -16
676  END IF
677  END IF
678  END IF
679 *
680  IF( info.NE.0 ) THEN
681  CALL xerbla( 'CGESVXX', -info )
682  RETURN
683  END IF
684 *
685  IF( equil ) THEN
686 *
687 * Compute row and column scalings to equilibrate the matrix A.
688 *
689  CALL cgeequb( n, n, a, lda, r, c, rowcnd, colcnd, amax,
690  $ infequ )
691  IF( infequ.EQ.0 ) THEN
692 *
693 * Equilibrate the matrix.
694 *
695  CALL claqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,
696  $ equed )
697  rowequ = lsame( equed, 'R' ) .OR. lsame( equed, 'B' )
698  colequ = lsame( equed, 'C' ) .OR. lsame( equed, 'B' )
699  END IF
700 *
701 * If the scaling factors are not applied, set them to 1.0.
702 *
703  IF ( .NOT.rowequ ) THEN
704  DO j = 1, n
705  r( j ) = 1.0
706  END DO
707  END IF
708  IF ( .NOT.colequ ) THEN
709  DO j = 1, n
710  c( j ) = 1.0
711  END DO
712  END IF
713  END IF
714 *
715 * Scale the right-hand side.
716 *
717  IF( notran ) THEN
718  IF( rowequ ) CALL clascl2( n, nrhs, r, b, ldb )
719  ELSE
720  IF( colequ ) CALL clascl2( n, nrhs, c, b, ldb )
721  END IF
722 *
723  IF( nofact .OR. equil ) THEN
724 *
725 * Compute the LU factorization of A.
726 *
727  CALL clacpy( 'Full', n, n, a, lda, af, ldaf )
728  CALL cgetrf( n, n, af, ldaf, ipiv, info )
729 *
730 * Return if INFO is non-zero.
731 *
732  IF( info.GT.0 ) THEN
733 *
734 * Pivot in column INFO is exactly 0
735 * Compute the reciprocal pivot growth factor of the
736 * leading rank-deficient INFO columns of A.
737 *
738  rpvgrw = cla_gerpvgrw( n, info, a, lda, af, ldaf )
739  RETURN
740  END IF
741  END IF
742 *
743 * Compute the reciprocal pivot growth factor RPVGRW.
744 *
745  rpvgrw = cla_gerpvgrw( n, n, a, lda, af, ldaf )
746 *
747 * Compute the solution matrix X.
748 *
749  CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
750  CALL cgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
751 *
752 * Use iterative refinement to improve the computed solution and
753 * compute error bounds and backward error estimates for it.
754 *
755  CALL cgerfsx( trans, equed, n, nrhs, a, lda, af, ldaf,
756  $ ipiv, r, c, b, ldb, x, ldx, rcond, berr,
757  $ n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params,
758  $ work, rwork, info )
759 *
760 * Scale solutions.
761 *
762  IF ( colequ .AND. notran ) THEN
763  CALL clascl2 ( n, nrhs, c, x, ldx )
764  ELSE IF ( rowequ .AND. .NOT.notran ) THEN
765  CALL clascl2 ( n, nrhs, r, x, ldx )
766  END IF
767 *
768  RETURN
769 *
770 * End of CGESVXX
771 *
772  END
subroutine cgerfsx(TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV, R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, RWORK, INFO)
CGERFSX
Definition: cgerfsx.f:416
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:123
real function cla_gerpvgrw(N, NCOLS, A, LDA, AF, LDAF)
CLA_GERPVGRW multiplies a square real matrix by a complex matrix.
Definition: cla_gerpvgrw.f:100
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgeequb(M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
CGEEQUB
Definition: cgeequb.f:149
subroutine claqge(M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, EQUED)
CLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ...
Definition: claqge.f:145
subroutine clascl2(M, N, D, X, LDX)
CLASCL2 performs diagonal scaling on a vector.
Definition: clascl2.f:93
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:110
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
subroutine cgesvxx(FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, 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)
CGESVXX computes the solution to system of linear equations A * X = B for GE matrices ...
Definition: cgesvxx.f:545
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55