Implement packed triangular solver.

8 files changed, 144 insertions, 1264 deletions

diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt
index 171b75aa1..3877e1285 100644
--- a/blas/CMakeLists.txt
+++ b/blas/CMakeLists.txt
@@ -18,10 +18,10 @@ if(EIGEN_Fortran_COMPILER_WORKS)
 set(EigenBlas_SRCS ${EigenBlas_SRCS}
     complexdots.f
     srotm.f srotmg.f drotm.f drotmg.f
-    lsame.f  dspmv.f  dtpsv.f ssbmv.f
-    chbmv.f  chpr.f   sspmv.f  stpsv.f
-    zhbmv.f  zhpr.f   chpmv.f  ctpsv.f  dsbmv.f
-    zhpmv.f  ztpsv.f
+    lsame.f  dspmv.f  ssbmv.f
+    chbmv.f  chpr.f   sspmv.f
+    zhbmv.f  zhpr.f   chpmv.f dsbmv.f
+    zhpmv.f
     dtbmv.f stbmv.f ctbmv.f ztbmv.f
 )
 else()
diff --git a/blas/PackedTriangularSolverVector.h b/blas/PackedTriangularSolverVector.h
new file mode 100644
index 000000000..5c0bb4bd6
--- /dev/null
+++ b/blas/PackedTriangularSolverVector.h
@@ -0,0 +1,88 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
+#define EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
+
+namespace internal {
+
+template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder>
+struct packed_triangular_solve_vector;
+
+// forward and backward substitution, row-major, rhs is a vector
+template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
+struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, RowMajor>
+{
+  enum {
+    IsLower = (Mode&Lower)==Lower
+  };
+  static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
+  {
+    internal::conj_if<Conjugate> cj;
+    typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap;
+    typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType;
+
+    lhs += IsLower ? 0 : (size*(size+1)>>1)-1;
+    for(Index pi=0; pi<size; ++pi)
+    {
+      Index i = IsLower ? pi : size-pi-1;
+      Index s = IsLower ? 0 : 1;
+      if (pi>0)
+	rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi))
+	    .cwiseProduct(Map<const Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower ? 0 : i+1),pi))).sum();
+      if (!(Mode & UnitDiag))
+	rhs[i] /= cj(lhs[IsLower ? i : 0]);
+      IsLower ? lhs += pi+1 : lhs -= pi+2;
+    }
+  }
+};
+
+// forward and backward substitution, column-major, rhs is a vector
+template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate>
+struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Conjugate, ColMajor>
+{
+  enum {
+    IsLower = (Mode&Lower)==Lower
+  };
+  static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
+  {
+    internal::conj_if<Conjugate> cj;
+    typedef Map<const Matrix<LhsScalar,Dynamic,1> > LhsMap;
+    typedef typename conj_expr_if<Conjugate,LhsMap>::type ConjLhsType;
+
+    lhs += IsLower ? 0 : size*(size-1)>>1;
+    for(Index pi=0; pi<size; ++pi)
+    {
+      Index i = IsLower ? pi : size-pi-1;
+      Index r = size - pi - 1;
+      if (!(Mode & UnitDiag))
+	rhs[i] /= cj(lhs[IsLower ? 0 : i]);
+      if (r>0)
+	Map<Matrix<RhsScalar,Dynamic,1> >(rhs+(IsLower? i+1 : 0),r) -=
+	    rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r));
+      IsLower ? lhs += size-pi : lhs -= r;
+    }
+  }
+};
+
+template<typename LhsScalar, typename RhsScalar, typename Index, int Mode, bool Conjugate, int StorageOrder>
+struct packed_triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheRight, Mode, Conjugate, StorageOrder>
+{
+  static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs)
+  {
+    packed_triangular_solve_vector<LhsScalar,RhsScalar,Index,OnTheLeft,
+	((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
+	Conjugate,StorageOrder==RowMajor?ColMajor:RowMajor
+      >::run(size, lhs, rhs);
+  }
+};
+
+} // end namespace internal
+
+#endif // EIGEN_PACKED_TRIANGULAR_SOLVER_VECTOR_H
diff --git a/blas/common.h b/blas/common.h
index 1019d8623..ee38b31d1 100644
--- a/blas/common.h
+++ b/blas/common.h
@@ -77,6 +77,7 @@ namespace Eigen {
 #include "GeneralRank1Update.h"
 #include "PackedSelfadjointProduct.h"
 #include "PackedTriangularMatrixVector.h"
+#include "PackedTriangularSolverVector.h"
 #include "Rank2Update.h"
 }
 
