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

NAME
  DPPSVX - 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 DPPSVX( 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

      DOUBLE	     PRECISION RCOND

      INTEGER	     IWORK( * )

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

PURPOSE
  DPPSVX 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
  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 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, AFP contains	the
	  factored form	of A.  If EQUED	= 'Y', the matrix A has	been equili-
	  brated 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 necessary, 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.

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

  AP	  (input/output) DOUBLE	PRECISION array, dimension (N*(N+1)/2)
	  On entry, the	upper or lower triangle	of the symmetric 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 overwritten by
	  diag(S)*A*diag(S).

  AFP	  (input or output) DOUBLE PRECISION array, dimension
	  (N*(N+1)/2) If FACT =	'F', then AFP is an input argument and on
	  entry	contains the triangular	factor U or L from the Cholesky	fac-
	  torization 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 equili-
	  brated 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/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.

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 ]