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

## ◆ sposvxx()

 subroutine sposvxx ( 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 RPVGRW, real, dimension( * ) BERR, integer N_ERR_BNDS, real, dimension( nrhs, * ) ERR_BNDS_NORM, real, dimension( nrhs, * ) ERR_BNDS_COMP, integer NPARAMS, real, dimension( * ) PARAMS, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO )

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

Purpose:
SPOSVXX 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.

If requested, both normwise and maximum componentwise error bounds
are returned. SPOSVXX will return a solution with a tiny
guaranteed error (O(eps) where eps is the working machine
precision) unless the matrix is very ill-conditioned, in which
case a warning is returned. Relevant condition numbers also are
calculated and returned.

SPOSVXX accepts user-provided factorizations and equilibration
factors; see the definitions of the FACT and EQUED options.
Solving with refinement and using a factorization from a previous
SPOSVXX call will also produce a solution with either O(eps)
errors or warnings, but we cannot make that claim for general
user-provided factorizations and equilibration factors if they
differ from what SPOSVXX would itself produce.
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 (see argument RCOND).  If the reciprocal of the condition number
is less than machine precision, 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. By default (unless PARAMS(LA_LINRX_ITREF_I) is set to zero),
the routine will use iterative refinement to try to get a small
error and error bounds.  Refinement calculates the residual to at
least twice the working precision.

6. If equilibration was used, the matrix X is premultiplied by
diag(S) so that it solves the original system before
equilibration.
Some optional parameters are bundled in the PARAMS array.  These
settings determine how refinement is performed, but often the
defaults are acceptable.  If the defaults are acceptable, users
can pass NPARAMS = 0 which prevents the source code from accessing
the PARAMS argument.
Parameters

Definition at line 493 of file sposvxx.f.