diff --git a/blas/ctpsv.f b/blas/ctpsv.f
deleted file mode 100644
index 1804797ea..000000000
--- a/blas/ctpsv.f
+++ /dev/null
@@ -1,332 +0,0 @@
-      SUBROUTINE CTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
-*     .. Scalar Arguments ..
-      INTEGER INCX,N
-      CHARACTER DIAG,TRANS,UPLO
-*     ..
-*     .. Array Arguments ..
-      COMPLEX AP(*),X(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  CTPSV  solves one of the systems of equations
-*
-*     A*x = b,   or   A'*x = b,   or   conjg( A' )*x = b,
-*
-*  where b and x are n element vectors and A is an n by n unit, or
-*  non-unit, upper or lower triangular matrix, supplied in packed form.
-*
-*  No test for singularity or near-singularity is included in this
-*  routine. Such tests must be performed before calling this routine.
-*
-*  Arguments
-*  ==========
-*
-*  UPLO   - CHARACTER*1.
-*           On entry, UPLO specifies whether the matrix is an upper or
-*           lower triangular matrix as follows:
-*
-*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
-*
-*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
-*
-*           Unchanged on exit.
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the equations to be solved as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   A*x = b.
-*
-*              TRANS = 'T' or 't'   A'*x = b.
-*
-*              TRANS = 'C' or 'c'   conjg( A' )*x = b.
-*
-*           Unchanged on exit.
-*
-*  DIAG   - CHARACTER*1.
-*           On entry, DIAG specifies whether or not A is unit
-*           triangular as follows:
-*
-*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
-*
-*              DIAG = 'N' or 'n'   A is not assumed to be unit
-*                                  triangular.
-*
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the order of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  AP     - COMPLEX          array of DIMENSION at least
-*           ( ( n*( n + 1 ) )/2 ).
-*           Before entry with  UPLO = 'U' or 'u', the array AP must
-*           contain the upper triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
-*           respectively, and so on.
-*           Before entry with UPLO = 'L' or 'l', the array AP must
-*           contain the lower triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
-*           respectively, and so on.
-*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
-*           A are not referenced, but are assumed to be unity.
-*           Unchanged on exit.
-*
-*  X      - COMPLEX          array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ).
-*           Before entry, the incremented array X must contain the n
-*           element right-hand side vector b. On exit, X is overwritten
-*           with the solution vector x.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
-*     ..
-*     .. Local Scalars ..
-      COMPLEX TEMP
-      INTEGER I,INFO,IX,J,JX,K,KK,KX
-      LOGICAL NOCONJ,NOUNIT
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC CONJG
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
-          INFO = 1
-      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +         .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 2
-      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
-          INFO = 3
-      ELSE IF (N.LT.0) THEN
-          INFO = 4
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('CTPSV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF (N.EQ.0) RETURN
-*
-      NOCONJ = LSAME(TRANS,'T')
-      NOUNIT = LSAME(DIAG,'N')
-*
-*     Set up the start point in X if the increment is not unity. This
-*     will be  ( N - 1 )*INCX  too small for descending loops.
-*
-      IF (INCX.LE.0) THEN
-          KX = 1 - (N-1)*INCX
-      ELSE IF (INCX.NE.1) THEN
-          KX = 1
-      END IF
-*
-*     Start the operations. In this version the elements of AP are
-*     accessed sequentially with one pass through AP.
-*
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  x := inv( A )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK - 1
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K - 1
-   10                     CONTINUE
-                      END IF
-                      KK = KK - J
-   20             CONTINUE
-              ELSE
-                  JX = KX + (N-1)*INCX
-                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 K = KK - 1,KK - J + 1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   30                     CONTINUE
-                      END IF
-                      JX = JX - INCX
-                      KK = KK - J
-   40             CONTINUE
-              END IF
-          ELSE
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK + 1
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K + 1
-   50                     CONTINUE
-                      END IF
-                      KK = KK + (N-J+1)
-   60             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 K = KK + 1,KK + N - J
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   70                     CONTINUE
-                      END IF
-                      JX = JX + INCX
-                      KK = KK + (N-J+1)
-   80             CONTINUE
-              END IF
-          END IF
-      ELSE
-*
-*        Form  x := inv( A' )*x  or  x := inv( conjg( A' ) )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 110 J = 1,N
-                      TEMP = X(J)
-                      K = KK
-                      IF (NOCONJ) THEN
-                          DO 90 I = 1,J - 1
-                              TEMP = TEMP - AP(K)*X(I)
-                              K = K + 1
-   90                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      ELSE
-                          DO 100 I = 1,J - 1
-                              TEMP = TEMP - CONJG(AP(K))*X(I)
-                              K = K + 1
-  100                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK+J-1))
-                      END IF
-                      X(J) = TEMP
-                      KK = KK + J
-  110             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 140 J = 1,N
-                      TEMP = X(JX)
-                      IX = KX
-                      IF (NOCONJ) THEN
-                          DO 120 K = KK,KK + J - 2
-                              TEMP = TEMP - AP(K)*X(IX)
-                              IX = IX + INCX
-  120                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      ELSE
-                          DO 130 K = KK,KK + J - 2
-                              TEMP = TEMP - CONJG(AP(K))*X(IX)
-                              IX = IX + INCX
-  130                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK+J-1))
-                      END IF
-                      X(JX) = TEMP
-                      JX = JX + INCX
-                      KK = KK + J
-  140             CONTINUE
-              END IF
-          ELSE
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 170 J = N,1,-1
-                      TEMP = X(J)
-                      K = KK
-                      IF (NOCONJ) THEN
-                          DO 150 I = N,J + 1,-1
-                              TEMP = TEMP - AP(K)*X(I)
-                              K = K - 1
-  150                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      ELSE
-                          DO 160 I = N,J + 1,-1
-                              TEMP = TEMP - CONJG(AP(K))*X(I)
-                              K = K - 1
-  160                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK-N+J))
-                      END IF
-                      X(J) = TEMP
-                      KK = KK - (N-J+1)
-  170             CONTINUE
-              ELSE
-                  KX = KX + (N-1)*INCX
-                  JX = KX
-                  DO 200 J = N,1,-1
-                      TEMP = X(JX)
-                      IX = KX
-                      IF (NOCONJ) THEN
-                          DO 180 K = KK,KK - (N- (J+1)),-1
-                              TEMP = TEMP - AP(K)*X(IX)
-                              IX = IX - INCX
-  180                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      ELSE
-                          DO 190 K = KK,KK - (N- (J+1)),-1
-                              TEMP = TEMP - CONJG(AP(K))*X(IX)
-                              IX = IX - INCX
-  190                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/CONJG(AP(KK-N+J))
-                      END IF
-                      X(JX) = TEMP
-                      JX = JX - INCX
-                      KK = KK - (N-J+1)
-  200             CONTINUE
-              END IF
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of CTPSV .
-*
-      END
diff --git a/blas/dtpsv.f b/blas/dtpsv.f
deleted file mode 100644
index c7e58d32f..000000000
--- a/blas/dtpsv.f
+++ /dev/null
@@ -1,296 +0,0 @@
-      SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
-*     .. Scalar Arguments ..
-      INTEGER INCX,N
-      CHARACTER DIAG,TRANS,UPLO
-*     ..
-*     .. Array Arguments ..
-      DOUBLE PRECISION AP(*),X(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DTPSV  solves one of the systems of equations
-*
-*     A*x = b,   or   A'*x = b,
-*
-*  where b and x are n element vectors and A is an n by n unit, or
-*  non-unit, upper or lower triangular matrix, supplied in packed form.
-*
-*  No test for singularity or near-singularity is included in this
-*  routine. Such tests must be performed before calling this routine.
-*
-*  Arguments
-*  ==========
-*
-*  UPLO   - CHARACTER*1.
-*           On entry, UPLO specifies whether the matrix is an upper or
-*           lower triangular matrix as follows:
-*
-*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
-*
-*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
-*
-*           Unchanged on exit.
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the equations to be solved as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   A*x = b.
-*
-*              TRANS = 'T' or 't'   A'*x = b.
-*
-*              TRANS = 'C' or 'c'   A'*x = b.
-*
-*           Unchanged on exit.
-*
-*  DIAG   - CHARACTER*1.
-*           On entry, DIAG specifies whether or not A is unit
-*           triangular as follows:
-*
-*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
-*
-*              DIAG = 'N' or 'n'   A is not assumed to be unit
-*                                  triangular.
-*
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the order of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  AP     - DOUBLE PRECISION array of DIMENSION at least
-*           ( ( n*( n + 1 ) )/2 ).
-*           Before entry with  UPLO = 'U' or 'u', the array AP must
-*           contain the upper triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
-*           respectively, and so on.
-*           Before entry with UPLO = 'L' or 'l', the array AP must
-*           contain the lower triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
-*           respectively, and so on.
-*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
-*           A are not referenced, but are assumed to be unity.
-*           Unchanged on exit.
-*
-*  X      - DOUBLE PRECISION array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ).
-*           Before entry, the incremented array X must contain the n
-*           element right-hand side vector b. On exit, X is overwritten
-*           with the solution vector x.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
-*     ..
-*     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
-      INTEGER I,INFO,IX,J,JX,K,KK,KX
-      LOGICAL NOUNIT
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
-          INFO = 1
-      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +         .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 2
-      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
-          INFO = 3
-      ELSE IF (N.LT.0) THEN
-          INFO = 4
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('DTPSV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF (N.EQ.0) RETURN
-*
-      NOUNIT = LSAME(DIAG,'N')
-*
-*     Set up the start point in X if the increment is not unity. This
-*     will be  ( N - 1 )*INCX  too small for descending loops.
-*
-      IF (INCX.LE.0) THEN
-          KX = 1 - (N-1)*INCX
-      ELSE IF (INCX.NE.1) THEN
-          KX = 1
-      END IF
-*
-*     Start the operations. In this version the elements of AP are
-*     accessed sequentially with one pass through AP.
-*
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  x := inv( A )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK - 1
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K - 1
-   10                     CONTINUE
-                      END IF
-                      KK = KK - J
-   20             CONTINUE
-              ELSE
-                  JX = KX + (N-1)*INCX
-                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 K = KK - 1,KK - J + 1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   30                     CONTINUE
-                      END IF
-                      JX = JX - INCX
-                      KK = KK - J
-   40             CONTINUE
-              END IF
-          ELSE
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK + 1
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K + 1
-   50                     CONTINUE
-                      END IF
-                      KK = KK + (N-J+1)
-   60             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 K = KK + 1,KK + N - J
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   70                     CONTINUE
-                      END IF
-                      JX = JX + INCX
-                      KK = KK + (N-J+1)
-   80             CONTINUE
-              END IF
-          END IF
-      ELSE
-*
-*        Form  x := inv( A' )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 100 J = 1,N
-                      TEMP = X(J)
-                      K = KK
-                      DO 90 I = 1,J - 1
-                          TEMP = TEMP - AP(K)*X(I)
-                          K = K + 1
-   90                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      X(J) = TEMP
-                      KK = KK + J
-  100             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 120 J = 1,N
-                      TEMP = X(JX)
-                      IX = KX
-                      DO 110 K = KK,KK + J - 2
-                          TEMP = TEMP - AP(K)*X(IX)
-                          IX = IX + INCX
-  110                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      X(JX) = TEMP
-                      JX = JX + INCX
-                      KK = KK + J
-  120             CONTINUE
-              END IF
-          ELSE
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 140 J = N,1,-1
-                      TEMP = X(J)
-                      K = KK
-                      DO 130 I = N,J + 1,-1
-                          TEMP = TEMP - AP(K)*X(I)
-                          K = K - 1
-  130                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      X(J) = TEMP
-                      KK = KK - (N-J+1)
-  140             CONTINUE
-              ELSE
-                  KX = KX + (N-1)*INCX
-                  JX = KX
-                  DO 160 J = N,1,-1
-                      TEMP = X(JX)
-                      IX = KX
-                      DO 150 K = KK,KK - (N- (J+1)),-1
-                          TEMP = TEMP - AP(K)*X(IX)
-                          IX = IX - INCX
-  150                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      X(JX) = TEMP
-                      JX = JX - INCX
-                      KK = KK - (N-J+1)
-  160             CONTINUE
-              END IF
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of DTPSV .
-*
-      END
diff --git a/blas/level2_impl.h b/blas/level2_impl.h
index 997ad016f..bd41f7e60 100644
--- a/blas/level2_impl.h
+++ b/blas/level2_impl.h
@@ -470,8 +470,55 @@ int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar
   *  No test for singularity or near-singularity is included in this
   *  routine. Such tests must be performed before calling this routine.
   */
-// int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx)
-// {
-//   return 1;
-// }
+int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx)
+{
+  typedef void (*functype)(int, const Scalar*, Scalar*);
+  static functype func[16];
+
+  static bool init = false;
+  if(!init)
+  {
+    for(int k=0; k<16; ++k)
+      func[k] = 0;
+
+    func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run);
+    func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run);
+    func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run);
+
+    func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run);
+    func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run);
+    func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run);
+
+    func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
+    func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
+    func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
+
+    func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
+    func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
+    func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::packed_triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
+
+    init = true;
+  }
+
+  Scalar* ap = reinterpret_cast<Scalar*>(pap);
+  Scalar* x = reinterpret_cast<Scalar*>(px);
+
+  int info = 0;
+  if(UPLO(*uplo)==INVALID)                                            info = 1;
+  else if(OP(*opa)==INVALID)                                          info = 2;
+  else if(DIAG(*diag)==INVALID)                                       info = 3;
+  else if(*n<0)                                                       info = 4;
+  else if(*incx==0)                                                   info = 7;
+  if(info)
+    return xerbla_(SCALAR_SUFFIX_UP"TPSV ",&info,6);
+
+  Scalar* actual_x = get_compact_vector(x,*n,*incx);
+
+  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
+  func[code](*n, ap, actual_x);
+
+  if(actual_x!=x) delete[] copy_back(actual_x,x,*n,*incx);
+
+  return 1;
+}
 
