From 3a4c78b588ff523cb07bd7068cbe857b9b6a7ded Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 1 Dec 2011 18:06:28 +0100 Subject: add code for band triangular problems: - currently available from the BLAS interface only - and for vectors only --- blas/BandTriangularSolver.h | 112 ++++++++++++++++++++++++++++++++++++++++ blas/CMakeLists.txt | 9 ++-- blas/common.h | 6 +++ blas/level2_impl.h | 121 +++++++++++++++++++++++++++++++++++++++++--- 4 files changed, 236 insertions(+), 12 deletions(-) create mode 100644 blas/BandTriangularSolver.h (limited to 'blas') diff --git a/blas/BandTriangularSolver.h b/blas/BandTriangularSolver.h new file mode 100644 index 000000000..420221aea --- /dev/null +++ b/blas/BandTriangularSolver.h @@ -0,0 +1,112 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_BAND_TRIANGULARSOLVER_H +#define EIGEN_BAND_TRIANGULARSOLVER_H + +namespace internal { + + /* \internal + * Solve Ax=b with A a band triangular matrix + * TODO: extend it to matrices for x abd b */ +template +struct band_solve_triangular_selector; + + +template +struct band_solve_triangular_selector +{ + typedef Map, 0, OuterStride<> > LhsMap; + typedef Map > RhsMap; + enum { IsLower = (Mode&Lower) ? 1 : 0 }; + static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) + { + const LhsMap lhs(_lhs,size,k+1,OuterStride<>(lhsStride)); + RhsMap other(_other,size,1); + typename internal::conditional< + ConjLhs, + const CwiseUnaryOp,LhsMap>, + const LhsMap&> + ::type cjLhs(lhs); + + for(int col=0 ; col0) + other.coeffRef(i,col) -= cjLhs.row(i).segment(actual_start,actual_k).transpose() + .cwiseProduct(other.col(col).segment(IsLower ? i-actual_k : i+1,actual_k)).sum(); + + if((Mode&UnitDiag)==0) + other.coeffRef(i,col) /= cjLhs(i,IsLower ? k : 0); + } + } + } + +}; + +template +struct band_solve_triangular_selector +{ + typedef Map, 0, OuterStride<> > LhsMap; + typedef Map > RhsMap; + enum { IsLower = (Mode&Lower) ? 1 : 0 }; + static void run(Index size, Index k, const LhsScalar* _lhs, Index lhsStride, RhsScalar* _other) + { + const LhsMap lhs(_lhs,k+1,size,OuterStride<>(lhsStride)); + RhsMap other(_other,size,1); + typename internal::conditional< + ConjLhs, + const CwiseUnaryOp,LhsMap>, + const LhsMap&> + ::type cjLhs(lhs); + + for(int col=0 ; col0) + other.col(col).segment(IsLower ? i+1 : i-actual_k, actual_k) + -= other.coeff(i,col) * cjLhs.col(i).segment(actual_start,actual_k); + + } + } + } +}; + + +} // end namespace internal + +#endif // EIGEN_BAND_TRIANGULARSOLVER_H diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt index 55440f4dc..453d5874c 100644 --- a/blas/CMakeLists.txt +++ b/blas/CMakeLists.txt @@ -18,10 +18,11 @@ if(EIGEN_Fortran_COMPILER_WORKS) set(EigenBlas_SRCS ${EigenBlas_SRCS} complexdots.f srotm.f srotmg.f drotm.f drotmg.f - lsame.f chpr2.f ctbsv.f dspmv.f dtbmv.f dtpsv.f ssbmv.f sspr.f stpmv.f - zhpr2.f ztbsv.f chbmv.f chpr.f ctpmv.f dspr2.f dtbsv.f sspmv.f stbmv.f stpsv.f - zhbmv.f zhpr.f ztpmv.f chpmv.f ctbmv.f ctpsv.f dsbmv.f dspr.f dtpmv.f sspr2.f - stbsv.f zhpmv.f ztbmv.f ztpsv.f + lsame.f chpr2.f dspmv.f dtpsv.f ssbmv.f sspr.f stpmv.f + zhpr2.f chbmv.f chpr.f ctpmv.f dspr2.f sspmv.f stpsv.f + zhbmv.f zhpr.f ztpmv.f chpmv.f ctpsv.f dsbmv.f dspr.f dtpmv.f sspr2.f + zhpmv.f ztpsv.f + dtbmv.f stbmv.f ctbmv.f ztbmv.f ) else() diff --git a/blas/common.h b/blas/common.h index a63a5d286..714bce574 100644 --- a/blas/common.h +++ b/blas/common.h @@ -93,6 +93,12 @@ inline bool check_uplo(const char* uplo) #include #include + + +namespace Eigen { +#include "BandTriangularSolver.h" +} + using namespace Eigen; typedef SCALAR Scalar; diff --git a/blas/level2_impl.h b/blas/level2_impl.h index 8cbc2f424..0781fa56a 100644 --- a/blas/level2_impl.h +++ b/blas/level2_impl.h @@ -271,6 +271,7 @@ int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealSca return 0; } +#if 0 /** TBMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, @@ -278,10 +279,56 @@ int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealSca * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular band matrix, with ( k + 1 ) diagonals. */ -// int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx) -// { -// return 1; -// } +int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) +{ + Scalar* a = reinterpret_cast(pa); + Scalar* x = reinterpret_cast(px); + int coeff_rows = *k + 1; + + 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(*k<0) info = 5; + else if(*lda::run); + func[TR | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[ADJ | (UP << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector::run); + + func[NOTR | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[TR | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[ADJ | (LO << 2) | (NUNIT << 3)] = (internal::band_solve_triangular_selector::run); + + func[NOTR | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[TR | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[ADJ | (UP << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector::run); + + func[NOTR | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[TR | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector::run); + func[ADJ | (LO << 2) | (UNIT << 3)] = (internal::band_solve_triangular_selector::run); + + init = true; + } + + Scalar* a = reinterpret_cast(pa); + Scalar* x = reinterpret_cast(px); + int coeff_rows = *k+1; + + int info = 0; + if(UPLO(*uplo)==INVALID) info = 1; + else if(OP(*op)==INVALID) info = 2; + else if(DIAG(*diag)==INVALID) info = 3; + else if(*n<0) info = 4; + else if(*k<0) info = 5; + else if(*lda=16 || func[code]==0) + return 0; + + func[code](*n, *k, a, *lda, actual_x); + + if(actual_x!=x) delete[] copy_back(actual_x,x,actual_n,*incx); + + return 0; +} /** DTPMV performs one of the matrix-vector operations * -- cgit v1.2.3