MODULE CPTSVX_MODULE
   INTERFACE CPTSVX
      MODULE PROCEDURE CPTSVX_ORIGINAL
      MODULE PROCEDURE CPTSVX_EXTENDED
   END INTERFACE
CONTAINS
      SUBROUTINE CPTSVX_ORIGINAL( FACT, N, NRHS, D, E, DF, EF, B, LDB,  &
     &                   X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
!
!  -- LAPACK routine (version 3.0) --
!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
!     Courant Institute, Argonne National Lab, and Rice University
!     June 30, 1999
!
!     .. Scalar Arguments ..
      CHARACTER          FACT
      INTEGER            INFO, LDB, LDX, N, NRHS
      REAL               RCOND
!     ..
!     .. Array Arguments ..
      REAL               BERR( * ), D( * ), DF( * ), FERR( * ),         &
     &                   RWORK( * )
      COMPLEX            B( LDB, * ), E( * ), EF( * ), WORK( * ),       &
     &                   X( LDX, * )
!     ..
!
!  Purpose
!  =======
!
!  CPTSVX uses the factorization A = L*D*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 tridiagonal 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 = 'N', the matrix A is factored as A = L*D*L**H, where L
!     is a unit lower bidiagonal matrix and D is diagonal.  The
!     factorization can also be regarded as having the form
!     A = U**H*D*U.
!
!  2. 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.
!
!  3. The system of equations is solved for X using the factored form
!     of A.
!
!  4. Iterative refinement is applied to improve the computed solution
!     matrix and calculate error bounds and backward error estimates
!     for it.
!
!  Arguments
!  =========
!
!  FACT    (input) CHARACTER*1
!          Specifies whether or not the factored form of the matrix
!          A is supplied on entry.
!          = 'F':  On entry, DF and EF contain the factored form of A.
!                  D, E, DF, and EF will not be modified.
!          = 'N':  The matrix A will be copied to DF and EF and
!                  factored.
!
!  N       (input) INTEGER
!          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.
!
!  D       (input) REAL array, dimension (N)
!          The n diagonal elements of the tridiagonal matrix A.
!
!  E       (input) COMPLEX array, dimension (N-1)
!          The (n-1) subdiagonal elements of the tridiagonal matrix A.
!
!  DF      (input or output) REAL array, dimension (N)
!          If FACT = 'F', then DF is an input argument and on entry
!          contains the n diagonal elements of the diagonal matrix D
!          from the L*D*L**H factorization of A.
!          If FACT = 'N', then DF is an output argument and on exit
!          contains the n diagonal elements of the diagonal matrix D
!          from the L*D*L**H factorization of A.
!
!  EF      (input or output) COMPLEX array, dimension (N-1)
!          If FACT = 'F', then EF is an input argument and on entry
!          contains the (n-1) subdiagonal elements of the unit
!          bidiagonal factor L from the L*D*L**H factorization of A.
!          If FACT = 'N', then EF is an output argument and on exit
!          contains the (n-1) subdiagonal elements of the unit
!          bidiagonal factor L from the L*D*L**H factorization of A.
!
!  B       (input) COMPLEX array, dimension (LDB,NRHS)
!          The N-by-NRHS right hand side matrix B.
!
!  LDB     (input) INTEGER
!          The leading dimension of the array B.  LDB >= max(1,N).
!
!  X       (output) COMPLEX array, dimension (LDX,NRHS)
!          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
!
!  LDX     (input) INTEGER
!          The leading dimension of the array X.  LDX >= max(1,N).
!
!  RCOND   (output) REAL
!          The reciprocal condition number of the matrix A.  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) REAL array, dimension (NRHS)
!          The 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).
!
!  BERR    (output) REAL 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 array, dimension (N)
!
!  RWORK   (workspace) REAL 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.
!
!  =====================================================================
!
!     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
!     ..
!     .. Local Scalars ..
      LOGICAL            NOFACT
      REAL               ANORM
