From 65caa40a3d792dbbd28b3f47de2a87efea58bb24 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Mon, 10 Sep 2012 06:29:02 +0800 Subject: Implement packed triangular solver. --- blas/CMakeLists.txt | 8 +- blas/PackedTriangularSolverVector.h | 88 ++++++++++ blas/common.h | 1 + blas/ctpsv.f | 332 ------------------------------------ blas/dtpsv.f | 296 -------------------------------- blas/level2_impl.h | 55 +++++- blas/stpsv.f | 296 -------------------------------- blas/ztpsv.f | 332 ------------------------------------ 8 files changed, 144 insertions(+), 1264 deletions(-) create mode 100644 blas/PackedTriangularSolverVector.h delete mode 100644 blas/ctpsv.f delete mode 100644 blas/dtpsv.f delete mode 100644 blas/stpsv.f delete mode 100644 blas/ztpsv.f (limited to 'blas') 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 +// +// 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 +struct packed_triangular_solve_vector; + +// forward and backward substitution, row-major, rhs is a vector +template +struct packed_triangular_solve_vector +{ + enum { + IsLower = (Mode&Lower)==Lower + }; + static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) + { + internal::conj_if cj; + typedef Map > LhsMap; + typedef typename conj_expr_if::type ConjLhsType; + + lhs += IsLower ? 0 : (size*(size+1)>>1)-1; + for(Index pi=0; pi0) + rhs[i] -= (ConjLhsType(LhsMap(lhs+s,pi)) + .cwiseProduct(Map >(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 +struct packed_triangular_solve_vector +{ + enum { + IsLower = (Mode&Lower)==Lower + }; + static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) + { + internal::conj_if cj; + typedef Map > LhsMap; + typedef typename conj_expr_if::type ConjLhsType; + + lhs += IsLower ? 0 : size*(size-1)>>1; + for(Index pi=0; pi0) + Map >(rhs+(IsLower? i+1 : 0),r) -= + rhs[i] * ConjLhsType(LhsMap(lhs+(IsLower? 1 : 0),r)); + IsLower ? lhs += size-pi : lhs -= r; + } + } +}; + +template +struct packed_triangular_solve_vector +{ + static void run(Index size, const LhsScalar* lhs, RhsScalar* rhs) + { + packed_triangular_solve_vector::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::run); + func[TR | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector::run); + + func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[TR | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::packed_triangular_solve_vector::run); + + func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[TR | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector::run); + + func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[TR | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector::run); + func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::packed_triangular_solve_vector::run); + + init = true; + } + + Scalar* ap = reinterpret_cast(pap); + Scalar* x = reinterpret_cast(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 -- cgit v1.2.3