DPOSVX  - 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 DPOSVX( FACT,  UPLO, N, NRHS, A, LDA, AF, LDAF,
                          EQUED, S, B, LDB, X, LDX, RCOND,  FERR,
                          BERR, WORK, IWORK, INFO )

           CHARACTER      EQUED, FACT, UPLO

           INTEGER        INFO, LDA, LDAF, LDB, LDX, N, NRHS

           DOUBLE         PRECISION RCOND

           INTEGER        IWORK( * )

           DOUBLE         PRECISION  A(  LDA, * ), AF( LDAF, * ),
                          B( LDB, * ), BERR( * ), FERR( * ), S( *
                          ), WORK( * ), X( LDX, * )


PURPOSE

       DPOSVX  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 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


       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, AF contains the fac­
               tored 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 nec­
               essary, then copied to AF 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.

       A       (input/output) DOUBLE PRECISION array, dimension
               (LDA,N)
               On entry, the symmetric matrix A, except if FACT =
               'F' and EQUED = 'Y', then A must contain the equi­
               librated matrix diag(S)*A*diag(S).  If UPLO = 'U',
               the leading N-by-N upper triangular part of A con­
               tains  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 tri­
               angular  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 over­
               written by diag(S)*A*diag(S).

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

       AF      (input  or  output)  DOUBLE  PRECISION  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).

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

       EQUED   (input or output) CHARACTER*1
               Specifies the form of equilibration that was done.
               =  'N':   No  equilibration (always true if FACT =

               = '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) DOUBLE PRECISION 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.

       B         (input/output)   DOUBLE   PRECISION   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.

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

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

               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) DOUBLE PRECISION 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
               <= 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  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.


Man(1) output converted with man2html