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