Purpose
 =======
 
      LA_SYSVX computes the solution to a linear system of equations 
 A*X = B, where A is a real or complex symmetric matrix and X and B are
 rectangular matrices or vectors.
      LA_HESVX computes the solution to a linear system of equations 
 A*X = B, where A is a complex Hermitian matrix and X and B are
 rectangular matrices or vectors.
      LA_SYSVX and LA_HESVX can also optionally estimate the condition 
 number of A and compute error bounds.
 
 =========
 
         SUBROUTINE LA_SYSVX / LA HESVX( A, B, X, UPLO=uplo, AF=af, &
                        IPIV=ipiv, FACT=fact, FERR=ferr, BERR=berr, &
                        RCOND=rcond, INFO=info )
              <type>(<wp>), INTENT(IN) :: A(:,:), <rhs>
              <type>(<wp>), INTENT(OUT) :: <sol>
              CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO
              <type>(<wp>), INTENT(INOUT), OPTIONAL :: AF(:,:)
              INTEGER, INTENT(INOUT), OPTIONAL :: IPIV(:)
              CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: FACT
              REAL(<wp>), INTENT(OUT), OPTIONAL :: <err>, RCOND
              INTEGER, INTENT(OUT), OPTIONAL :: INFO
         where
              <type> ::= REAL | COMPLEX
              <wp>   ::= KIND(1.0) | KIND(1.0D0)
              <rhs>  ::= B(:,:) | B(:)
              <sol>  ::= X(:,:) | X(:)
              <err>  ::= FERR(:), BERR(:) | FERR, BERR
 
 Arguments
 =========
 
 A       (input) REAL or COMPLEX square array, shape (:,:).
         The symmetric or Hermitian matrix A.
         If UPLO = 'U', the 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 
 	 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.
 B       (input) REAL or COMPLEX array, shape (:,:) with size(B,1) = 
         size(A,1) or shape (:) with size(B) = size(A,1).
         The matrix B.
 X       (output) REAL or COMPLEX array, shape (:,:) with size(X,1) = 
         size(A,1) and size(X,2) = size(B,2), or shape (:) with size(X)
 	  = size(A,1).
         The solution matrix X.
 UPLO    Optional (input) CHARACTER(LEN=1).
            = 'U': Upper triangle of A is stored;
            = 'L': Lower triangle of A is stored.
         Default value: 'U'.
 AF      Optional (input or output) REAL or COMPLEX array, shape (:,:)
         with the same size as A.
         If FACT = 'F', then AF is an input argument that contains the 
 	 block diagonal matrix D and the multipliers used to obtain the
 	 factor L or U from the factorization of A, returned by a
         previous call to LA_SYSVX or LA_HESVX.
         If FACT = 'N', then AF is an output argument that contains the
 	 block diagonal matrix D and the multipliers used to obtain the
 	 factor L or U from the factorization of A.
 IPIV    Optional (input or output) INTEGER array, shape (:) with 
         size(IPIV) = size(A,1).
         If FACT = 'F', then IPIV is an input argument that contains 
 	 details of the row and column interchanges and the block 
 	 structure of D.
         If IPIV(k) > 0 , then rows and columns k and IPIV(k) were 
 	 interchanged and D(k,k) is a 1 by 1 diagonal block.
         If IPIV(k) < 0 , then there are two cases:
           1. If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
 	     columns k-1 and -IPIV(k) were interchanged and 
 	     D(k-1:k,k-1:k) is a 2 by 2 diagonal block.
           2. If UPLO = 'L' and IPIV(k) = IPIV(k+1) < 0, then rows and
 	     columns k+1 and -IPIV(k) were interchanged and 
 	     D(k:k+1,k:k+1) is a 2 by 2 diagonal block.
         If FACT = 'N', then IPIV is an output argument that contains 
 	 details of the row and column interchanges and the block 
 	 structure of D; as described above.
 FACT    Optional (input) CHARACTER(LEN=1).
         Specifies whether the factored form of the matrix A has been 
 	 supplied on entry.
           = 'N': The matrix A will be copied to AF and factored.
           = 'F': AF and IPIV contain the factored form of A.
         Default value: 'N'.
 FERR    Optional (output) REAL array of shape (:), with 
         size(FERR) = size(X,2), or REAL scalar.
         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    Optional (output) REAL array of shape (:), with size(BERR) = 
         size(X,2), or REAL scalar.
         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).
 RCOND   Optional (output) REAL
         The estimate of the reciprocal condition number of A. If RCOND 
 	 is less than the machine
         precision, the matrix is singular to working precision. This 
 	 condition is indicated by a return code of INFO > 0.
 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: D(i,i) = 0. The factorization has been completed, but
 	         the block diagonal matrix D is singular, so the 
 		 solution could not be computed.
             = n+1: D is nonsingular, but RCOND is less than machine 
 	         precision, so the matrix is singular to working 
 		 precision. Nevertheless, the solution and error bounds
 		 are computed because the computed solution can be more
 		 accurate than the value of RCOND would suggest.
             n is the order of A.
         If INFO is not present and an error occurs, then the program is
         terminated with an error message.