From f95b77be6216db7b5448d7d728339cae81129bc9 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 27 Jul 2009 11:42:54 +0200 Subject: [PATCH] trmm is now fully working and available via TriangularView::operator* --- Eigen/src/Core/TriangularMatrix.h | 23 +- .../Core/products/TriangularMatrixMatrix.h | 392 ++++++++++++++++++ .../Core/products/TriangularMatrixVector.h | 112 +++-- test/CMakeLists.txt | 11 +- test/product_trmm.cpp | 69 +++ ...roduct_triangular.cpp => product_trmv.cpp} | 0 test/product_trsm.cpp | 1 + 7 files changed, 575 insertions(+), 33 deletions(-) create mode 100644 Eigen/src/Core/products/TriangularMatrixMatrix.h create mode 100644 test/product_trmm.cpp rename test/{product_triangular.cpp => product_trmv.cpp} (100%) diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index d2f1e6c28..861b738cb 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -142,8 +142,10 @@ struct ei_traits > : ei_traits }; }; -template -struct ei_triangular_vector_product_returntype; +template +struct ei_triangular_product_returntype; template class TriangularView : public TriangularBase > @@ -247,11 +249,24 @@ template class TriangularView return res; } + /** Efficient triangular matrix times vector/matrix product */ template - ei_triangular_vector_product_returntype + ei_triangular_product_returntype operator*(const MatrixBase& rhs) const { - return ei_triangular_vector_product_returntype(*this, rhs.derived(), 1); + return ei_triangular_product_returntype + + (m_matrix, rhs.derived()); + } + + /** Efficient vector/matrix times triangular matrix product */ + template friend + ei_triangular_product_returntype + operator*(const MatrixBase& lhs, const TriangularView& rhs) + { + return ei_triangular_product_returntype + + (lhs.derived(),rhs.m_matrix); } template diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h new file mode 100644 index 000000000..43a4c3d18 --- /dev/null +++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h @@ -0,0 +1,392 @@ +// 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_TRIANGULAR_MATRIX_MATRIX_H +#define EIGEN_TRIANGULAR_MATRIX_MATRIX_H + +// template +// struct ei_gemm_pack_lhs_triangular +// { +// Matrix::IsComplex && Conjugate> cj; +// ei_const_blas_data_mapper lhs(_lhs,lhsStride); +// int count = 0; +// const int peeled_mc = (rows/mr)*mr; +// for(int i=0; i +struct ei_product_triangular_matrix_matrix; + +template +struct ei_product_triangular_matrix_matrix +{ + static EIGEN_STRONG_INLINE void run( + int size, int otherSize, + const Scalar* lhs, int lhsStride, + const Scalar* rhs, int rhsStride, + Scalar* res, int resStride, + Scalar alpha) + { + ei_product_triangular_matrix_matrix + ::run(size, otherSize, rhs, rhsStride, lhs, lhsStride, res, resStride, alpha); + } +}; + +// implements col-major += alpha * op(triangular) * op(general) +template +struct ei_product_triangular_matrix_matrix +{ + + static EIGEN_DONT_INLINE void run( + int size, int cols, + const Scalar* _lhs, int lhsStride, + const Scalar* _rhs, int rhsStride, + Scalar* res, int resStride, + Scalar alpha) + { + int rows = size; + + ei_const_blas_data_mapper lhs(_lhs,lhsStride); + ei_const_blas_data_mapper rhs(_rhs,rhsStride); + + if (ConjugateRhs) + alpha = ei_conj(alpha); + + typedef ei_product_blocking_traits Blocking; + enum { + SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), + IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + }; + + int kc = std::min(Blocking::Max_kc/4,size); // cache block size along the K direction + int mc = std::min(Blocking::Max_mc,rows); // cache block size along the M direction + + Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc); + Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*Blocking::PacketSize); + + Matrix triangularBuffer; + triangularBuffer.setZero(); + triangularBuffer.diagonal().setOnes(); + + ei_gebp_kernel > gebp_kernel; + ei_gemm_pack_lhs pack_lhs; + ei_gemm_pack_rhs pack_rhs; + + for(int k2=IsLowerTriangular ? size : 0; + IsLowerTriangular ? k2>0 : k2 skip it + // 2 - the diagonal block => special kernel + // 3 - the panel below the diagonal block => GEPP + // the block diagonal + { + // for each small vertical panels of lhs + for (int k1=0; k1(actual_kc-k1, SmallPanelWidth); + int lengthTarget = IsLowerTriangular ? actual_kc-k1-actualPanelWidth : k1; + int startBlock = actual_k2+k1; + int blockBOffset = k1; + + // => GEBP with the micro triangular block + // The trick is to pack this micro block while filling the opposite triangular part with zeros. + // To this end we do an extra triangular copy to small temporary buffer + for (int k=0;k0) + { + int startTarget = IsLowerTriangular ? actual_k2+k1+actualPanelWidth : actual_k2; + + pack_lhs(blockA, &lhs(startTarget,startBlock), lhsStride, actualPanelWidth, lengthTarget); + + gebp_kernel(res+startTarget, resStride, blockA, blockB, lengthTarget, actualPanelWidth, cols, + actualPanelWidth, actual_kc, 0, blockBOffset*Blocking::PacketSize); + } + } + } + // the part below the diagonal => GEPP + { + int start = IsLowerTriangular ? k2 : 0; + int end = IsLowerTriangular ? size : actual_k2; + for(int i2=start; i2() + (blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc); + + gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols); + } + } + } + + ei_aligned_stack_delete(Scalar, blockA, kc*mc); + ei_aligned_stack_delete(Scalar, blockB, kc*cols*Blocking::PacketSize); + } +}; + +// implements col-major += alpha * op(general) * op(triangular) +template +struct ei_product_triangular_matrix_matrix +{ + + static EIGEN_DONT_INLINE void run( + int size, int rows, + const Scalar* _lhs, int lhsStride, + const Scalar* _rhs, int rhsStride, + Scalar* res, int resStride, + Scalar alpha) + { + int cols = size; + + ei_const_blas_data_mapper lhs(_lhs,lhsStride); + ei_const_blas_data_mapper rhs(_rhs,rhsStride); + + if (ConjugateRhs) + alpha = ei_conj(alpha); + + typedef ei_product_blocking_traits Blocking; + enum { + SmallPanelWidth = EIGEN_ENUM_MAX(Blocking::mr,Blocking::nr), + IsLowerTriangular = (Mode&LowerTriangular) == LowerTriangular + }; + + int kc = std::min(Blocking::Max_kc/4,size); // cache block size along the K direction + int mc = std::min(Blocking::Max_mc,rows); // cache block size along the M direction + + Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc); + Scalar* blockB = ei_aligned_stack_new(Scalar, kc*cols*Blocking::PacketSize); + + Matrix triangularBuffer; + triangularBuffer.setZero(); + triangularBuffer.diagonal().setOnes(); + + ei_gebp_kernel > gebp_kernel; + ei_gemm_pack_lhs pack_lhs; + ei_gemm_pack_rhs pack_rhs; + ei_gemm_pack_rhs pack_rhs_panel; + + for(int k2=IsLowerTriangular ? 0 : size; + IsLowerTriangular ? k20; + IsLowerTriangular ? k2+=kc : k2-=kc) + { + const int actual_kc = std::min(IsLowerTriangular ? size-k2 : k2, kc); + int actual_k2 = IsLowerTriangular ? k2 : k2-actual_kc; + int rs = IsLowerTriangular ? actual_k2 : size - k2; + Scalar* geb = blockB+actual_kc*actual_kc*Blocking::PacketSize; + + pack_rhs(geb, &rhs(actual_k2,IsLowerTriangular ? 0 : k2), rhsStride, alpha, actual_kc, rs); + + // pack the triangular part of the rhs padding the unrolled blocks with zeros + { + for (int j2=0; j2(actual_kc-j2, SmallPanelWidth); + int actual_j2 = actual_k2 + j2; + int panelOffset = IsLowerTriangular ? j2+actualPanelWidth : 0; + int panelLength = IsLowerTriangular ? actual_kc-j2-actualPanelWidth : j2; + // general part + pack_rhs_panel(blockB+j2*actual_kc*Blocking::PacketSize, + &rhs(actual_k2+panelOffset, actual_j2), rhsStride, alpha, + panelLength, actualPanelWidth, + actual_kc, panelOffset); + + // append the triangular part via a temporary buffer + for (int j=0;j(actual_kc-j2, SmallPanelWidth); + int panelLength = IsLowerTriangular ? actual_kc-j2 : j2+actualPanelWidth; + int blockOffset = IsLowerTriangular ? j2 : 0; + + gebp_kernel(res+i2+(actual_k2+j2)*resStride, resStride, + blockA, blockB+j2*actual_kc*Blocking::PacketSize, + actual_mc, panelLength, actualPanelWidth, + actual_kc, actual_kc, // strides + blockOffset, blockOffset*Blocking::PacketSize);// offsets + } + } + gebp_kernel(res+i2+(IsLowerTriangular ? 0 : k2)*resStride, resStride, + blockA, geb, actual_mc, actual_kc, rs); + } + } + + ei_aligned_stack_delete(Scalar, blockA, kc*mc); + ei_aligned_stack_delete(Scalar, blockB, kc*cols*Blocking::PacketSize); + } +}; + +/*************************************************************************** +* Wrapper to ei_product_triangular_matrix_matrix +***************************************************************************/ + +template +struct ei_triangular_product_returntype + : public ReturnByValue, + Matrix::Scalar, + Lhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> > +{ + ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + {} + + typedef typename Lhs::Scalar Scalar; + + typedef typename Lhs::Nested LhsNested; + typedef typename ei_cleantype::type _LhsNested; + typedef ei_blas_traits<_LhsNested> LhsBlasTraits; + typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType; + typedef typename ei_cleantype::type _ActualLhsType; + + typedef typename Rhs::Nested RhsNested; + typedef typename ei_cleantype::type _RhsNested; + typedef ei_blas_traits<_RhsNested> RhsBlasTraits; + typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType; + typedef typename ei_cleantype::type _ActualRhsType; + +// enum { +// LhsUpLo = LhsMode&(UpperTriangularBit|LowerTriangularBit), +// LhsIsTriangular = (LhsMode&SelfAdjointBit)==SelfAdjointBit, +// RhsUpLo = RhsMode&(UpperTriangularBit|LowerTriangularBit), +// RhsIsSelfAdjoint = (RhsMode&SelfAdjointBit)==SelfAdjointBit +// }; + + template inline void _addTo(Dest& dst) const + { evalTo(dst,1); } + template inline void _subTo(Dest& dst) const + { evalTo(dst,-1); } + + template void evalTo(Dest& dst) const + { + dst.resize(m_lhs.rows(), m_rhs.cols()); + dst.setZero(); + evalTo(dst,1); + } + + template void evalTo(Dest& dst, Scalar alpha) const + { + const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs); + const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs); + + Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs) + * RhsBlasTraits::extractScalarFactor(m_rhs); + + ei_product_triangular_matrix_matrix::Flags&RowMajorBit) ? RowMajor : ColMajor, LhsBlasTraits::NeedToConjugate, + (ei_traits<_ActualRhsType>::Flags&RowMajorBit) ? RowMajor : ColMajor, RhsBlasTraits::NeedToConjugate, + (ei_traits::Flags&RowMajorBit) ? RowMajor : ColMajor> + ::run( + lhs.rows(), LhsIsTriangular ? rhs.cols() : lhs.rows(), // sizes + &lhs.coeff(0,0), lhs.stride(), // lhs info + &rhs.coeff(0,0), rhs.stride(), // rhs info + &dst.coeffRef(0,0), dst.stride(), // result info + actualAlpha // alpha + ); + } + + const LhsNested m_lhs; + const RhsNested m_rhs; +}; + +#endif // EIGEN_TRIANGULAR_MATRIX_MATRIX_H diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index 533aad170..0fbbb50d2 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -113,49 +113,113 @@ struct ei_product_triangular_vector_selector -struct ei_triangular_vector_product_returntype - : public ReturnByValue, +// 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; +// }; + + +/*************************************************************************** +* Wrapper to ei_product_triangular_vector +***************************************************************************/ + +template +struct ei_triangular_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) + + typedef typename Lhs::Nested LhsNested; + typedef typename ei_cleantype::type _LhsNested; + typedef ei_blas_traits<_LhsNested> LhsBlasTraits; + typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType; + typedef typename ei_cleantype::type _ActualLhsType; + + typedef typename Rhs::Nested RhsNested; + typedef typename ei_cleantype::type _RhsNested; + typedef ei_blas_traits<_RhsNested> RhsBlasTraits; + typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType; + typedef typename ei_cleantype::type _ActualRhsType; + + ei_triangular_product_returntype(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) {} + template inline void _addTo(Dest& dst) const + { evalTo(dst,1); } + template inline void _subTo(Dest& dst) const + { evalTo(dst,-1); } + 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()); + dst.resize(m_lhs.rows(), m_rhs.cols()); + dst.setZero(); + evalTo(dst,1); + } - typedef ei_blas_traits RhsBlasTraits; - typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType; - typedef typename ei_cleantype::type _ActualRhsType; - const ActualRhsType actualRhs = RhsBlasTraits::extract(m_rhs); + template void evalTo(Dest& dst, Scalar alpha) const + { + const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs); + const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs); + + Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs) + * RhsBlasTraits::extractScalarFactor(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, + Mode, LhsBlasTraits::NeedToConjugate, RhsBlasTraits::NeedToConjugate, ei_traits::Flags&RowMajorBit> - ::run(actualLhs,actualRhs,dst,actualAlpha); + ::run(lhs,rhs,dst,actualAlpha); } - const Lhs m_lhs; - const typename Rhs::Nested m_rhs; - const Scalar m_alpha; + const LhsNested m_lhs; + const RhsNested m_rhs; }; #endif // EIGEN_TRIANGULARMATRIXVECTOR_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 462032453..99224ff60 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -98,10 +98,6 @@ ei_add_test(redux) ei_add_test(product_small) ei_add_test(product_large ${EI_OFLAG}) ei_add_test(product_extra ${EI_OFLAG}) -ei_add_test(product_selfadjoint ${EI_OFLAG}) -ei_add_test(product_symm ${EI_OFLAG}) -ei_add_test(product_syrk ${EI_OFLAG}) -ei_add_test(product_trsm ${EI_OFLAG}) ei_add_test(diagonalmatrices) ei_add_test(adjoint) ei_add_test(submatrices) @@ -113,7 +109,12 @@ 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(product_selfadjoint ${EI_OFLAG}) +ei_add_test(product_symm ${EI_OFLAG}) +ei_add_test(product_syrk ${EI_OFLAG}) +ei_add_test(product_trmv ${EI_OFLAG}) +ei_add_test(product_trmm ${EI_OFLAG}) +ei_add_test(product_trsm ${EI_OFLAG}) ei_add_test(bandmatrix) ei_add_test(cholesky " " "${GSL_LIBRARIES}") ei_add_test(lu ${EI_OFLAG}) diff --git a/test/product_trmm.cpp b/test/product_trmm.cpp new file mode 100644 index 000000000..47ffb4af3 --- /dev/null +++ b/test/product_trmm.cpp @@ -0,0 +1,69 @@ +// This file is part 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 trmm(int size,int othersize) +{ + typedef typename NumTraits::Real RealScalar; + + Matrix tri(size,size), upTri(size,size), loTri(size,size); + Matrix ge1(size,othersize), ge2(10,size), ge3; + Matrix rge3; + + Scalar s1 = ei_random(), + s2 = ei_random(); + + tri.setRandom(); + loTri = tri.template triangularView(); + upTri = tri.template triangularView(); + ge1.setRandom(); + ge2.setRandom(); + + VERIFY_IS_APPROX( ge3 = tri.template triangularView() * ge1, loTri * ge1); + VERIFY_IS_APPROX(rge3 = tri.template triangularView() * ge1, loTri * ge1); + VERIFY_IS_APPROX( ge3 = tri.template triangularView() * ge1, upTri * ge1); + VERIFY_IS_APPROX(rge3 = tri.template triangularView() * ge1, upTri * ge1); + VERIFY_IS_APPROX( ge3 = (s1*tri.adjoint()).template triangularView() * (s2*ge1), s1*loTri.adjoint() * (s2*ge1)); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView() * ge1, loTri.adjoint() * ge1); + VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView() * ge1, upTri.adjoint() * ge1); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView() * ge1, upTri.adjoint() * ge1); + VERIFY_IS_APPROX( ge3 = tri.template triangularView() * ge2.adjoint(), loTri * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.template triangularView() * ge2.adjoint(), loTri * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = tri.template triangularView() * ge2.adjoint(), upTri * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.template triangularView() * ge2.adjoint(), upTri * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView() * ge2.adjoint(), loTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX( ge3 = tri.adjoint().template triangularView() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); + VERIFY_IS_APPROX(rge3 = tri.adjoint().template triangularView() * ge2.adjoint(), upTri.adjoint() * ge2.adjoint()); +} + +void test_product_trmm() +{ + for(int i = 0; i < g_repeat ; i++) + { + trmm(ei_random(1,320),ei_random(1,320)); + trmm >(ei_random(1,320),ei_random(1,320)); + } +} diff --git a/test/product_triangular.cpp b/test/product_trmv.cpp similarity index 100% rename from test/product_triangular.cpp rename to test/product_trmv.cpp diff --git a/test/product_trsm.cpp b/test/product_trsm.cpp index 80df5861e..bda158048 100644 --- a/test/product_trsm.cpp +++ b/test/product_trsm.cpp @@ -85,6 +85,7 @@ template void trsm(int size,int cols) solve_ref(rmLhs.template triangularView(),rmRef); VERIFY_IS_APPROX(rmRhs, rmRef); } + void test_product_trsm() { for(int i = 0; i < g_repeat ; i++)