LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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

Download SPOSVX + dependencies [TGZ] [ZIP] [TXT]

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 principal minor of order i is not positive,
    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 principal minor of order i of A
                       is not positive, 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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 xerbla(srname, info)
Definition cblat2.f:3285
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine spocon(uplo, n, a, lda, anorm, rcond, work, iwork, info)
SPOCON
Definition spocon.f:121
subroutine spoequ(n, a, lda, s, scond, amax, info)
SPOEQU
Definition spoequ.f:112
subroutine sporfs(uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx, ferr, berr, work, iwork, info)
SPORFS
Definition sporfs.f:183
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:107
subroutine spotrs(uplo, n, nrhs, a, lda, b, ldb, info)
SPOTRS
Definition spotrs.f:110
Here is the call graph for this function:
Here is the caller graph for this function: