Purpose
 =======
 
     LA_PBSVX computes the solution to a linear system of equations
 A*X = B, where A has band form and is real symmetric or complex 
 Hermitian and, in either case, positive definite, and where X and B 
 are rectangular matrices or vectors.
     LA_PBSVX can also optionally equilibrate the system if A is poorly
 scaled, estimate the condition number of (the equilibrated) A, and 
 compute error bounds.
 
 =========
 
       SUBROUTINE LA_PBSVX( AB, B, X, UPLO=uplo, AFB=afb, FACT=fact, &
                             EQUED=equed, S=s, FERR=ferr, BERR=berr, &
                             RCOND=rcond, INFO=info )
             <type>(<wp>), INTENT(INOUT) :: AB(:,:), <rhs>
             <type>(<wp>), INTENT(OUT) :: <sol>
             CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: UPLO
             <type>(<wp>), INTENT(INOUT), OPTIONAL :: AFB(:,:)
             CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: FACT
             CHARACTER(LEN=1), INTENT(INOUT), OPTIONAL :: EQUED
             REAL(<wp>), INTENT(INOUT), OPTIONAL :: S(:)
             REAL(<wp>), INTENT(OUT), OPTIONAL :: <err>
             REAL(<wp>), INTENT(OUT), OPTIONAL :: 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
 =========
 
 AB       (input/output) REAL or COMPLEX array, shape (:,:) with 
          size(AB,1) = kd + 1 and size(AB,2) = n, where kd is the number
          of superdiagonals or subdiagonals and n is the order of A.
          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') 
          triangle of matrix A, or its equilibration, in band storage. 
          The (kd + 1) diagonals of A are stored in the rows of AB so
          that the j-th column of A is stored in the j-th column of AB
   	  as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j
                                                   1<=j<=n
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd)
                                                   1<=j<=n.
          On exit, if FACT = 'E', then the equilibrated version of A is
          stored in AB; otherwise, AB is unchanged.
 B        (input/output) REAL or COMPLEX array, shape (:,:) with 
          size(B,1) = n or shape (:) with size(B) = n.
          On entry, the matrix B.
          On exit, the scaled version of B if the system has been 
          equilibrated; otherwise, B is unchanged.
 X        (output) REAL or COMPLEX array, shape (:,:) with size(X,1)=n
          and size(X,2) = size(B,2), or shape (:) with size(X) = n.
          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'.
 AFB      Optional (input or output) REAL or COMPLEX array, shape (:) 
          with the same size as AB.
          If FACT = 'F' then AFB is an input argument that contains the
          factor U or L from the Cholesky factorization of (the
 	  equilibrated) A, in the same storage format as A, returned by 
          a previous call to LA_PBSVX.
          If FACT /= 'F' then AFB is an output argument that contains 
 	  the factor U or L from the Cholesky factorization of (the 
          equilibrated) A in the same storage format as A.
 FACT     Optional (input) CHARACTER(LEN=1).
          Specifies whether the factored form of the matrix A is supplied
          on entry, and if not, whether A should be equilibrated before 
          it is factored.
            = 'N': The matrix A will be copied to AFB and factored (no 
 	          equilibration).
            = 'E': The matrix A will be equilibrated, then copied to 
 	          AFB and factored.
            = 'F': AFB contains the factored form of (the equilibrated)
 	          A.
          Default value: 'N'.
 EQUED    Optional (input or output) CHARACTER(LEN=1).
          Specifies the form of equilibration that was done.
          EQUED is an input argument if FACT = 'F', otherwise it is an
          output argument
            = 'N': No equilibration (always true if FACT = 'N').
            = 'Y': Equilibration, i.e., A has been premultiplied and 
 	          postmultiplied by diag(S).
          Default value: 'N'.
 S        Optional (input or output) REAL array, shape (:) with size(S) 
          = size(A,1).
          The scaling factors for A.
          S is an input argument if FACT = 'F' and EQUED = 'Y'.
          S is an output argument if FACT = 'E' and EQUED = 'Y'.
 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 (the 
          equilibrated) 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     Optional (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 (the equilibrated) A 
 	         is not positive definite, so the factorization could 
 		 not be completed and the solution and error bounds 
 		 could not be computed. RCOND= 0 is returned.
             = n+1: U or L 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.
          If INFO is not present and an error occurs, then the program 
          is terminated with an error message.