!     ..
!     .. External Functions ..
      LOGICAL            LSAME
      REAL               CLANHT, SLAMCH
      EXTERNAL           LSAME, CLANHT, SLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL           CCOPY, CLACPY, CPTCON, CPTRFS, CPTTRF, CPTTRS, &
     &                   SCOPY, XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC          MAX
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      INFO = 0
      NOFACT = LSAME( FACT, 'N' )
      IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -9
      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
         INFO = -11
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CPTSVX', -INFO )
         RETURN
      END IF
!
      IF( NOFACT ) THEN
!
!        Compute the L*D*L' (or U'*D*U) factorization of A.
!
         CALL SCOPY( N, D, 1, DF, 1 )
         IF( N.GT.1 ) CALL CCOPY( N-1, E, 1, EF, 1 )
         CALL CPTTRF( N, DF, EF, INFO )
!
!        Return if INFO is non-zero.
!
         IF( INFO.NE.0 ) THEN
            IF( INFO.GT.0 ) RCOND = ZERO
            RETURN
         END IF
      END IF
!
!     Compute the norm of the matrix A.
!
      ANORM = CLANHT( '1', N, D, E )
!
!     Compute the reciprocal of the condition number of A.
!
      CALL CPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
!
!     Set INFO = N+1 if the matrix is singular to working precision.
!
      IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) INFO = N + 1
!
!     Compute the solution vectors X.
!
      CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
      CALL CPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO )
!
!     Use iterative refinement to improve the computed solutions and
!     compute error bounds and backward error estimates for them.
!
      CALL CPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,&
     &             BERR, WORK, RWORK, INFO )
!
      RETURN
!
!     End of CPTSVX_ORIGINAL
!
      END SUBROUTINE CPTSVX_ORIGINAL
      SUBROUTINE CPTSVX_EXTENDED( FACT, UPLO, N, NRHS, D, E, DF, EF, B, &
     &           LDB, X, LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
!
!  -- LAPACK routine (version 3.0) --
!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
!     Courant Institute, Argonne National Lab, and Rice University
!     June 30, 1999
!
!     .. Scalar Arguments ..
      CHARACTER          FACT, UPLO
      INTEGER            INFO, LDB, LDX, N, NRHS
      REAL               RCOND
!     ..
!     .. Array Arguments ..
      REAL               BERR( * ), D( * ), DF( * ), FERR( * ),         &
     &                   RWORK( * )
      COMPLEX            B( LDB, * ), E( * ), EF( * ), WORK( * ),       &
     &                   X( LDX, * )
!     ..
!
!  Purpose
!  =======
!
!  CPTSVX_EXTENDED uses the factorization A = L*D*L*' or A = U'*D*U to
!  compute the solution to a complex system of linear equations A*X = B,
!  where A is an N-by-N Hermitian positive definite tridiagonal 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 = 'N' and UPLO = 'L', the matrix A is factored as
!     A = L*D*L', where L is a unit lower bidiagonal matrix and D is
!     diagonal.  If FACT = 'N' and UPLO = 'U', the matrix A is factored
!     as A = U'*D*U, where U is a unit upper bidiagonal matrix.
!
!  2. 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.
!
!  3. The system of equations is solved for X using the factored form
!     of A.
!
!  4. Iterative refinement is applied to improve the computed solution
!     matrix and calculate error bounds and backward error estimates
!     for it.
!
!  Arguments
!  =========
!
!  FACT    (input) CHARACTER*1
!          Specifies whether or not the factored form of the matrix
!          A is supplied on entry.
!          = 'F':  On entry, DF and EF contain the factored form of A.
!                  D, E, DF, and EF will not be modified.
!          = 'N':  The matrix A will be copied to DF and EF and
!                  factored.
!
!  UPLO    (input) CHARACTER*1
!          Specifies whether the subdiagonal or the superdiagonal of the
!          tridiagonal matrix A is stored and the form of the
!          factorization.
!          = 'L':  E is the subdiagonal of A, form is A = L*D*L'
!          = 'U':  E is the superdiagonal of A, form is A = U'*D*U
!
!  N       (input) INTEGER
!          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.
!
!  D       (input) REAL array, dimension (N)
!          The n diagonal elements of the tridiagonal matrix A.
!
!  E       (input) COMPLEX array, dimension (N-1)
!          If UPLO = 'L', the (N-1) subdiagonal elements of the
!          tridiagonal matrix A; if UPLO = 'U', the (N-1) superdiagonal
!          elements of A.
!
!  DF      (input or output) REAL array, dimension (N)
!          If FACT = 'F', then DF is an input argument and on entry
!          contains the n diagonal elements of the diagonal matrix D
!          from the factorization A = L*D*L' or A = U'*D*U.
!          If FACT = 'N', then DF is an output argument and on exit
!          contains the n diagonal elements of the diagonal matrix D
!          from the factorization A = L*D*L' or A = U'*D*U.
!
!  EF      (input or output) COMPLEX array, dimension (N-1)
!          If FACT = 'F', then EF is an input argument and on entry
!          contains, if UPLO = 'L', the (N-1) subdiagonal elements of
!          the unit bidiagonal factor L from the factorization
!          A = L*D*L', or if UPLO = 'U', the (N-1) superdiagonal
!          elements of the unit bidiagonal factor U from the
!          factorization A = U'*D*U.
!
!          If FACT = 'N', then EF is an output argument and on exit
!          contains, if UPLO = 'L', the (N-1) subdiagonal elements of
!          the unit bidiagonal factor L from the factorization
!          A = L*D*L', or if UPLO = 'U', the (N-1) superdiagonal
!          elements of the unit bidiagonal factor U from the
!          factorization A = U'*D*U.
!
!  B       (input) COMPLEX array, dimension (LDB,NRHS)
!          The N-by-NRHS right hand side matrix B.
!
!  LDB     (input) INTEGER
!          The leading dimension of the array B.  LDB >= max(1,N).
!
!  X       (output) COMPLEX array, dimension (LDX,NRHS)
!          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
!
!  LDX     (input) INTEGER
!          The leading dimension of the array X.  LDX >= max(1,N).
!
!  RCOND   (output) REAL
!          The reciprocal condition number of the matrix A.  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) REAL array, dimension (NRHS)
!          The 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).
!
!  BERR    (output) REAL 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 array, dimension (N)
!
!  RWORK   (workspace) REAL 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.
!
!  =====================================================================
!
!     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E+0 )
!     ..
!     .. Local Scalars ..
      LOGICAL            NOFACT
      REAL               ANORM
