DPOSVX(l)	     LAPACK driver routine (version 1.1)	    DPOSVX(l)

NAME
  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 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. 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, steps	4-6 are	skipped.

  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.

ARGUMENTS

  FACT	  (input) CHARACTER*1
	  Specifies whether or not the factored	form of	the matrix A is	sup-
	  plied	on entry, and if not, whether the matrix A should be equili-
	  brated before	it is factored.	 = 'F':	 On entry, AF contains the
	  factored form	of A.  If EQUED	= 'Y', the matrix A has	been equili-
	  brated with scaling factors given by S.  A and AF will not be	modi-
	  fied.	 = '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.

  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.

  NRHS	  (input) INTEGER
	  The number of	righthand 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 equilibrated 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 overwritten 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 descrip-
	  tion 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/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 variable if	FACT = 'F';
	  otherwise, it	is an output variable.

  S	  (input/output) DOUBLE	PRECISION array, dimension (N)
	  The scale factors for	A; not accessed	if EQUED = 'N'.	 S is an
	  input	variable if FACT = 'F';	otherwise, S is	an output variable.
	  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 righthand 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, 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	particular, if RCOND = 0), the matrix is singular to
	  working precision.  This condition is	indicated by a return code of
	  INFO > 0, and	the solution and error bounds are not computed.

  FERR	  (output) DOUBLE PRECISION array, dimension (NRHS)
	  The estimated	forward	error bounds for each solution vector X(j)
	  (the j-th column of the solution matrix X).  If XTRUE	is the true
	  solution, FERR(j) bounds the magnitude of the	largest	entry in
	  (X(j)	- XTRUE) divided by the	magnitude of the largest entry in
	  X(j).	 The quality of	the error bound	depends	on the quality of the
	  estimate of norm(inv(A)) computed in the code; if the	estimate of
	  norm(inv(A)) is accurate, the	error bound is guaranteed.

  BERR	  (output) DOUBLE PRECISION array, dimension (NRHS)
	  The componentwise relative backward error of each solution vector
	  X(j) (i.e., the smallest relative change in any entry	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 illegal value
	  > 0: if INFO = i, and	i is
	  <= N:	if INFO	= i, the leading minor of order	i of A is not posi-
	  tive definite, so the	factorization could not	be completed, and the
	  solution and error bounds could not be computed.  = N+1: RCOND is
	  less than machine precision.	The factorization has been completed,
	  but the matrix is singular to	working	precision, and the solution
	  and error bounds have	not been computed.