diff --git a/blas/stpsv.f b/blas/stpsv.f
deleted file mode 100644
index 7d95efbde..000000000
--- a/blas/stpsv.f
+++ /dev/null
@@ -1,296 +0,0 @@
-      SUBROUTINE STPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
-*     .. Scalar Arguments ..
-      INTEGER INCX,N
-      CHARACTER DIAG,TRANS,UPLO
-*     ..
-*     .. Array Arguments ..
-      REAL AP(*),X(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  STPSV  solves one of the systems of equations
-*
-*     A*x = b,   or   A'*x = b,
-*
-*  where b and x are n element vectors and A is an n by n unit, or
-*  non-unit, upper or lower triangular matrix, supplied in packed form.
-*
-*  No test for singularity or near-singularity is included in this
-*  routine. Such tests must be performed before calling this routine.
-*
-*  Arguments
-*  ==========
-*
-*  UPLO   - CHARACTER*1.
-*           On entry, UPLO specifies whether the matrix is an upper or
-*           lower triangular matrix as follows:
-*
-*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
-*
-*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
-*
-*           Unchanged on exit.
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the equations to be solved as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   A*x = b.
-*
-*              TRANS = 'T' or 't'   A'*x = b.
-*
-*              TRANS = 'C' or 'c'   A'*x = b.
-*
-*           Unchanged on exit.
-*
-*  DIAG   - CHARACTER*1.
-*           On entry, DIAG specifies whether or not A is unit
-*           triangular as follows:
-*
-*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
-*
-*              DIAG = 'N' or 'n'   A is not assumed to be unit
-*                                  triangular.
-*
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the order of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  AP     - REAL             array of DIMENSION at least
-*           ( ( n*( n + 1 ) )/2 ).
-*           Before entry with  UPLO = 'U' or 'u', the array AP must
-*           contain the upper triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
-*           respectively, and so on.
-*           Before entry with UPLO = 'L' or 'l', the array AP must
-*           contain the lower triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
-*           respectively, and so on.
-*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
-*           A are not referenced, but are assumed to be unity.
-*           Unchanged on exit.
-*
-*  X      - REAL             array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ).
-*           Before entry, the incremented array X must contain the n
-*           element right-hand side vector b. On exit, X is overwritten
-*           with the solution vector x.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      REAL ZERO
-      PARAMETER (ZERO=0.0E+0)
-*     ..
-*     .. Local Scalars ..
-      REAL TEMP
-      INTEGER I,INFO,IX,J,JX,K,KK,KX
-      LOGICAL NOUNIT
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
-          INFO = 1
-      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +         .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 2
-      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
-          INFO = 3
-      ELSE IF (N.LT.0) THEN
-          INFO = 4
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('STPSV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF (N.EQ.0) RETURN
-*
-      NOUNIT = LSAME(DIAG,'N')
-*
-*     Set up the start point in X if the increment is not unity. This
-*     will be  ( N - 1 )*INCX  too small for descending loops.
-*
-      IF (INCX.LE.0) THEN
-          KX = 1 - (N-1)*INCX
-      ELSE IF (INCX.NE.1) THEN
-          KX = 1
-      END IF
-*
-*     Start the operations. In this version the elements of AP are
-*     accessed sequentially with one pass through AP.
-*
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  x := inv( A )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK - 1
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K - 1
-   10                     CONTINUE
-                      END IF
-                      KK = KK - J
-   20             CONTINUE
-              ELSE
-                  JX = KX + (N-1)*INCX
-                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 K = KK - 1,KK - J + 1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   30                     CONTINUE
-                      END IF
-                      JX = JX - INCX
-                      KK = KK - J
-   40             CONTINUE
-              END IF
-          ELSE
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK + 1
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K + 1
-   50                     CONTINUE
-                      END IF
-                      KK = KK + (N-J+1)
-   60             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 K = KK + 1,KK + N - J
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   70                     CONTINUE
-                      END IF
-                      JX = JX + INCX
-                      KK = KK + (N-J+1)
-   80             CONTINUE
-              END IF
-          END IF
-      ELSE
-*
-*        Form  x := inv( A' )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 100 J = 1,N
-                      TEMP = X(J)
-                      K = KK
-                      DO 90 I = 1,J - 1
-                          TEMP = TEMP - AP(K)*X(I)
-                          K = K + 1
-   90                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      X(J) = TEMP
-                      KK = KK + J
-  100             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 120 J = 1,N
-                      TEMP = X(JX)
-                      IX = KX
-                      DO 110 K = KK,KK + J - 2
-                          TEMP = TEMP - AP(K)*X(IX)
-                          IX = IX + INCX
-  110                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      X(JX) = TEMP
-                      JX = JX + INCX
-                      KK = KK + J
-  120             CONTINUE
-              END IF
-          ELSE
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 140 J = N,1,-1
-                      TEMP = X(J)
-                      K = KK
-                      DO 130 I = N,J + 1,-1
-                          TEMP = TEMP - AP(K)*X(I)
-                          K = K - 1
-  130                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      X(J) = TEMP
-                      KK = KK - (N-J+1)
-  140             CONTINUE
-              ELSE
-                  KX = KX + (N-1)*INCX
-                  JX = KX
-                  DO 160 J = N,1,-1
-                      TEMP = X(JX)
-                      IX = KX
-                      DO 150 K = KK,KK - (N- (J+1)),-1
-                          TEMP = TEMP - AP(K)*X(IX)
-                          IX = IX - INCX
-  150                 CONTINUE
-                      IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      X(JX) = TEMP
-                      JX = JX - INCX
-                      KK = KK - (N-J+1)
-  160             CONTINUE
-              END IF
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of STPSV .
-*
-      END
diff --git a/blas/ztpsv.f b/blas/ztpsv.f
deleted file mode 100644
index b56e1d8c4..000000000
--- a/blas/ztpsv.f
+++ /dev/null
@@ -1,332 +0,0 @@
-      SUBROUTINE ZTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
-*     .. Scalar Arguments ..
-      INTEGER INCX,N
-      CHARACTER DIAG,TRANS,UPLO
-*     ..
-*     .. Array Arguments ..
-      DOUBLE COMPLEX AP(*),X(*)
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  ZTPSV  solves one of the systems of equations
-*
-*     A*x = b,   or   A'*x = b,   or   conjg( A' )*x = b,
-*
-*  where b and x are n element vectors and A is an n by n unit, or
-*  non-unit, upper or lower triangular matrix, supplied in packed form.
-*
-*  No test for singularity or near-singularity is included in this
-*  routine. Such tests must be performed before calling this routine.
-*
-*  Arguments
-*  ==========
-*
-*  UPLO   - CHARACTER*1.
-*           On entry, UPLO specifies whether the matrix is an upper or
-*           lower triangular matrix as follows:
-*
-*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
-*
-*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
-*
-*           Unchanged on exit.
-*
-*  TRANS  - CHARACTER*1.
-*           On entry, TRANS specifies the equations to be solved as
-*           follows:
-*
-*              TRANS = 'N' or 'n'   A*x = b.
-*
-*              TRANS = 'T' or 't'   A'*x = b.
-*
-*              TRANS = 'C' or 'c'   conjg( A' )*x = b.
-*
-*           Unchanged on exit.
-*
-*  DIAG   - CHARACTER*1.
-*           On entry, DIAG specifies whether or not A is unit
-*           triangular as follows:
-*
-*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
-*
-*              DIAG = 'N' or 'n'   A is not assumed to be unit
-*                                  triangular.
-*
-*           Unchanged on exit.
-*
-*  N      - INTEGER.
-*           On entry, N specifies the order of the matrix A.
-*           N must be at least zero.
-*           Unchanged on exit.
-*
-*  AP     - COMPLEX*16       array of DIMENSION at least
-*           ( ( n*( n + 1 ) )/2 ).
-*           Before entry with  UPLO = 'U' or 'u', the array AP must
-*           contain the upper triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
-*           respectively, and so on.
-*           Before entry with UPLO = 'L' or 'l', the array AP must
-*           contain the lower triangular matrix packed sequentially,
-*           column by column, so that AP( 1 ) contains a( 1, 1 ),
-*           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
-*           respectively, and so on.
-*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
-*           A are not referenced, but are assumed to be unity.
-*           Unchanged on exit.
-*
-*  X      - COMPLEX*16       array of dimension at least
-*           ( 1 + ( n - 1 )*abs( INCX ) ).
-*           Before entry, the incremented array X must contain the n
-*           element right-hand side vector b. On exit, X is overwritten
-*           with the solution vector x.
-*
-*  INCX   - INTEGER.
-*           On entry, INCX specifies the increment for the elements of
-*           X. INCX must not be zero.
-*           Unchanged on exit.
-*
-*  Further Details
-*  ===============
-*
-*  Level 2 Blas routine.
-*
-*  -- Written on 22-October-1986.
-*     Jack Dongarra, Argonne National Lab.
-*     Jeremy Du Croz, Nag Central Office.
-*     Sven Hammarling, Nag Central Office.
-*     Richard Hanson, Sandia National Labs.
-*
-*  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE COMPLEX ZERO
-      PARAMETER (ZERO= (0.0D+0,0.0D+0))
-*     ..
-*     .. Local Scalars ..
-      DOUBLE COMPLEX TEMP
-      INTEGER I,INFO,IX,J,JX,K,KK,KX
-      LOGICAL NOCONJ,NOUNIT
-*     ..
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC DCONJG
-*     ..
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
-          INFO = 1
-      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
-     +         .NOT.LSAME(TRANS,'C')) THEN
-          INFO = 2
-      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
-          INFO = 3
-      ELSE IF (N.LT.0) THEN
-          INFO = 4
-      ELSE IF (INCX.EQ.0) THEN
-          INFO = 7
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('ZTPSV ',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF (N.EQ.0) RETURN
-*
-      NOCONJ = LSAME(TRANS,'T')
-      NOUNIT = LSAME(DIAG,'N')
-*
-*     Set up the start point in X if the increment is not unity. This
-*     will be  ( N - 1 )*INCX  too small for descending loops.
-*
-      IF (INCX.LE.0) THEN
-          KX = 1 - (N-1)*INCX
-      ELSE IF (INCX.NE.1) THEN
-          KX = 1
-      END IF
-*
-*     Start the operations. In this version the elements of AP are
-*     accessed sequentially with one pass through AP.
-*
-      IF (LSAME(TRANS,'N')) THEN
-*
-*        Form  x := inv( A )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK - 1
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K - 1
-   10                     CONTINUE
-                      END IF
-                      KK = KK - J
-   20             CONTINUE
-              ELSE
-                  JX = KX + (N-1)*INCX
-                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 K = KK - 1,KK - J + 1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   30                     CONTINUE
-                      END IF
-                      JX = JX - INCX
-                      KK = KK - J
-   40             CONTINUE
-              END IF
-          ELSE
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK + 1
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K + 1
-   50                     CONTINUE
-                      END IF
-                      KK = KK + (N-J+1)
-   60             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 K = KK + 1,KK + N - J
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   70                     CONTINUE
-                      END IF
-                      JX = JX + INCX
-                      KK = KK + (N-J+1)
-   80             CONTINUE
-              END IF
-          END IF
-      ELSE
-*
-*        Form  x := inv( A' )*x  or  x := inv( conjg( A' ) )*x.
-*
-          IF (LSAME(UPLO,'U')) THEN
-              KK = 1
-              IF (INCX.EQ.1) THEN
-                  DO 110 J = 1,N
-                      TEMP = X(J)
-                      K = KK
-                      IF (NOCONJ) THEN
-                          DO 90 I = 1,J - 1
-                              TEMP = TEMP - AP(K)*X(I)
-                              K = K + 1
-   90                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      ELSE
-                          DO 100 I = 1,J - 1
-                              TEMP = TEMP - DCONJG(AP(K))*X(I)
-                              K = K + 1
-  100                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK+J-1))
-                      END IF
-                      X(J) = TEMP
-                      KK = KK + J
-  110             CONTINUE
-              ELSE
-                  JX = KX
-                  DO 140 J = 1,N
-                      TEMP = X(JX)
-                      IX = KX
-                      IF (NOCONJ) THEN
-                          DO 120 K = KK,KK + J - 2
-                              TEMP = TEMP - AP(K)*X(IX)
-                              IX = IX + INCX
-  120                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
-                      ELSE
-                          DO 130 K = KK,KK + J - 2
-                              TEMP = TEMP - DCONJG(AP(K))*X(IX)
-                              IX = IX + INCX
-  130                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK+J-1))
-                      END IF
-                      X(JX) = TEMP
-                      JX = JX + INCX
-                      KK = KK + J
-  140             CONTINUE
-              END IF
-          ELSE
-              KK = (N* (N+1))/2
-              IF (INCX.EQ.1) THEN
-                  DO 170 J = N,1,-1
-                      TEMP = X(J)
-                      K = KK
-                      IF (NOCONJ) THEN
-                          DO 150 I = N,J + 1,-1
-                              TEMP = TEMP - AP(K)*X(I)
-                              K = K - 1
-  150                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      ELSE
-                          DO 160 I = N,J + 1,-1
-                              TEMP = TEMP - DCONJG(AP(K))*X(I)
-                              K = K - 1
-  160                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK-N+J))
-                      END IF
-                      X(J) = TEMP
-                      KK = KK - (N-J+1)
-  170             CONTINUE
-              ELSE
-                  KX = KX + (N-1)*INCX
-                  JX = KX
-                  DO 200 J = N,1,-1
-                      TEMP = X(JX)
-                      IX = KX
-                      IF (NOCONJ) THEN
-                          DO 180 K = KK,KK - (N- (J+1)),-1
-                              TEMP = TEMP - AP(K)*X(IX)
-                              IX = IX - INCX
-  180                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
-                      ELSE
-                          DO 190 K = KK,KK - (N- (J+1)),-1
-                              TEMP = TEMP - DCONJG(AP(K))*X(IX)
-                              IX = IX - INCX
-  190                     CONTINUE
-                          IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK-N+J))
-                      END IF
-                      X(JX) = TEMP
-                      JX = JX - INCX
-                      KK = KK - (N-J+1)
-  200             CONTINUE
-              END IF
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of ZTPSV .
-*
-      END