SPPSVX  - use 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,


SYNOPSIS

       SUBROUTINE SPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S,
                          B, LDB,  X,  LDX,  RCOND,  FERR,  BERR,
                          WORK, IWORK, INFO )

           CHARACTER      EQUED, FACT, UPLO

           INTEGER        INFO, LDB, LDX, N, NRHS

           REAL           RCOND

           INTEGER        IWORK( * )

           REAL           AFP( * ), AP( * ), B( LDB, * ), BERR( *
                          ), FERR( * ), S( * ),  WORK(  *  ),  X(
                          LDX, * )


PURPOSE

       SPPSVX  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 posi­
       tive definite matrix stored in packed format 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

          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 fac­
       tored 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  premulti­
       plied by
          diag(S) so that it solves the original system before
          equilibration.


ARGUMENTS

       FACT    (input) 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, AFP contains the fac­
               tored form of A.  If EQUED = 'Y', the matrix A has
               been equilibrated with scaling factors given by S.
               AP  and  AFP  will  not  be modified.  = 'N':  The
               matrix A will be copied to AFP and factored.
               = 'E':  The matrix A will be equilibrated if  nec­
               essary, then copied to AFP and factored.

       UPLO    (input) CHARACTER*1
               = 'U':  Upper triangle of A is stored;
               = 'L':  Lower triangle of A is stored.

       N       (input) INTEGER
               The number of linear equations, i.e., the order of
               the matrix A.  N >= 0.

               The number of right hand sides, i.e.,  the  number
               of columns of the matrices B and X.  NRHS >= 0.

       AP      (input/output) REAL array, dimension (N*(N+1)/2)
               On  entry, the upper or lower triangle of the sym­
               metric matrix A, packed  columnwise  in  a  linear
               array,  except if FACT = 'F' and EQUED = 'Y', then
               A   must   contain   the    equilibrated    matrix
               diag(S)*A*diag(S).  The j-th column of A is stored
               in the array AP as follows: if UPLO = 'U', AP(i  +
               (j-1)*j/2)  =  A(i,j)  for 1<=i<=j; if UPLO = 'L',
               AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.   See
               below  for  further details.  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 over­
               written by diag(S)*A*diag(S).

       AFP     (input or output) REAL array, dimension
               (N*(N+1)/2) If FACT = 'F', then AFP  is  an  input
               argument and on entry contains the triangular fac­
               tor U or L from the  Cholesky  factorization  A  =
               U'*U or A = L*L', in the same storage format as A.
               If EQUED .ne. 'N', then AFP is the  factored  form
               of the equilibrated matrix A.

               If  FACT = 'N', then AFP is an output argument and
               on exit returns the triangular factor U or L  from
               the Cholesky factorization A = U'*U or A = L*L' of
               the original matrix A.

               If FACT = 'E', then AFP is an output argument  and
               on  exit returns the triangular factor U or L from
               the Cholesky factorization A = U'*U or A = L*L' of
               the  equilibrated matrix A (see the description of
               AP for the form of the equilibrated matrix).

       EQUED   (input or output) 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.

       S       (input or output) REAL array, dimension (N)
               The  scale  factors for A; not accessed if EQUED =
               'N'.  S is an input argument if FACT = 'F'; other­
               wise,  S is an output argument.  If FACT = 'F' and
               EQUED = 'Y', each element of S must be positive.

               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.

       LDB     (input) INTEGER
               The leading dimension of  the  array  B.   LDB  >=
               max(1,N).

       X       (output) 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.

       LDX     (input) INTEGER
               The  leading  dimension  of  the  array X.  LDX >=
               max(1,N).

       RCOND   (output) 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  par­
               ticular,  if RCOND = 0), the matrix is singular to
               working precision.  This condition is indicated by
               a return code of INFO > 0.

       FERR    (output) REAL array, dimension (NRHS)
               The  estimated  forward error bound for each solu­
               tion vector X(j) (the j-th column of the  solution
               matrix  X).   If XTRUE is the true solution corre­
               sponding 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 reli­
               able as the estimate  for  RCOND,  and  is  almost
               always a slight overestimate of the true error.

       BERR    (output) 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).

       WORK    (workspace) REAL array, dimension (3*N)

       IWORK   (workspace) INTEGER array, dimension (N)

       INFO    (output) INTEGER
               = 0:  successful exit
               < 0:  if INFO = -i, the i-th argument had an ille­
               gal value
               > 0:  if INFO = i, and i is

               positive definite, so the factorization could  not
               be  completed,  and the solution has not been com­
               puted. RCOND = 0 is returned.  = N+1: U is nonsin­
               gular,  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.


FURTHER DETAILS

       The packed storage scheme is illustrated by the  following
       example when N = 4, UPLO = 'U':

       Two-dimensional storage of the symmetric matrix A:

          a11 a12 a13 a14
              a22 a23 a24
                  a33 a34     (aij = conjg(aji))
                      a44

       Packed storage of the upper triangle of A:

       AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]


Man(1) output converted with man2html