497*
498* -- LAPACK driver routine --
499* -- LAPACK is a software package provided by Univ. of Tennessee, --
500* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
501*
502* .. Scalar Arguments ..
503 CHARACTER EQUED, FACT, UPLO
504 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
505 \$ N_ERR_BNDS
506 REAL RCOND, RPVGRW
507* ..
508* .. Array Arguments ..
509 INTEGER IWORK( * )
510 REAL A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
511 \$ X( LDX, * ), WORK( * )
512 REAL S( * ), PARAMS( * ), BERR( * ),
513 \$ ERR_BNDS_NORM( NRHS, * ),
514 \$ ERR_BNDS_COMP( NRHS, * )
515* ..
516*
517* ==================================================================
518*
519* .. Parameters ..
520 REAL ZERO, ONE
521 parameter( zero = 0.0e+0, one = 1.0e+0 )
522 INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
523 INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
524 INTEGER CMP_ERR_I, PIV_GROWTH_I
525 parameter( final_nrm_err_i = 1, final_cmp_err_i = 2,
526 \$ berr_i = 3 )
527 parameter( rcond_i = 4, nrm_rcond_i = 5, nrm_err_i = 6 )
528 parameter( cmp_rcond_i = 7, cmp_err_i = 8,
529 \$ piv_growth_i = 9 )
530* ..
531* .. Local Scalars ..
532 LOGICAL EQUIL, NOFACT, RCEQU
533 INTEGER INFEQU, J
534 REAL AMAX, BIGNUM, SMIN, SMAX,
535 \$ SCOND, SMLNUM
536* ..
537* .. External Functions ..
538 EXTERNAL lsame, slamch, sla_porpvgrw
539 LOGICAL LSAME
540 REAL SLAMCH, SLA_PORPVGRW
541* ..
542* .. External Subroutines ..
543 EXTERNAL spoequb, spotrf, spotrs, slacpy, slaqsy,
545* ..
546* .. Intrinsic Functions ..
547 INTRINSIC max, min
548* ..
549* .. Executable Statements ..
550*
551 info = 0
552 nofact = lsame( fact, 'N' )
553 equil = lsame( fact, 'E' )
554 smlnum = slamch( 'Safe minimum' )
555 bignum = one / smlnum
556 IF( nofact .OR. equil ) THEN
557 equed = 'N'
558 rcequ = .false.
559 ELSE
560 rcequ = lsame( equed, 'Y' )
561 ENDIF
562*
563* Default is failure. If an input parameter is wrong or
564* factorization fails, make everything look horrible. Only the
565* pivot growth is set here, the rest is initialized in SPORFSX.
566*
567 rpvgrw = zero
568*
569* Test the input parameters. PARAMS is not tested until SPORFSX.
570*
571 IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.
572 \$ lsame( fact, 'F' ) ) THEN
573 info = -1
574 ELSE IF( .NOT.lsame( uplo, 'U' ) .AND.
575 \$ .NOT.lsame( uplo, 'L' ) ) THEN
576 info = -2
577 ELSE IF( n.LT.0 ) THEN
578 info = -3
579 ELSE IF( nrhs.LT.0 ) THEN
580 info = -4
581 ELSE IF( lda.LT.max( 1, n ) ) THEN
582 info = -6
583 ELSE IF( ldaf.LT.max( 1, n ) ) THEN
584 info = -8
585 ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
586 \$ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
587 info = -9
588 ELSE
589 IF ( rcequ ) THEN
590 smin = bignum
591 smax = zero
592 DO 10 j = 1, n
593 smin = min( smin, s( j ) )
594 smax = max( smax, s( j ) )
595 10 CONTINUE
596 IF( smin.LE.zero ) THEN
597 info = -10
598 ELSE IF( n.GT.0 ) THEN
599 scond = max( smin, smlnum ) / min( smax, bignum )
600 ELSE
601 scond = one
602 END IF
603 END IF
604 IF( info.EQ.0 ) THEN
605 IF( ldb.LT.max( 1, n ) ) THEN
606 info = -12
607 ELSE IF( ldx.LT.max( 1, n ) ) THEN
608 info = -14
609 END IF
610 END IF
611 END IF
612*
613 IF( info.NE.0 ) THEN
614 CALL xerbla( 'SPOSVXX', -info )
615 RETURN
616 END IF
617*
618 IF( equil ) THEN
619*
620* Compute row and column scalings to equilibrate the matrix A.
621*
622 CALL spoequb( n, a, lda, s, scond, amax, infequ )
623 IF( infequ.EQ.0 ) THEN
624*
625* Equilibrate the matrix.
626*
627 CALL slaqsy( uplo, n, a, lda, s, scond, amax, equed )
628 rcequ = lsame( equed, 'Y' )
629 END IF
630 END IF
631*
632* Scale the right-hand side.
633*
634 IF( rcequ ) CALL slascl2( n, nrhs, s, b, ldb )
635*
636 IF( nofact .OR. equil ) THEN
637*
638* Compute the Cholesky factorization of A.
639*
640 CALL slacpy( uplo, n, n, a, lda, af, ldaf )
641 CALL spotrf( uplo, n, af, ldaf, info )
642*
643* Return if INFO is non-zero.
644*
645 IF( info.NE.0 ) THEN
646*
647* Pivot in column INFO is exactly 0
648* Compute the reciprocal pivot growth factor of the
649* leading rank-deficient INFO columns of A.
650*
651 rpvgrw = sla_porpvgrw( uplo, info, a, lda, af, ldaf, work )
652 RETURN
653 ENDIF
654 END IF
655*
656* Compute the reciprocal growth factor RPVGRW.
657*
658 rpvgrw = sla_porpvgrw( uplo, n, a, lda, af, ldaf, work )
659*
660* Compute the solution matrix X.
661*
662 CALL slacpy( 'Full', n, nrhs, b, ldb, x, ldx )
663 CALL spotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
664*
665* Use iterative refinement to improve the computed solution and
666* compute error bounds and backward error estimates for it.
667*
668 CALL sporfsx( uplo, equed, n, nrhs, a, lda, af, ldaf,
669 \$ s, b, ldb, x, ldx, rcond, berr, n_err_bnds, err_bnds_norm,
670 \$ err_bnds_comp, nparams, params, work, iwork, info )
671
672*
673* Scale solutions.
674*
675 IF ( rcequ ) THEN
676 CALL slascl2 ( n, nrhs, s, x, ldx )
677 END IF
678*
679 RETURN
680*
681* End of SPOSVXX
682*
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 slascl2(M, N, D, X, LDX)
SLASCL2 performs diagonal scaling on a matrix.
Definition: slascl2.f:90
subroutine sporfsx(UPLO, EQUED, N, NRHS, A, LDA, AF, LDAF, S, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS, ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS, WORK, IWORK, INFO)
SPORFSX
Definition: sporfsx.f:394
subroutine spotrf(UPLO, N, A, LDA, INFO)
SPOTRF
Definition: spotrf.f:107
real function sla_porpvgrw(UPLO, NCOLS, A, LDA, AF, LDAF, WORK)
SLA_PORPVGRW computes the reciprocal pivot growth factor norm(A)/norm(U) for a symmetric or Hermitian...
Definition: sla_porpvgrw.f:104
subroutine spotrs(UPLO, N, NRHS, A, LDA, B, LDB, INFO)
SPOTRS
Definition: spotrs.f:110
subroutine spoequb(N, A, LDA, S, SCOND, AMAX, INFO)
SPOEQUB
Definition: spoequb.f:118
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: