LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ sposvx()

 subroutine sposvx ( character FACT, character UPLO, integer N, integer NRHS, real, dimension( lda, * ) A, integer LDA, real, dimension( ldaf, * ) AF, integer LDAF, character EQUED, real, dimension( * ) S, real, dimension( ldb, * ) B, integer LDB, real, dimension( ldx, * ) X, integer LDX, real RCOND, real, dimension( * ) FERR, real, dimension( * ) BERR, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO )

SPOSVX computes the solution to system of linear equations A * X = B for PO matrices

Purpose:
``` SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
compute the solution to a real system of linear equations
A * X = B,
where A is an N-by-N symmetric positive definite matrix and X and B
are N-by-NRHS matrices.

Error bounds on the solution and a condition estimate are also
provided.```
Description:
``` The following steps are performed:

1. If FACT = 'E', real scaling factors are computed to equilibrate
the system:
diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
Whether or not the system will be equilibrated depends on the
scaling of the matrix A, but if equilibration is used, A is
overwritten by diag(S)*A*diag(S) and B by diag(S)*B.

2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
factor the matrix A (after equilibration if FACT = 'E') as
A = U**T* U,  if UPLO = 'U', or
A = L * L**T,  if UPLO = 'L',
where U is an upper triangular matrix and L is a lower triangular
matrix.

3. If the leading i-by-i principal minor is not positive definite,
then the routine returns with INFO = i. Otherwise, the factored
form of A is used to estimate the condition number of the matrix
A.  If the reciprocal of the condition number is less than machine
precision, INFO = N+1 is returned as a warning, but the routine
still goes on to solve for X and compute error bounds as
described below.

4. The system of equations is solved for X using the factored form
of A.

5. Iterative refinement is applied to improve the computed solution
matrix and calculate error bounds and backward error estimates
for it.

6. If equilibration was used, the matrix X is premultiplied by
diag(S) so that it solves the original system before
equilibration.```
Parameters
 [in] FACT ``` FACT is CHARACTER*1 Specifies whether or not the factored form of the matrix A is supplied on entry, and if not, whether the matrix A should be equilibrated before it is factored. = 'F': On entry, AF contains the factored form of A. If EQUED = 'Y', the matrix A has been equilibrated with scaling factors given by S. A and AF will not be modified. = 'N': The matrix A will be copied to AF and factored. = 'E': The matrix A will be equilibrated if necessary, then copied to AF and factored.``` [in] UPLO ``` UPLO is CHARACTER*1 = 'U': Upper triangle of A is stored; = 'L': Lower triangle of A is stored.``` [in] N ``` N is INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0.``` [in] NRHS ``` NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrices B and X. NRHS >= 0.``` [in,out] A ``` A is REAL array, dimension (LDA,N) On entry, the symmetric matrix A, except if FACT = 'F' and EQUED = 'Y', then A must contain the equilibrated matrix diag(S)*A*diag(S). If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. A is not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit. On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by diag(S)*A*diag(S).``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [in,out] AF ``` AF is REAL array, dimension (LDAF,N) If FACT = 'F', then AF is an input argument and on entry contains the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T, in the same storage format as A. If EQUED .ne. 'N', then AF is the factored form of the equilibrated matrix diag(S)*A*diag(S). If FACT = 'N', then AF is an output argument and on exit returns the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T of the original matrix A. If FACT = 'E', then AF is an output argument and on exit returns the triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T of the equilibrated matrix A (see the description of A for the form of the equilibrated matrix).``` [in] LDAF ``` LDAF is INTEGER The leading dimension of the array AF. LDAF >= max(1,N).``` [in,out] EQUED ``` EQUED is CHARACTER*1 Specifies the form of equilibration that was done. = 'N': No equilibration (always true if FACT = 'N'). = 'Y': Equilibration was done, i.e., A has been replaced by diag(S) * A * diag(S). EQUED is an input argument if FACT = 'F'; otherwise, it is an output argument.``` [in,out] S ``` S is REAL array, dimension (N) The scale factors for A; not accessed if EQUED = 'N'. S is an input argument if FACT = 'F'; otherwise, S is an output argument. If FACT = 'F' and EQUED = 'Y', each element of S must be positive.``` [in,out] B ``` B is REAL array, dimension (LDB,NRHS) On entry, the N-by-NRHS right hand side matrix B. On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten by diag(S) * B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [out] X ``` X is REAL array, dimension (LDX,NRHS) If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to the original system of equations. Note that if EQUED = 'Y', A and B are modified on exit, and the solution to the equilibrated system is inv(diag(S))*X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] RCOND ``` RCOND is REAL The estimate of the reciprocal condition number of the matrix A after equilibration (if done). If RCOND is less than the machine precision (in particular, if RCOND = 0), the matrix is singular to working precision. This condition is indicated by a return code of INFO > 0.``` [out] FERR ``` FERR is REAL array, dimension (NRHS) The estimated forward error bound for each solution vector X(j) (the j-th column of the solution matrix X). If XTRUE is the true solution corresponding to X(j), FERR(j) is an estimated upper bound for the magnitude of the largest element in (X(j) - XTRUE) divided by the magnitude of the largest element in X(j). The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error.``` [out] BERR ``` BERR is REAL array, dimension (NRHS) The componentwise relative backward error of each solution vector X(j) (i.e., the smallest relative change in any element of A or B that makes X(j) an exact solution).``` [out] WORK ` WORK is REAL array, dimension (3*N)` [out] IWORK ` IWORK is INTEGER array, dimension (N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, and i is <= N: the leading minor of order i of A is not positive definite, so the factorization could not be completed, and the solution has not been computed. RCOND = 0 is returned. = N+1: U is nonsingular, but RCOND is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest.```

Definition at line 304 of file sposvx.f.

307*
308* -- LAPACK driver routine --
309* -- LAPACK is a software package provided by Univ. of Tennessee, --
310* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
311*
312* .. Scalar Arguments ..
313 CHARACTER EQUED, FACT, UPLO
314 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
315 REAL RCOND
316* ..
317* .. Array Arguments ..
318 INTEGER IWORK( * )
319 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
320 \$ BERR( * ), FERR( * ), S( * ), WORK( * ),
321 \$ X( LDX, * )
322* ..
323*
324* =====================================================================
325*
326* .. Parameters ..
327 REAL ZERO, ONE
328 parameter( zero = 0.0e+0, one = 1.0e+0 )
329* ..
330* .. Local Scalars ..
331 LOGICAL EQUIL, NOFACT, RCEQU
332 INTEGER I, INFEQU, J
333 REAL AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
334* ..
335* .. External Functions ..
336 LOGICAL LSAME
337 REAL SLAMCH, SLANSY
338 EXTERNAL lsame, slamch, slansy
339* ..
340* .. External Subroutines ..
341 EXTERNAL slacpy, slaqsy, spocon, spoequ, sporfs, spotrf,
342 \$ spotrs, xerbla
343* ..
344* .. Intrinsic Functions ..
345 INTRINSIC max, min
346* ..
347* .. Executable Statements ..
348*
349 info = 0
350 nofact = lsame( fact, 'N' )
351 equil = lsame( fact, 'E' )
352 IF( nofact .OR. equil ) THEN
353 equed = 'N'
354 rcequ = .false.
355 ELSE
356 rcequ = lsame( equed, 'Y' )
357 smlnum = slamch( 'Safe minimum' )
358 bignum = one / smlnum
359 END IF
360*
361* Test the input parameters.
362*
363 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
364 \$ THEN
365 info = -1
366 ELSE IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) )
367 \$ THEN
368 info = -2
369 ELSE IF( n.LT.0 ) THEN
370 info = -3
371 ELSE IF( nrhs.LT.0 ) THEN
372 info = -4
373 ELSE IF( lda.LT.max( 1, n ) ) THEN
374 info = -6
375 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
376 info = -8
377 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
378 \$ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
379 info = -9
380 ELSE
381 IF( rcequ ) THEN
382 smin = bignum
383 smax = zero
384 DO 10 j = 1, n
385 smin = min( smin, s( j ) )
386 smax = max( smax, s( j ) )
387 10 CONTINUE
388 IF( smin.LE.zero ) THEN
389 info = -10
390 ELSE IF( n.GT.0 ) THEN
391 scond = max( smin, smlnum ) / min( smax, bignum )
392 ELSE
393 scond = one
394 END IF
395 END IF
396 IF( info.EQ.0 ) THEN
397 IF( ldb.LT.max( 1, n ) ) THEN
398 info = -12
399 ELSE IF( ldx.LT.max( 1, n ) ) THEN
400 info = -14
401 END IF
402 END IF
403 END IF
404*
405 IF( info.NE.0 ) THEN
406 CALL xerbla( 'SPOSVX', -info )
407 RETURN
408 END IF
409*
410 IF( equil ) THEN
411*
412* Compute row and column scalings to equilibrate the matrix A.
413*
414 CALL spoequ( n, a, lda, s, scond, amax, infequ )
415 IF( infequ.EQ.0 ) THEN
416*
417* Equilibrate the matrix.
418*
419 CALL slaqsy( uplo, n, a, lda, s, scond, amax, equed )
420 rcequ = lsame( equed, 'Y' )
421 END IF
422 END IF
423*
424* Scale the right hand side.
425*
426 IF( rcequ ) THEN
427 DO 30 j = 1, nrhs
428 DO 20 i = 1, n
429 b( i, j ) = s( i )*b( i, j )
430 20 CONTINUE
431 30 CONTINUE
432 END IF
433*
434 IF( nofact .OR. equil ) THEN
435*
436* Compute the Cholesky factorization A = U**T *U or A = L*L**T.
437*
438 CALL slacpy( uplo, n, n, a, lda, af, ldaf )
439 CALL spotrf( uplo, n, af, ldaf, info )
440*
441* Return if INFO is non-zero.
442*
443 IF( info.GT.0 )THEN
444 rcond = zero
445 RETURN
446 END IF
447 END IF
448*
449* Compute the norm of the matrix A.
450*
451 anorm = slansy( '1', uplo, n, a, lda, work )
452*
453* Compute the reciprocal of the condition number of A.
454*
455 CALL spocon( uplo, n, af, ldaf, anorm, rcond, work, iwork, info )
456*
457* Compute the solution matrix X.
458*
459 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
460 CALL spotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
461*
462* Use iterative refinement to improve the computed solution and
463* compute error bounds and backward error estimates for it.
464*
465 CALL sporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,
466 \$ ferr, berr, work, iwork, info )
467*
468* Transform the solution matrix X to a solution of the original
469* system.
470*
471 IF( rcequ ) THEN
472 DO 50 j = 1, nrhs
473 DO 40 i = 1, n
474 x( i, j ) = s( i )*x( i, j )
475 40 CONTINUE
476 50 CONTINUE
477 DO 60 j = 1, nrhs
478 ferr( j ) = ferr( j ) / scond
479 60 CONTINUE
480 END IF
481*
482* Set INFO = N+1 if the matrix is singular to working precision.
483*
484 IF( rcond.LT.slamch( 'Epsilon' ) )
485 \$ info = n + 1
486*
487 RETURN
488*
489* End of SPOSVX
490*
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine spotrf(UPLO, N, A, LDA, INFO)
SPOTRF
Definition: spotrf.f:107
subroutine spocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SPOCON
Definition: spocon.f:121
subroutine sporfs(UPLO, N, NRHS, A, LDA, AF, LDAF, B, LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO)
SPORFS
Definition: sporfs.f:183
subroutine spoequ(N, A, LDA, S, SCOND, AMAX, INFO)
SPOEQU
Definition: spoequ.f:112
subroutine spotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
SPOTRS
Definition: spotrs.f:110
real function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: slansy.f:122
subroutine slaqsy(UPLO, N, A, LDA, S, SCOND, AMAX, EQUED)
SLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
Definition: slaqsy.f:133
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: