mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-11-27 06:30:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix selfadjoint to dense

Browse Source
This commit is contained in:
Gael Guennebaud 2010-05-19 16:35:34 +02:00
parent 08fbfa93e0
commit bf09fe55ed
4 changed files with 129 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -39,14 +39,14 @@
*
* \sa class TriangularBase, MatrixBase::selfAdjointView()
*/
template<typename MatrixType, unsigned int TriangularPart>
struct ei_traits<SelfAdjointView<MatrixType, TriangularPart> > : ei_traits<MatrixType>
template<typename MatrixType, unsigned int UpLo>
struct ei_traits<SelfAdjointView<MatrixType, UpLo> > : ei_traits<MatrixType>
{
typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
typedef MatrixType ExpressionType;
enum {
Mode = TriangularPart | SelfAdjoint,
Mode = UpLo | SelfAdjoint,
Flags = _MatrixTypeNested::Flags & (HereditaryBits)
& (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved
CoeffReadCost = _MatrixTypeNested::CoeffReadCost
@ -77,7 +77,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
inline int cols() const { return m_matrix.cols(); }
inline int outerStride() const { return m_matrix.outerStride(); }
inline int innerStride() const { return m_matrix.innerStride(); }
/** \sa MatrixBase::coeff()
* \warning the coordinates must fit into the referenced triangular part
*/
@ -165,8 +165,10 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
// return ei_matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
// }
// selfadjoint to dense matrix
template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, UnrollCount, ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
{
enum {
col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
@ -175,7 +177,7 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Unroll
inline static void run(Derived1 &dst, const Derived2 &src)
{
ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, UnrollCount-1, ClearOpposite>::run(dst, src);
ei_triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src);
if(row == col)
dst.coeffRef(row, col) = ei_real(src.coeff(row, col));
@ -184,17 +186,67 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Unroll
}
};
// selfadjoint to dense matrix
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Dynamic, ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite>
{
inline static void run(Derived1 &, const Derived2 &) {}
};
template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite>
{
enum {
col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
row = (UnrollCount-1) % Derived1::RowsAtCompileTime
};
inline static void run(Derived1 &dst, const Derived2 &src)
{
ei_triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src);
if(row == col)
dst.coeffRef(row, col) = ei_real(src.coeff(row, col));
else if(row > col)
dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col));
}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite>
{
inline static void run(Derived1 &, const Derived2 &) {}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
for(int j = 0; j < dst.cols(); ++j)
{
for(int i = 0; i < j; ++i)
dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j));
dst.coeffRef(j, j) = ei_real(src.coeff(j, j));
{
dst.copyCoeff(i, j, src);
dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j));
}
dst.copyCoeff(j, j, src);
}
}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
for(int i = 0; i < dst.rows(); ++i)
{
for(int j = 0; j < i; ++j)
{
dst.copyCoeff(i, j, src);
dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j));
}
dst.copyCoeff(i, i, src);
}
}
};

34
Eigen/src/Core/TriangularMatrix.h
View File

@ -472,40 +472,6 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic,
}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
for(int j = 0; j < dst.cols(); ++j)
{
for(int i = 0; i < j; ++i)
{
dst.copyCoeff(i, j, src);
dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j));
}
dst.copyCoeff(j, j, src);
}
}
};
template<typename Derived1, typename Derived2, bool ClearOpposite>
struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
{
inline static void run(Derived1 &dst, const Derived2 &src)
{
for(int i = 0; i < dst.rows(); ++i)
{
for(int j = 0; j < i; ++j)
{
dst.copyCoeff(i, j, src);
dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j));
}
dst.copyCoeff(i, i, src);
}
}
};
// FIXME should we keep that possibility
template<typename MatrixType, unsigned int Mode>
template<typename OtherDerived>

1
test/CMakeLists.txt
View File

@ -116,6 +116,7 @@ ei_add_test(array_for_matrix)
ei_add_test(array_replicate)
ei_add_test(array_reverse)
ei_add_test(triangular)
ei_add_test(selfadjoint)
ei_add_test(product_selfadjoint)
ei_add_test(product_symm)
ei_add_test(product_syrk)

66
test/selfadjoint.cpp Normal file
View File

@ -0,0 +1,66 @@
// This file is triangularView of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 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"
// This file tests the basic selfadjointView API,
// the related products and decompositions are tested in specific files.
template<typename MatrixType> void selfadjoint(const MatrixType& m)
{
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
int rows = m.rows();
int cols = m.cols();
MatrixType m1 = MatrixType::Random(rows, cols),
m3(rows, cols);
m1.diagonal() = m1.diagonal().real().template cast<Scalar>();
// check selfadjoint to dense
m3 = m1.template selfadjointView<Upper>();
VERIFY_IS_APPROX(MatrixType(m3.template triangularView<Upper>()), MatrixType(m1.template triangularView<Upper>()));
VERIFY_IS_APPROX(m3, m3.adjoint());
m3 = m1.template selfadjointView<Lower>();
VERIFY_IS_APPROX(MatrixType(m3.template triangularView<Lower>()), MatrixType(m1.template triangularView<Lower>()));
VERIFY_IS_APPROX(m3, m3.adjoint());
}
void test_selfadjoint()
{
for(int i = 0; i < g_repeat ; i++)
{
EIGEN_UNUSED int s = ei_random<int>(1,20);
CALL_SUBTEST_1( selfadjoint(Matrix<float, 1, 1>()) );
CALL_SUBTEST_2( selfadjoint(Matrix<float, 2, 2>()) );
CALL_SUBTEST_3( selfadjoint(Matrix3cf()) );
CALL_SUBTEST_4( selfadjoint(MatrixXcd(s,s)) );
CALL_SUBTEST_5( selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(s, s)) );
}
}