!     ..
!     .. External Functions ..
      LOGICAL            LSAME
      REAL               CLANHT, SLAMCH
      EXTERNAL           LSAME, CLANHT, SLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL           CCOPY, CLACPY, CPTCON, CPTRFS, CPTTRF, CPTTRS, &
     &                   SCOPY, XERBLA
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC          MAX
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      INFO = 0
      NOFACT = LSAME( FACT, 'N' )
      IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )  &
     &   THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -10
      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
         INFO = -12
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CPTSVX', -INFO )
         RETURN
      END IF
!
      IF( NOFACT ) THEN
!
!        Compute the L*D*L' (or U'*D*U) factorization of A.
!
         CALL SCOPY( N, D, 1, DF, 1 )
         IF( N.GT.1 ) CALL CCOPY( N-1, E, 1, EF, 1 )
         CALL CPTTRF( N, DF, EF, INFO )
!
!        Return if INFO is non-zero.
!
         IF( INFO.NE.0 ) THEN
            IF( INFO.GT.0 ) RCOND = ZERO
            RETURN
         END IF
      END IF
!
!     Compute the norm of the matrix A.
!
      ANORM = CLANHT( '1', N, D, E )
!
!     Compute the reciprocal of the condition number of A.
!
      CALL CPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
!
!     Set INFO = N+1 if the matrix is singular to working precision.
!
      IF( RCOND.LT.SLAMCH( 'Epsilon' ) ) INFO = N + 1
!
!     Compute the solution vectors X.
!
      CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
      CALL CPTTRS( UPLO, N, NRHS, DF, EF, X, LDX, INFO )
!
!     Use iterative refinement to improve the computed solutions and
!     compute error bounds and backward error estimates for them.
!
      CALL CPTRFS( UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,   &
     &             BERR, WORK, RWORK, INFO )
!
      RETURN
!
!     End of CPTSVX_EXTENDED
!
      END SUBROUTINE CPTSVX_EXTENDED
END MODULE CPTSVX_MODULE