1 *> \brief <b> SPOSVX computes the solution to system of linear equations A * X = B for PO matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SPOSVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sposvx.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sposvx.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sposvx.f">
21 * SUBROUTINE SPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
22 * S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
25 * .. Scalar Arguments ..
26 * CHARACTER EQUED, FACT, UPLO
27 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
30 * .. Array Arguments ..
32 * REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
33 * $ BERR( * ), FERR( * ), S( * ), WORK( * ),
43 *> SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
44 *> compute the solution to a real system of linear equations
46 *> where A is an N-by-N symmetric positive definite matrix and X and B
47 *> are N-by-NRHS matrices.
49 *> Error bounds on the solution and a condition estimate are also
58 *> The following steps are performed:
60 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate
62 *> diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
63 *> Whether or not the system will be equilibrated depends on the
64 *> scaling of the matrix A, but if equilibration is used, A is
65 *> overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
67 *> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
68 *> factor the matrix A (after equilibration if FACT = 'E') as
69 *> A = U**T* U, if UPLO = 'U', or
70 *> A = L * L**T, if UPLO = 'L',
71 *> where U is an upper triangular matrix and L is a lower triangular
74 *> 3. If the leading i-by-i principal minor is not positive definite,
75 *> then the routine returns with INFO = i. Otherwise, the factored
76 *> form of A is used to estimate the condition number of the matrix
77 *> A. If the reciprocal of the condition number is less than machine
78 *> precision, INFO = N+1 is returned as a warning, but the routine
79 *> still goes on to solve for X and compute error bounds as
82 *> 4. The system of equations is solved for X using the factored form
85 *> 5. Iterative refinement is applied to improve the computed solution
86 *> matrix and calculate error bounds and backward error estimates
89 *> 6. If equilibration was used, the matrix X is premultiplied by
90 *> diag(S) so that it solves the original system before
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 contains the factored form of A.
104 *> If EQUED = 'Y', the matrix A has been equilibrated
105 *> with scaling factors given by S. A and AF will not
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.
114 *> UPLO is CHARACTER*1
115 *> = 'U': Upper triangle of A is stored;
116 *> = 'L': Lower triangle of A is stored.
122 *> The number of linear equations, i.e., the order of the
129 *> The number of right hand sides, i.e., the number of columns
130 *> of the matrices B and X. NRHS >= 0.
135 *> A is REAL array, dimension (LDA,N)
136 *> On entry, the symmetric matrix A, except if FACT = 'F' and
137 *> EQUED = 'Y', then A must contain the equilibrated matrix
138 *> diag(S)*A*diag(S). If UPLO = 'U', the leading
139 *> N-by-N upper triangular part of A contains the upper
140 *> triangular part of the matrix A, and the strictly lower
141 *> triangular part of A is not referenced. If UPLO = 'L', the
142 *> leading N-by-N lower triangular part of A contains the lower
143 *> triangular part of the matrix A, and the strictly upper
144 *> triangular part of A is not referenced. A is not modified if
145 *> FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
147 *> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
148 *> diag(S)*A*diag(S).
154 *> The leading dimension of the array A. LDA >= max(1,N).
159 *> AF is REAL array, dimension (LDAF,N)
160 *> If FACT = 'F', then AF is an input argument and on entry
161 *> contains the triangular factor U or L from the Cholesky
162 *> factorization A = U**T*U or A = L*L**T, in the same storage
163 *> format as A. If EQUED .ne. 'N', then AF is the factored form
164 *> of the equilibrated matrix diag(S)*A*diag(S).
166 *> If FACT = 'N', then AF is an output argument and on exit
167 *> returns the triangular factor U or L from the Cholesky
168 *> factorization A = U**T*U or A = L*L**T of the original
171 *> If FACT = 'E', then AF is an output argument and on exit
172 *> returns the triangular factor U or L from the Cholesky
173 *> factorization A = U**T*U or A = L*L**T of the equilibrated
174 *> matrix A (see the description of A for the form of the
175 *> equilibrated matrix).
181 *> The leading dimension of the array AF. LDAF >= max(1,N).
184 *> \param[in,out] EQUED
186 *> EQUED is CHARACTER*1
187 *> Specifies the form of equilibration that was done.
188 *> = 'N': No equilibration (always true if FACT = 'N').
189 *> = 'Y': Equilibration was done, i.e., A has been replaced by
190 *> diag(S) * A * diag(S).
191 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
197 *> S is REAL array, dimension (N)
198 *> The scale factors for A; not accessed if EQUED = 'N'. S is
199 *> an input argument if FACT = 'F'; otherwise, S is an output
200 *> argument. If FACT = 'F' and EQUED = 'Y', each element of S
206 *> B is REAL array, dimension (LDB,NRHS)
207 *> On entry, the N-by-NRHS right hand side matrix B.
208 *> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
209 *> B is overwritten by diag(S) * B.
215 *> The leading dimension of the array B. LDB >= max(1,N).
220 *> X is REAL array, dimension (LDX,NRHS)
221 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
222 *> the original system of equations. Note that if EQUED = 'Y',
223 *> A and B are modified on exit, and the solution to the
224 *> equilibrated system is inv(diag(S))*X.
230 *> The leading dimension of the array X. LDX >= max(1,N).
236 *> The estimate of the reciprocal condition number of the matrix
237 *> A after equilibration (if done). If RCOND is less than the
238 *> machine precision (in particular, if RCOND = 0), the matrix
239 *> is singular to working precision. This condition is
240 *> indicated by a return code of INFO > 0.
245 *> FERR is REAL array, dimension (NRHS)
246 *> The estimated forward error bound for each solution vector
247 *> X(j) (the j-th column of the solution matrix X).
248 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
249 *> is an estimated upper bound for the magnitude of the largest
250 *> element in (X(j) - XTRUE) divided by the magnitude of the
251 *> largest element in X(j). The estimate is as reliable as
252 *> the estimate for RCOND, and is almost always a slight
253 *> overestimate of the true error.
258 *> BERR is REAL array, dimension (NRHS)
259 *> The componentwise relative backward error of each solution
260 *> vector X(j) (i.e., the smallest relative change in
261 *> any element of A or B that makes X(j) an exact solution).
266 *> WORK is REAL array, dimension (3*N)
271 *> IWORK is INTEGER array, dimension (N)
277 *> = 0: successful exit
278 *> < 0: if INFO = -i, the i-th argument had an illegal value
279 *> > 0: if INFO = i, and i is
280 *> <= N: the leading minor of order i of A is
281 *> not positive definite, so the factorization
282 *> could not be completed, and the solution has not
283 *> been computed. RCOND = 0 is returned.
284 *> = N+1: U is nonsingular, but RCOND is less than machine
285 *> precision, meaning that the matrix is singular
286 *> to working precision. Nevertheless, the
287 *> solution and error bounds are computed because
288 *> there are a number of situations where the
289 *> computed solution can be more accurate than the
290 *> value of RCOND would suggest.
296 *> \author Univ. of Tennessee
297 *> \author Univ. of California Berkeley
298 *> \author Univ. of Colorado Denver
303 *> \ingroup realPOsolve
305 * =====================================================================
306 SUBROUTINE SPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED,
307 $ S, B, LDB, X, LDX, RCOND, FERR, BERR, WORK,
310 * -- LAPACK driver routine (version 3.4.1) --
311 * -- LAPACK is a software package provided by Univ. of Tennessee, --
312 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
315 * .. Scalar Arguments ..
316 CHARACTER EQUED, FACT, UPLO
317 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
320 * .. Array Arguments ..
322 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
323 $ BERR( * ), FERR( * ), S( * ), WORK( * ),
327 * =====================================================================
331 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
333 * .. Local Scalars ..
334 LOGICAL EQUIL, NOFACT, RCEQU
336 REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
338 * .. External Functions ..
341 EXTERNAL LSAME, SLAMCH, SLANSY
343 * .. External Subroutines ..
344 EXTERNAL SLACPY, SLAQSY, SPOCON, SPOEQU, SPORFS, SPOTRF,
347 * .. Intrinsic Functions ..
350 * .. Executable Statements ..
353 NOFACT = LSAME( FACT, 'N' )
354 EQUIL = LSAME( FACT, 'E' )
355 IF( NOFACT .OR. EQUIL ) THEN
359 RCEQU = LSAME( EQUED, 'Y' )
360 SMLNUM = SLAMCH( 'Safe minimum' )
361 BIGNUM = ONE / SMLNUM
364 * Test the input parameters.
366 IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
369 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
372 ELSE IF( N.LT.0 ) THEN
374 ELSE IF( NRHS.LT.0 ) THEN
376 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
378 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
380 ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
381 $ ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
388 SMIN = MIN( SMIN, S( J ) )
389 SMAX = MAX( SMAX, S( J ) )
391 IF( SMIN.LE.ZERO ) THEN
393 ELSE IF( N.GT.0 ) THEN
394 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
400 IF( LDB.LT.MAX( 1, N ) ) THEN
402 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
409 CALL XERBLA( 'SPOSVX', -INFO )
415 * Compute row and column scalings to equilibrate the matrix A.
417 CALL SPOEQU( N, A, LDA, S, SCOND, AMAX, INFEQU )
418 IF( INFEQU.EQ.0 ) THEN
420 * Equilibrate the matrix.
422 CALL SLAQSY( UPLO, N, A, LDA, S, SCOND, AMAX, EQUED )
423 RCEQU = LSAME( EQUED, 'Y' )
427 * Scale the right hand side.
432 B( I, J ) = S( I )*B( I, J )
437 IF( NOFACT .OR. EQUIL ) THEN
439 * Compute the Cholesky factorization A = U**T *U or A = L*L**T.
441 CALL SLACPY( UPLO, N, N, A, LDA, AF, LDAF )
442 CALL SPOTRF( UPLO, N, AF, LDAF, INFO )
444 * Return if INFO is non-zero.
452 * Compute the norm of the matrix A.
454 ANORM = SLANSY( '1', UPLO, N, A, LDA, WORK )
456 * Compute the reciprocal of the condition number of A.
458 CALL SPOCON( UPLO, N, AF, LDAF, ANORM, RCOND, WORK, IWORK, INFO )
460 * Compute the solution matrix X.
462 CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
463 CALL SPOTRS( UPLO, N, NRHS, AF, LDAF, X, LDX, INFO )
465 * Use iterative refinement to improve the computed solution and
466 * compute error bounds and backward error estimates for it.
468 CALL SPORFS( UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX,
469 $ FERR, BERR, WORK, IWORK, INFO )
471 * Transform the solution matrix X to a solution of the original
477 X( I, J ) = S( I )*X( I, J )
481 FERR( J ) = FERR( J ) / SCOND
485 * Set INFO = N+1 if the matrix is singular to working precision.
487 IF( RCOND.LT.SLAMCH( 'Epsilon' ) )