--- blas2.hpp.orig	2009-03-03 12:47:04.000000000 +0000
+++ blas2.hpp	2009-03-03 12:48:22.000000000 +0000
@@ -5,7 +5,7 @@
 // (See accompanying file LICENSE_1_0.txt or copy at
 // http://www.boost.org/LICENSE_1_0.txt)
 //
-// Modified to include dtrsv by Kian Ming A. Chai (3 Mar 2009)
+// Modified to include xtrsv and xtpsv by Kian Ming A. Chai (3 Mar 2009)
 
 #ifndef BOOST_BINDINGS_BLAS_BLAS2_HPP
 #define BOOST_BINDINGS_BLAS_BLAS2_HPP
@@ -106,6 +106,46 @@
     detail::trsv( UPLO, TRANS, DIAG, n, a_ptr, lda, x_ptr, stride_x);
   }
 
+  /*
+   *  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.
+   *
+   */
+  // ! CAUTION this function assumes that all matrices involved are column-major matrices
+  template < typename matrix_type, typename vector_type >
+  inline
+  void tpsv(const char TRANS, const matrix_type &ap, const vector_type &x)
+  {
+    typedef typename traits::matrix_traits<matrix_type> mtraits;
+    typedef typename mtraits::value_type value_type;
+    BOOST_STATIC_ASSERT( ( boost::is_same< typename traits::vector_traits<vector_type>::value_type,
+			   value_type>::value ) );
+    // precondition: matrix_type must be dense or dense_proxy
+    /* not all compilers can handle the traits
+    BOOST_STATIC_ASSERT( ( boost::is_same< typename mtraits::matrix_structure,
+			   boost::numeric::bindings::traits::general_t
+                           >::value ) ) ;
+    */
+    const int m = traits::matrix_size1( ap ) ;
+    const int n = traits::matrix_size2( ap ) ;
+    assert( n == m ); // must be a square matrix also
+    assert ( traits::vector_size( x ) == n  ) ; 
+    char uplo = traits::matrix_uplo_tag (ap);
+    const int stride_x = traits::vector_stride( x ) ;
+
+    const value_type *ap_ptr = traits::matrix_storage( ap ) ;
+    const value_type *x_ptr = traits::vector_storage( x ) ;
+
+    detail::tpsv( uplo, TRANS, 'N', n, ap_ptr, x_ptr, stride_x );
+  }
+
 /*
   // A <- alpha * x * trans(y) ( outer product ), alpha, x and y are complex-valued 
   template < typename vector_type_x, typename vector_type_y, typename value_type, typename matrix_type >