mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-12 19:20:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement packed triangular solver.

Browse Source
This commit is contained in:
Chen-Pang He 2012-09-10 06:29:02 +08:00
parent 3642ca4d46
commit 65caa40a3d
8 changed files with 144 additions and 1264 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
blas/CMakeLists.txt
View File

@ -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()

88
blas/PackedTriangularSolverVector.h Normal file
View File

@ -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

1
blas/common.h
View File

@ -77,6 +77,7 @@ namespace Eigen {
#include "GeneralRank1Update.h"
#include "PackedSelfadjointProduct.h"
#include "PackedTriangularMatrixVector.h"
#include "PackedTriangularSolverVector.h"
#include "Rank2Update.h"
}

332
blas/ctpsv.f
View File

@ -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

296
blas/dtpsv.f
View File

@ -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

55
blas/level2_impl.h
View File

@ -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;
}

296
blas/stpsv.f
View File

@ -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

332
blas/ztpsv.f
View File

@ -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