mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-19 16:19:37 +08:00
add triangular * vector product
This commit is contained in:
Gael Guennebaud
parent
a2087cd7a3
commit
a2cf7ba955
@ -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
|
||||
|
||||
|
||||
@ -45,7 +45,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,RowMajor>
|
||||
};
|
||||
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<Lhs,Rhs,Mode,NoUnrolling,ColMajor>
|
||||
|
||||
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();
|
||||
|
||||
@ -142,13 +142,17 @@ struct ei_traits<TriangularView<MatrixType, _Mode> > : ei_traits<MatrixType>
|
||||
};
|
||||
};
|
||||
|
||||
template<typename MatrixType, unsigned int _Mode> class TriangularView
|
||||
: public TriangularBase<TriangularView<MatrixType, _Mode> >
|
||||
template<typename Lhs,typename Rhs>
|
||||
struct ei_triangular_vector_product_returntype;
|
||||
|
||||
template<typename _MatrixType, unsigned int _Mode> class TriangularView
|
||||
: public TriangularBase<TriangularView<_MatrixType, _Mode> >
|
||||
{
|
||||
public:
|
||||
|
||||
typedef TriangularBase<TriangularView> Base;
|
||||
typedef typename ei_traits<TriangularView>::Scalar Scalar;
|
||||
typedef _MatrixType MatrixType;
|
||||
typedef typename MatrixType::PlainMatrixType PlainMatrixType;
|
||||
|
||||
enum {
|
||||
@ -243,6 +247,13 @@ template<typename MatrixType, unsigned int _Mode> class TriangularView
|
||||
return res;
|
||||
}
|
||||
|
||||
template<typename OtherDerived>
|
||||
ei_triangular_vector_product_returntype<TriangularView,OtherDerived>
|
||||
operator*(const MatrixBase<OtherDerived>& rhs) const
|
||||
{
|
||||
return ei_triangular_vector_product_returntype<TriangularView,OtherDerived>(*this, rhs.derived(), 1);
|
||||
}
|
||||
|
||||
template<typename OtherDerived>
|
||||
typename ei_plain_matrix_type_column_major<OtherDerived>::type
|
||||
solve(const MatrixBase<OtherDerived>& other) const;
|
||||
|
||||
@ -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.
|
||||
|
||||
161
Eigen/src/Core/products/TriangularMatrixVector.h
Normal file
161
Eigen/src/Core/products/TriangularMatrixVector.h
Normal file
@ -0,0 +1,161 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
|
||||
//
|
||||
// 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
#ifndef EIGEN_TRIANGULARMATRIXVECTOR_H
|
||||
#define EIGEN_TRIANGULARMATRIXVECTOR_H
|
||||
|
||||
template<typename MatrixType, typename Rhs, typename Result,
|
||||
int Mode, bool ConjLhs, bool ConjRhs, int StorageOrder>
|
||||
struct ei_product_triangular_vector_selector;
|
||||
|
||||
template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs>
|
||||
struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,ColMajor>
|
||||
{
|
||||
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<Lhs>::Scalar alpha)
|
||||
{
|
||||
static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
|
||||
typename ei_conj_expr_if<ConjLhs,Lhs>::ret cjLhs(lhs);
|
||||
typename ei_conj_expr_if<ConjRhs,Rhs>::ret cjRhs(rhs);
|
||||
|
||||
int size = lhs.cols();
|
||||
for (int pi=0; pi<size; pi+=PanelWidth)
|
||||
{
|
||||
int actualPanelWidth = std::min(PanelWidth, size-pi);
|
||||
for (int k=0; k<actualPanelWidth; ++k)
|
||||
{
|
||||
int i = pi + k;
|
||||
int s = IsLowerTriangular ? (HasUnitDiag ? i+1 : i ) : pi;
|
||||
int r = IsLowerTriangular ? actualPanelWidth-k : k+1;
|
||||
if ((!HasUnitDiag) || (--r)>0)
|
||||
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<ConjLhs,ConjRhs>(
|
||||
r,
|
||||
&(lhs.const_cast_derived().coeffRef(s,pi)), lhs.stride(),
|
||||
rhs.segment(pi, actualPanelWidth),
|
||||
&(res.coeffRef(s)),
|
||||
alpha);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs>
|
||||
struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,RowMajor>
|
||||
{
|
||||
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<Lhs>::Scalar alpha)
|
||||
{
|
||||
static const int PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
|
||||
typename ei_conj_expr_if<ConjLhs,Lhs>::ret cjLhs(lhs);
|
||||
typename ei_conj_expr_if<ConjRhs,Rhs>::ret cjRhs(rhs);
|
||||
int size = lhs.cols();
|
||||
for (int pi=0; pi<size; pi+=PanelWidth)
|
||||
{
|
||||
int actualPanelWidth = std::min(PanelWidth, size-pi);
|
||||
for (int k=0; k<actualPanelWidth; ++k)
|
||||
{
|
||||
int i = pi + k;
|
||||
int s = IsLowerTriangular ? pi : (HasUnitDiag ? i+1 : i);
|
||||
int r = IsLowerTriangular ? k+1 : actualPanelWidth-k;
|
||||
if ((!HasUnitDiag) || (--r)>0)
|
||||
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<Result,Dynamic,1> target(res,pi,0,actualPanelWidth,1);
|
||||
ei_cache_friendly_product_rowmajor_times_vector<ConjLhs,ConjRhs>(
|
||||
&(lhs.const_cast_derived().coeffRef(pi,s)), lhs.stride(),
|
||||
&(rhs.const_cast_derived().coeffRef(s)), r,
|
||||
target, alpha);
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs,typename Rhs>
|
||||
struct ei_triangular_vector_product_returntype
|
||||
: public ReturnByValue<ei_triangular_vector_product_returntype<Lhs,Rhs>,
|
||||
Matrix<typename ei_traits<Rhs>::Scalar,
|
||||
Rhs::RowsAtCompileTime,Rhs::ColsAtCompileTime> >
|
||||
{
|
||||
typedef typename Lhs::Scalar Scalar;
|
||||
typedef typename ei_cleantype<typename Rhs::Nested>::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<typename Dest> void evalTo(Dest& dst) const
|
||||
{
|
||||
typedef typename Lhs::MatrixType MatrixType;
|
||||
|
||||
typedef ei_blas_traits<MatrixType> LhsBlasTraits;
|
||||
typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
|
||||
typedef typename ei_cleantype<ActualLhsType>::type _ActualLhsType;
|
||||
const ActualLhsType actualLhs = LhsBlasTraits::extract(m_lhs._expression());
|
||||
|
||||
typedef ei_blas_traits<Rhs> RhsBlasTraits;
|
||||
typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhsType;
|
||||
typedef typename ei_cleantype<ActualRhsType>::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<Lhs>::Mode,
|
||||
LhsBlasTraits::NeedToConjugate,
|
||||
RhsBlasTraits::NeedToConjugate,
|
||||
ei_traits<Lhs>::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
|
||||
@ -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.
|
||||
|
||||
@ -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)
|
||||
|
||||
92
test/product_triangular.cpp
Normal file
92
test/product_triangular.cpp
Normal file
@ -0,0 +1,92 @@
|
||||
// This file is triangularView of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@gmail.com>
|
||||
//
|
||||
// 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
#include "main.h"
|
||||
|
||||
template<typename MatrixType> void product_triangular(const MatrixType& m)
|
||||
{
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename NumTraits<Scalar>::Real RealScalar;
|
||||
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
|
||||
|
||||
RealScalar largerEps = 10*test_precision<RealScalar>();
|
||||
|
||||
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<Scalar>();
|
||||
|
||||
m1 = MatrixType::Random(rows, cols);
|
||||
|
||||
// check with a column-major matrix
|
||||
m3 = m1.template triangularView<Eigen::LowerTriangular>();
|
||||
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UpperTriangular>();
|
||||
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UpperTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UnitLowerTriangular>();
|
||||
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
|
||||
VERIFY((m3 * v1).isApprox(m1.template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps));
|
||||
|
||||
// check conjugated and scalar multiple expressions (col-major)
|
||||
m3 = m1.template triangularView<Eigen::LowerTriangular>();
|
||||
VERIFY(((s1*m3).conjugate() * v1).isApprox((s1*m1).conjugate().template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UpperTriangular>();
|
||||
VERIFY((m3.conjugate() * v1.conjugate()).isApprox(m1.conjugate().template triangularView<Eigen::UpperTriangular>() * v1.conjugate(), largerEps));
|
||||
|
||||
// check with a row-major matrix
|
||||
m3 = m1.template triangularView<Eigen::UpperTriangular>();
|
||||
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::LowerTriangular>();
|
||||
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UpperTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
|
||||
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitLowerTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UnitLowerTriangular>();
|
||||
VERIFY((m3.transpose() * v1).isApprox(m1.transpose().template triangularView<Eigen::UnitUpperTriangular>() * v1, largerEps));
|
||||
|
||||
// check conjugated and scalar multiple expressions (row-major)
|
||||
m3 = m1.template triangularView<Eigen::UpperTriangular>();
|
||||
VERIFY((m3.adjoint() * v1).isApprox(m1.adjoint().template triangularView<Eigen::LowerTriangular>() * v1, largerEps));
|
||||
m3 = m1.template triangularView<Eigen::LowerTriangular>();
|
||||
VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::UpperTriangular>() * (s1*v1.conjugate()), largerEps));
|
||||
m3 = m1.template triangularView<Eigen::UnitUpperTriangular>();
|
||||
|
||||
// TODO check with sub-matrices
|
||||
}
|
||||
|
||||
void test_product_triangular()
|
||||
{
|
||||
for(int i = 0; i < g_repeat ; i++) {
|
||||
CALL_SUBTEST( product_triangular(Matrix<float, 1, 1>()) );
|
||||
CALL_SUBTEST( product_triangular(Matrix<float, 2, 2>()) );
|
||||
CALL_SUBTEST( product_triangular(Matrix3d()) );
|
||||
CALL_SUBTEST( product_triangular(Matrix<std::complex<float>,23, 23>()) );
|
||||
CALL_SUBTEST( product_triangular(MatrixXcd(17,17)) );
|
||||
CALL_SUBTEST( product_triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(19, 19)) );
|
||||
}
|
||||
}
|
||||
@ -51,6 +51,8 @@ template<typename MatrixType> void triangular(const MatrixType& m)
|
||||
v2 = VectorType::Random(rows),
|
||||
vzero = VectorType::Zero(rows);
|
||||
|
||||
Scalar s1 = ei_random<Scalar>();
|
||||
|
||||
MatrixType m1up = m1.template triangularView<Eigen::UpperTriangular>();
|
||||
MatrixType m2up = m2.template triangularView<Eigen::UpperTriangular>();
|
||||
|
||||
@ -142,7 +144,7 @@ void test_triangular()
|
||||
CALL_SUBTEST( triangular(Matrix3d()) );
|
||||
CALL_SUBTEST( triangular(MatrixXcf(4, 4)) );
|
||||
CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) );
|
||||
CALL_SUBTEST( triangular(MatrixXd(17,17)) );
|
||||
CALL_SUBTEST( triangular(MatrixXcd(17,17)) );
|
||||
CALL_SUBTEST( triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(5, 5)) );
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user