Implement packed triangular solver.
diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt
index 171b75a..3877e12 100644
--- a/blas/CMakeLists.txt
+++ b/blas/CMakeLists.txt
@@ -18,10 +18,10 @@
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 0000000..5c0bb4b
--- /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 1019d86..ee38b31 100644
--- a/blas/common.h
+++ b/blas/common.h
@@ -77,6 +77,7 @@
#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 1804797..0000000
--- 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 c7e58d3..0000000
--- 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 997ad01..bd41f7e 100644
--- a/blas/level2_impl.h
+++ b/blas/level2_impl.h
@@ -470,8 +470,55 @@
* 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 7d95efb..0000000
--- 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 b56e1d8..0000000
--- 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
