From a2cf7ba9552bbe6b9371e1a2678914bab185968a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 13 Jul 2009 13:17:55 +0200 Subject: [PATCH] add triangular * vector product --- Eigen/Core | 1 + Eigen/src/Core/SolveTriangular.h | 4 +- Eigen/src/Core/TriangularMatrix.h | 15 +- .../Core/products/SelfadjointMatrixVector.h | 8 +- .../Core/products/TriangularMatrixVector.h | 161 ++++++++++++++++++ Eigen/src/Core/util/Macros.h | 8 +- test/CMakeLists.txt | 1 + test/product_triangular.cpp | 92 ++++++++++ test/triangular.cpp | 4 +- 9 files changed, 281 insertions(+), 13 deletions(-) create mode 100644 Eigen/src/Core/products/TriangularMatrixVector.h create mode 100644 test/product_triangular.cpp diff --git a/Eigen/Core b/Eigen/Core index 18e6a6045..8627ada53 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -181,6 +181,7 @@ namespace Eigen { #include "src/Core/SelfAdjointView.h" #include "src/Core/SolveTriangular.h" #include "src/Core/products/SelfadjointRank2Update.h" +#include "src/Core/products/TriangularMatrixVector.h" } // namespace Eigen diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 200b4a325..76e4289de 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -45,7 +45,7 @@ struct ei_triangular_solver_selector }; static void run(const Lhs& lhs, Rhs& other) { - static const int PanelWidth = EIGEN_TUNE_TRSV_PANEL_WIDTH; + static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH; const ActualLhsType& actualLhs = LhsProductTraits::extract(lhs); const int size = lhs.cols(); @@ -104,7 +104,7 @@ struct ei_triangular_solver_selector static void run(const Lhs& lhs, Rhs& other) { - static const int PanelWidth = EIGEN_TUNE_TRSV_PANEL_WIDTH; + static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH; const ActualLhsType& actualLhs = LhsProductTraits::extract(lhs); const int size = lhs.cols(); diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 81621dbcc..506be7965 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -142,13 +142,17 @@ struct ei_traits > : ei_traits }; }; -template class TriangularView - : public TriangularBase > +template +struct ei_triangular_vector_product_returntype; + +template class TriangularView + : public TriangularBase > { public: typedef TriangularBase Base; typedef typename ei_traits::Scalar Scalar; + typedef _MatrixType MatrixType; typedef typename MatrixType::PlainMatrixType PlainMatrixType; enum { @@ -243,6 +247,13 @@ template class TriangularView return res; } + template + ei_triangular_vector_product_returntype + operator*(const MatrixBase& rhs) const + { + return ei_triangular_vector_product_returntype(*this, rhs.derived(), 1); + } + template typename ei_plain_matrix_type_column_major::type solve(const MatrixBase& other) const; diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index aa3187a07..8ac83772c 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -25,7 +25,7 @@ #ifndef EIGEN_SELFADJOINT_MATRIX_VECTOR_H #define EIGEN_SELFADJOINT_MATRIX_VECTOR_H -/* Optimized col-major selfadjoint matrix * vector product: +/* Optimized selfadjoint matrix * vector product: * This algorithm processes 2 columns at onces that allows to both reduce * the number of load/stores of the result by a factor 2 and to reduce * the instruction dependency. @@ -78,11 +78,11 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector( size_t alignedStart = (starti) + ei_alignmentOffset(&res[starti], endi-starti); alignedEnd = alignedStart + ((endi-alignedStart)/(PacketSize))*(PacketSize); - res[j] += t0 * conj0(A0[j]); + res[j] += t0 * conj0(A0[j]); if(FirstTriangular) { - res[j+1] += t1 * conj0(A1[j+1]); - res[j] += t1 * conj0(A1[j]); + res[j+1] += t1 * conj0(A1[j+1]); + res[j] += t1 * conj0(A1[j]); t3 += conj1(A1[j]) * rhs[j]; } else diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h new file mode 100644 index 000000000..533aad170 --- /dev/null +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -0,0 +1,161 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 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_TRIANGULARMATRIXVECTOR_H +#define EIGEN_TRIANGULARMATRIXVECTOR_H + +template +struct ei_product_triangular_vector_selector; + +template +struct ei_product_triangular_vector_selector +{ + typedef typename Rhs::Scalar Scalar; + enum { + IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), + HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit + }; + static void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits::Scalar alpha) + { + static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH; + typename ei_conj_expr_if::ret cjLhs(lhs); + typename ei_conj_expr_if::ret cjRhs(rhs); + + int size = lhs.cols(); + for (int pi=0; pi0) + res.segment(s,r) += (alpha * cjRhs.coeff(i)) * cjLhs.col(i).segment(s,r); + if (HasUnitDiag) + res.coeffRef(i) += alpha * cjRhs.coeff(i); + } + int r = IsLowerTriangular ? size - pi - actualPanelWidth : pi; + if (r>0) + { + int s = IsLowerTriangular ? pi+actualPanelWidth : 0; + ei_cache_friendly_product_colmajor_times_vector( + r, + &(lhs.const_cast_derived().coeffRef(s,pi)), lhs.stride(), + rhs.segment(pi, actualPanelWidth), + &(res.coeffRef(s)), + alpha); + } + } + } +}; + +template +struct ei_product_triangular_vector_selector +{ + typedef typename Rhs::Scalar Scalar; + enum { + IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit), + HasUnitDiag = (Mode & UnitDiagBit)==UnitDiagBit + }; + static void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits::Scalar alpha) + { + static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH; + typename ei_conj_expr_if::ret cjLhs(lhs); + typename ei_conj_expr_if::ret cjRhs(rhs); + int size = lhs.cols(); + for (int pi=0; pi0) + res.coeffRef(i) += alpha * (cjLhs.row(i).segment(s,r).cwise() * cjRhs.segment(s,r).transpose()).sum(); + if (HasUnitDiag) + res.coeffRef(i) += alpha * cjRhs.coeff(i); + } + int r = IsLowerTriangular ? pi : size - pi - actualPanelWidth; + if (r>0) + { + int s = IsLowerTriangular ? 0 : pi + actualPanelWidth; + Block target(res,pi,0,actualPanelWidth,1); + ei_cache_friendly_product_rowmajor_times_vector( + &(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(), + &(rhs.const_cast_derived().coeffRef(s)), r, + target, alpha); + } + } + } +}; + +template +struct ei_triangular_vector_product_returntype + : public ReturnByValue, + Matrix::Scalar, + Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> > +{ + typedef typename Lhs::Scalar Scalar; + typedef typename ei_cleantype::type RhsNested; + ei_triangular_vector_product_returntype(const Lhs& lhs, const Rhs& rhs, Scalar alpha) + : m_lhs(lhs), m_rhs(rhs), m_alpha(alpha) + {} + + template void evalTo(Dest& dst) const + { + typedef typename Lhs::MatrixType MatrixType; + + typedef ei_blas_traits LhsBlasTraits; + typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType; + typedef typename ei_cleantype::type _ActualLhsType; + const ActualLhsType actualLhs = LhsBlasTraits::extract(m_lhs._expression()); + + typedef ei_blas_traits RhsBlasTraits; + typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType; + typedef typename ei_cleantype::type _ActualRhsType; + const ActualRhsType actualRhs = RhsBlasTraits::extract(m_rhs); + + Scalar actualAlpha = m_alpha * LhsBlasTraits::extractScalarFactor(m_lhs._expression()) + * RhsBlasTraits::extractScalarFactor(m_rhs); + + dst.resize(m_rhs.rows(), m_rhs.cols()); + dst.setZero(); + ei_product_triangular_vector_selector + <_ActualLhsType,_ActualRhsType,Dest, + ei_traits::Mode, + LhsBlasTraits::NeedToConjugate, + RhsBlasTraits::NeedToConjugate, + ei_traits::Flags&RowMajorBit> + ::run(actualLhs,actualRhs,dst,actualAlpha); + } + + const Lhs m_lhs; + const typename Rhs::Nested m_rhs; + const Scalar m_alpha; +}; + +#endif // EIGEN_TRIANGULARMATRIXVECTOR_H diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index bb73b03cc..0cc5ae9aa 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -94,11 +94,11 @@ #define EIGEN_TUNE_FOR_CPU_CACHE_SIZE (sizeof(float)*256*256) #endif -/** Defines the maximal width of the blocks used in the triangular solver - * for vectors (level 2 blas xTRSV). The default is 8. +/** Defines the maximal width of the blocks used in the triangular product and solver + * for vectors (level 2 blas xTRMV and xTRSV). The default is 8. */ -#ifndef EIGEN_TUNE_TRSV_PANEL_WIDTH -#define EIGEN_TUNE_TRSV_PANEL_WIDTH 8 +#ifndef EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH +#define EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH 8 #endif /** Allows to disable some optimizations which might affect the accuracy of the result. diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 229c7123b..886731d45 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -110,6 +110,7 @@ ei_add_test(array) ei_add_test(array_replicate) ei_add_test(array_reverse) ei_add_test(triangular) +ei_add_test(product_triangular) ei_add_test(cholesky " " "${GSL_LIBRARIES}") ei_add_test(lu ${EI_OFLAG}) ei_add_test(determinant) diff --git a/test/product_triangular.cpp b/test/product_triangular.cpp new file mode 100644 index 000000000..876fb4388 --- /dev/null +++ b/test/product_triangular.cpp @@ -0,0 +1,92 @@ +// This file is triangularView of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 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 . + +#include "main.h" + +template void product_triangular(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix VectorType; + + RealScalar largerEps = 10*test_precision(); + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m3(rows, cols); + VectorType v1 = VectorType::Random(rows); + + Scalar s1 = ei_random(); + + m1 = MatrixType::Random(rows, cols); + + // check with a column-major matrix + m3 = m1.template triangularView(); + VERIFY((m3 * v1).isApprox(m1.template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3 * v1).isApprox(m1.template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3 * v1).isApprox(m1.template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3 * v1).isApprox(m1.template triangularView() * v1, largerEps)); + + // check conjugated and scalar multiple expressions (col-major) + m3 = m1.template triangularView(); + VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView() * v1.conjugate(), largerEps)); + + // check with a row-major matrix + m3 = m1.template triangularView(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView() * v1, largerEps)); + + // check conjugated and scalar multiple expressions (row-major) + m3 = m1.template triangularView(); + VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView() * v1, largerEps)); + m3 = m1.template triangularView(); + VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView() * (s1*v1.conjugate()), largerEps)); + m3 = m1.template triangularView(); + + // TODO check with sub-matrices +} + +void test_product_triangular() +{ + for(int i = 0; i < g_repeat ; i++) { + CALL_SUBTEST( product_triangular(Matrix()) ); + CALL_SUBTEST( product_triangular(Matrix()) ); + CALL_SUBTEST( product_triangular(Matrix3d()) ); + CALL_SUBTEST( product_triangular(Matrix,23, 23>()) ); + CALL_SUBTEST( product_triangular(MatrixXcd(17,17)) ); + CALL_SUBTEST( product_triangular(Matrix(19, 19)) ); + } +} diff --git a/test/triangular.cpp b/test/triangular.cpp index 7c680a8ed..6385bffd1 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -51,6 +51,8 @@ template void triangular(const MatrixType& m) v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); + Scalar s1 = ei_random(); + MatrixType m1up = m1.template triangularView(); MatrixType m2up = m2.template triangularView(); @@ -142,7 +144,7 @@ void test_triangular() CALL_SUBTEST( triangular(Matrix3d()) ); CALL_SUBTEST( triangular(MatrixXcf(4, 4)) ); CALL_SUBTEST( triangular(Matrix,8, 8>()) ); - CALL_SUBTEST( triangular(MatrixXd(17,17)) ); + CALL_SUBTEST( triangular(MatrixXcd(17,17)) ); CALL_SUBTEST( triangular(Matrix(5, 5)) ); } }