CXML
        	    
LAPACK version 3.0
zposvx(3)

PURPOSE
  ZPOSVX - use the Cholesky factorization A = U**H*U or A = L*L**H to compute
  the solution to a complex system of linear equations A * X = B,

SYNTAX
  SUBROUTINE ZPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B, LDB,
		     X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )

      CHARACTER	     EQUED, FACT, UPLO

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

      DOUBLE	     PRECISION RCOND

      DOUBLE	     PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * )

      COMPLEX*16     A( LDA, * ), AF( LDAF, * ), B( LDB, * ), WORK( * ), X(
		     LDX, * )

DESCRIPTION
  ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to compute
  the solution to a complex system of linear equations A * X = B, where A is
  an N-by-N Hermitian 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**H* U,  if UPLO = 'U', or
	A = L * L**H,  if UPLO = 'L',
     where U is an upper triangular matrix and L is a lower triangular
     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 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
	  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.

  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.

  A	  (input/output) COMPLEX*16 array, dimension (LDA,N)
	  On entry, the Hermitian 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).

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

  AF	  (input or output) COMPLEX*16 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**H*U or A = L*L**H, 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**H*U or A = L*L**H 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**H*U or A = L*L**H 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 = '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) DOUBLE PRECISION 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.

  B	  (input/output) COMPLEX*16 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) COMPLEX*16 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 particular, 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 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.

  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 element of A or B
	  that makes X(j) an exact solution).

  WORK	  (workspace) COMPLEX*16 array, dimension (2*N)

  RWORK	  (workspace) DOUBLE PRECISION 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:	 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 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.