mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-24 14:45:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add an efficient rank2 update function (like the level2 blas xSYR2 routine).

Browse Source 
Note that it is already used in Tridiagonalization.
This commit is contained in:
Gael Guennebaud 2009-07-11 21:14:59 +02:00
parent b47dea8b7a
commit a2087cd7a3
11 changed files with 187 additions and 51 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Eigen/Core
View File

@ -180,6 +180,7 @@ namespace Eigen {
#include "src/Core/TriangularMatrix.h"
#include "src/Core/SelfAdjointView.h"
#include "src/Core/SolveTriangular.h"
#include "src/Core/products/SelfadjointRank2Update.h"
} // namespace Eigen

58
Eigen/src/Core/Product.h
View File

@ -76,7 +76,7 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
/* Helper class to analyze the factors of a Product expression.
* In particular it allows to pop out operator-, scalar multiples,
* and conjugate */
template<typename XprType> struct ei_product_factor_traits
template<typename XprType> struct ei_blas_traits
{
typedef typename ei_traits<XprType>::Scalar Scalar;
typedef XprType ActualXprType;
@ -85,15 +85,19 @@ template<typename XprType> struct ei_product_factor_traits
NeedToConjugate = false,
ActualAccess = int(ei_traits<XprType>::Flags)&DirectAccessBit ? HasDirectAccess : NoDirectAccess
};
typedef typename ei_meta_if<int(ActualAccess)==HasDirectAccess,
const ActualXprType&,
typename ActualXprType::PlainMatrixType
>::ret DirectLinearAccessType;
static inline const ActualXprType& extract(const XprType& x) { return x; }
static inline Scalar extractScalarFactor(const XprType&) { return Scalar(1); }
};
// pop conjugate
template<typename Scalar, typename NestedXpr> struct ei_product_factor_traits<CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestedXpr> >
: ei_product_factor_traits<NestedXpr>
template<typename Scalar, typename NestedXpr> struct ei_blas_traits<CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestedXpr> >
: ei_blas_traits<NestedXpr>
{
typedef ei_product_factor_traits<NestedXpr> Base;
typedef ei_blas_traits<NestedXpr> Base;
typedef CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestedXpr> XprType;
typedef typename Base::ActualXprType ActualXprType;
@ -106,10 +110,10 @@ template<typename Scalar, typename NestedXpr> struct ei_product_factor_traits<Cw
};
// pop scalar multiple
template<typename Scalar, typename NestedXpr> struct ei_product_factor_traits<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, NestedXpr> >
: ei_product_factor_traits<NestedXpr>
template<typename Scalar, typename NestedXpr> struct ei_blas_traits<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, NestedXpr> >
: ei_blas_traits<NestedXpr>
{
typedef ei_product_factor_traits<NestedXpr> Base;
typedef ei_blas_traits<NestedXpr> Base;
typedef CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, NestedXpr> XprType;
typedef typename Base::ActualXprType ActualXprType;
static inline const ActualXprType& extract(const XprType& x) { return Base::extract(x._expression()); }
@ -118,10 +122,10 @@ template<typename Scalar, typename NestedXpr> struct ei_product_factor_traits<Cw
};
// pop opposite
template<typename Scalar, typename NestedXpr> struct ei_product_factor_traits<CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, NestedXpr> >
: ei_product_factor_traits<NestedXpr>
template<typename Scalar, typename NestedXpr> struct ei_blas_traits<CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, NestedXpr> >
: ei_blas_traits<NestedXpr>
{
typedef ei_product_factor_traits<NestedXpr> Base;
typedef ei_blas_traits<NestedXpr> Base;
typedef CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, NestedXpr> XprType;
typedef typename Base::ActualXprType ActualXprType;
static inline const ActualXprType& extract(const XprType& x) { return Base::extract(x._expression()); }
@ -130,11 +134,11 @@ template<typename Scalar, typename NestedXpr> struct ei_product_factor_traits<Cw
};
// pop opposite
template<typename NestedXpr> struct ei_product_factor_traits<NestByValue<NestedXpr> >
: ei_product_factor_traits<NestedXpr>
template<typename NestedXpr> struct ei_blas_traits<NestByValue<NestedXpr> >
: ei_blas_traits<NestedXpr>
{
typedef typename NestedXpr::Scalar Scalar;
typedef ei_product_factor_traits<NestedXpr> Base;
typedef ei_blas_traits<NestedXpr> Base;
typedef NestByValue<NestedXpr> XprType;
typedef typename Base::ActualXprType ActualXprType;
static inline const ActualXprType& extract(const XprType& x) { return Base::extract(static_cast<const NestedXpr&>(x)); }
@ -148,8 +152,8 @@ template<typename NestedXpr> struct ei_product_factor_traits<NestByValue<NestedX
*/
template<typename Lhs, typename Rhs> struct ei_product_mode
{
typedef typename ei_product_factor_traits<Lhs>::ActualXprType ActualLhs;
typedef typename ei_product_factor_traits<Rhs>::ActualXprType ActualRhs;
typedef typename ei_blas_traits<Lhs>::ActualXprType ActualLhs;
typedef typename ei_blas_traits<Rhs>::ActualXprType ActualRhs;
enum{
value = Lhs::MaxColsAtCompileTime == Dynamic
@ -600,10 +604,10 @@ static void ei_cache_friendly_product_rowmajor_times_vector(
template<typename ProductType,
int LhsRows = ei_traits<ProductType>::RowsAtCompileTime,
int LhsOrder = int(ei_traits<ProductType>::LhsFlags)&RowMajorBit ? RowMajor : ColMajor,
int LhsHasDirectAccess = ei_product_factor_traits<typename ei_traits<ProductType>::_LhsNested>::ActualAccess,
int LhsHasDirectAccess = ei_blas_traits<typename ei_traits<ProductType>::_LhsNested>::ActualAccess,
int RhsCols = ei_traits<ProductType>::ColsAtCompileTime,
int RhsOrder = int(ei_traits<ProductType>::RhsFlags)&RowMajorBit ? RowMajor : ColMajor,
int RhsHasDirectAccess = ei_product_factor_traits<typename ei_traits<ProductType>::_RhsNested>::ActualAccess>
int RhsHasDirectAccess = ei_blas_traits<typename ei_traits<ProductType>::_RhsNested>::ActualAccess>
struct ei_cache_friendly_product_selector
{
template<typename DestDerived>
@ -633,8 +637,8 @@ template<typename ProductType, int LhsRows, int RhsOrder, int RhsAccess>
struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirectAccess,1,RhsOrder,RhsAccess>
{
typedef typename ProductType::Scalar Scalar;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
typedef typename RhsProductTraits::ActualXprType ActualRhsType;
@ -694,8 +698,8 @@ template<typename ProductType, int LhsOrder, int LhsAccess, int RhsCols>
struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCols,RowMajor,HasDirectAccess>
{
typedef typename ProductType::Scalar Scalar;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
typedef typename RhsProductTraits::ActualXprType ActualRhsType;
@ -740,8 +744,8 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirect
{
typedef typename ProductType::Scalar Scalar;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
typedef typename RhsProductTraits::ActualXprType ActualRhsType;
@ -783,8 +787,8 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
{
typedef typename ProductType::Scalar Scalar;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_product_factor_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_LhsNested> LhsProductTraits;
typedef ei_blas_traits<typename ei_traits<ProductType>::_RhsNested> RhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
typedef typename RhsProductTraits::ActualXprType ActualRhsType;
@ -903,8 +907,8 @@ template<typename Lhs, typename Rhs, int ProductMode>
template<typename DestDerived>
inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived& res, Scalar alpha) const
{
typedef ei_product_factor_traits<_LhsNested> LhsProductTraits;
typedef ei_product_factor_traits<_RhsNested> RhsProductTraits;
typedef ei_blas_traits<_LhsNested> LhsProductTraits;
typedef ei_blas_traits<_RhsNested> RhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
typedef typename RhsProductTraits::ActualXprType ActualRhsType;

10
Eigen/src/Core/SelfAdjointView.h
View File

@ -106,6 +106,16 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
return ei_selfadjoint_vector_product_returntype<SelfAdjointView,OtherDerived>(*this, rhs.derived());
}
/** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
* \f$ this = this + \alpha ( u v^* + v u^*) \f$
*
* The vectors \a u and \c v \b must be column vectors, however they can be
* a adjoint expression without any overhead. Only the meaningful triangular
* part of the matrix is updated, the rest is left unchanged.
*/
template<typename DerivedU, typename DerivedV>
void rank2update(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
/////////// Cholesky module ///////////
const LLT<PlainMatrixType, UpLo> llt() const;

17
Eigen/src/Core/SolveTriangular.h
View File

@ -33,12 +33,12 @@ template<typename Lhs, typename Rhs,
>
struct ei_triangular_solver_selector;
// forward substitution, row-major
// forward and backward substitution, row-major
template<typename Lhs, typename Rhs, int Mode>
struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,RowMajor>
{
typedef typename Rhs::Scalar Scalar;
typedef ei_product_factor_traits<Lhs> LhsProductTraits;
typedef ei_blas_traits<Lhs> LhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
enum {
IsLowerTriangular = ((Mode&LowerTriangularBit)==LowerTriangularBit)
@ -60,6 +60,9 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,RowMajor>
int r = IsLowerTriangular ? pi : size - pi; // remaining size
if (r > 0)
{
// let's directly call the low level product function because:
// 1 - it is faster to compile
// 2 - it is slighlty faster at runtime
int startRow = IsLowerTriangular ? pi : pi-actualPanelWidth;
int startCol = IsLowerTriangular ? 0 : pi;
Block<Rhs,Dynamic,1> target(other,startRow,c,actualPanelWidth,1);
@ -86,17 +89,13 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,RowMajor>
}
};
// Implements the following configurations:
// - inv(LowerTriangular, ColMajor) * Column vectors
// - inv(LowerTriangular,UnitDiag,ColMajor) * Column vectors
// - inv(UpperTriangular, ColMajor) * Column vectors
// - inv(UpperTriangular,UnitDiag,ColMajor) * Column vectors
// forward and backward substitution, column-major
template<typename Lhs, typename Rhs, int Mode>
struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,ColMajor>
{
typedef typename Rhs::Scalar Scalar;
typedef typename ei_packet_traits<Scalar>::type Packet;
typedef ei_product_factor_traits<Lhs> LhsProductTraits;
typedef ei_blas_traits<Lhs> LhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
enum {
PacketSize = ei_packet_traits<Scalar>::size,
@ -136,7 +135,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,NoUnrolling,ColMajor>
int r = IsLowerTriangular ? size - endBlock : startBlock; // remaining size
if (r > 0)
{
// let's directly call this function because:
// let's directly call the low level product function because:
// 1 - it is faster to compile
// 2 - it is slighlty faster at runtime
ei_cache_friendly_product_colmajor_times_vector<LhsProductTraits::NeedToConjugate,false>(

2
Eigen/src/Core/products/GeneralMatrixMatrix.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-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

2
Eigen/src/Core/products/GeneralMatrixVector.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-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

2
Eigen/src/Core/products/SelfadjointMatrixVector.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-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

96
Eigen/src/Core/products/SelfadjointRank2Update.h Normal file
View File

@ -0,0 +1,96 @@
// 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_SELFADJOINTRANK2UPTADE_H
#define EIGEN_SELFADJOINTRANK2UPTADE_H
/* Optimized selfadjoint matrix += alpha * uv' + vu'
* It corresponds to the Level2 syr2 BLAS routine
*/
template<typename Scalar, typename UType, typename VType, int UpLo>
struct ei_selfadjoint_rank2_update_selector;
template<typename Scalar, typename UType, typename VType>
struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,LowerTriangular>
{
static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha)
{
const int size = u.size();
// std::cerr << "lower \n" << u.transpose() << "\n" << v.transpose() << "\n\n";
for (int i=0; i<size; ++i)
{
// std::cerr <<
Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i+i, size-i) +=
(alpha * ei_conj(u.coeff(i))) * v.end(size-i)
+ (alpha * ei_conj(v.coeff(i))) * u.end(size-i);
}
}
};
template<typename Scalar, typename UType, typename VType>
struct ei_selfadjoint_rank2_update_selector<Scalar,UType,VType,UpperTriangular>
{
static void run(Scalar* mat, int stride, const UType& u, const VType& v, Scalar alpha)
{
const int size = u.size();
for (int i=0; i<size; ++i)
Map<Matrix<Scalar,Dynamic,1> >(mat+stride*i, i+1) +=
(alpha * ei_conj(u.coeff(i))) * v.start(i+1)
+ (alpha * ei_conj(v.coeff(i))) * u.start(i+1);
}
};
template<bool Cond, typename T> struct ei_conj_expr_if
: ei_meta_if<!Cond, T,
CwiseUnaryOp<ei_scalar_conjugate_op<typename ei_traits<T>::Scalar>,T> > {};
template<typename MatrixType, unsigned int UpLo>
template<typename DerivedU, typename DerivedV>
void SelfAdjointView<MatrixType,UpLo>
::rank2update(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha)
{
typedef ei_blas_traits<DerivedU> UBlasTraits;
typedef typename UBlasTraits::DirectLinearAccessType ActualUType;
typedef typename ei_cleantype<ActualUType>::type _ActualUType;
const ActualUType actualU = UBlasTraits::extract(u.derived());
typedef ei_blas_traits<DerivedV> VBlasTraits;
typedef typename VBlasTraits::DirectLinearAccessType ActualVType;
typedef typename ei_cleantype<ActualVType>::type _ActualVType;
const ActualVType actualV = VBlasTraits::extract(v.derived());
Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived())
* VBlasTraits::extractScalarFactor(v.derived());
enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit)?1:0 };
ei_selfadjoint_rank2_update_selector<Scalar,
typename ei_conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::ret,
typename ei_conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::ret,
(IsRowMajor ? (UpLo==UpperTriangular ? LowerTriangular : UpperTriangular) : UpLo)>
::run(const_cast<Scalar*>(_expression().data()),_expression().stride(),actualU,actualV,actualAlpha);
}
#endif // EIGEN_SELFADJOINTRANK2UPTADE_H

6
Eigen/src/QR/Tridiagonalization.h
View File

@ -236,10 +236,8 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
+ (h*ei_conj(h)*Scalar(-0.5)*(matA.col(i).end(n-i-1).dot(hCoeffs.end(n-i-1)))) *
matA.col(i).end(n-i-1);
// symmetric rank-2 update
for (int j1=i+1; j1<n; ++j1)
matA.col(j1).end(n-j1) -= matA.col(i).end(n-j1) * ei_conj(hCoeffs.coeff(j1-1))
+ hCoeffs.end(n-j1) * ei_conj(matA.coeff(j1,i));
matA.corner(BottomRight, n-i-1, n-i-1).template selfadjointView<LowerTriangular>()
.rank2update(matA.col(i).end(n-i-1), hCoeffs.end(n-i-1), -1);
// note: at that point matA(i+1,i+1) is the (i+1)-th element of the final diagonal
// note: the sequence of the beta values leads to the subdiagonal entries

4
test/eigensolver_selfadjoint.cpp
View File

@ -119,8 +119,8 @@ void test_eigensolver_selfadjoint()
// very important to test a 3x3 matrix since we provide a special path for it
CALL_SUBTEST( selfadjointeigensolver(Matrix3f()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix4d()) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXf(7,7)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(5,5)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXf(4,4)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(7,7)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
// some trivial but implementation-wise tricky cases

40
test/product_selfadjoint.cpp
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@gmail.com>
// 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
@ -29,20 +29,29 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
typedef Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> RowVectorType;
int rows = m.rows();
int cols = m.cols();
MatrixType m1 = MatrixType::Random(rows, cols),
m2 = MatrixType::Random(rows, cols);
m2 = MatrixType::Random(rows, cols),
m3;
VectorType v1 = VectorType::Random(rows),
v2 = VectorType::Random(rows);
RowVectorType r1 = RowVectorType::Random(rows),
r2 = RowVectorType::Random(rows);
Scalar s1 = ei_random<Scalar>(),
s2 = ei_random<Scalar>(),
s3 = ei_random<Scalar>();
m1 = m1.adjoint()*m1;
// lower
m2.setZero();
m2.template part<LowerTriangular>() = m1;
m2.template triangularView<LowerTriangular>() = m1;
ei_product_selfadjoint_vector<Scalar,MatrixType::Flags&RowMajorBit,LowerTriangularBit>
(cols,m2.data(),cols, v1.data(), v2.data());
VERIFY_IS_APPROX(v2, m1 * v1);
@ -50,11 +59,30 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
// upper
m2.setZero();
m2.template part<UpperTriangular>() = m1;
m2.template triangularView<UpperTriangular>() = m1;
ei_product_selfadjoint_vector<Scalar,MatrixType::Flags&RowMajorBit,UpperTriangularBit>(cols,m2.data(),cols, v1.data(), v2.data());
VERIFY_IS_APPROX(v2, m1 * v1);
VERIFY_IS_APPROX((m2.template selfadjointView<UpperTriangular>() * v1).eval(), m1 * v1);
// rank2 update
m2 = m1.template triangularView<LowerTriangular>();
m2.template selfadjointView<LowerTriangular>().rank2update(v1,v2);
VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<LowerTriangular>().toDense());
m2 = m1.template triangularView<UpperTriangular>();
m2.template selfadjointView<UpperTriangular>().rank2update(-v1,s2*v2,s3);
VERIFY_IS_APPROX(m2, (m1 + (-s2*s3) * (v1 * v2.adjoint()+ v2 * v1.adjoint())).template triangularView<UpperTriangular>().toDense());
m2 = m1.template triangularView<UpperTriangular>();
m2.template selfadjointView<UpperTriangular>().rank2update(-r1.adjoint(),r2.adjoint()*s3,s1);
VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<UpperTriangular>().toDense());
m2 = m1.template triangularView<LowerTriangular>();
m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rank2update(v1.end(rows-1),v2.start(cols-1));
m3 = m1;
m3.block(1,1,rows-1,cols-1) += v1.end(rows-1) * v2.start(cols-1).adjoint()+ v2.start(cols-1) * v1.end(rows-1).adjoint();
VERIFY_IS_APPROX(m2, m3.template triangularView<LowerTriangular>().toDense());
}
void test_product_selfadjoint()
@ -65,8 +93,8 @@ void test_product_selfadjoint()
CALL_SUBTEST( product_selfadjoint(Matrix3d()) );
CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) );
CALL_SUBTEST( product_selfadjoint(MatrixXcd(21,21)) );
CALL_SUBTEST( product_selfadjoint(MatrixXd(17,17)) );
CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(18,18)) );
CALL_SUBTEST( product_selfadjoint(MatrixXd(4,4)) );
CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(17,17)) );
CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
}
}