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