Rasmus Munk LarsenandRasmus Munk Larsen
8a3d61c790
libeigen/eigen!2509 Closes #2868 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
539 lines
24 KiB
C++
539 lines
24 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla
|
|
// Public License v. 2.0. If a copy of the MPL was not distributed
|
|
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
// SPDX-License-Identifier: MPL-2.0
|
|
|
|
// discard stack allocation as that too bypasses malloc
|
|
#define EIGEN_STACK_ALLOCATION_LIMIT 0
|
|
// heap allocation will raise an assert if enabled at runtime
|
|
#define EIGEN_RUNTIME_NO_MALLOC
|
|
|
|
#include "main.h"
|
|
using namespace std;
|
|
template <typename MatrixType>
|
|
void diagonalmatrices(const MatrixType& m) {
|
|
typedef typename MatrixType::Scalar Scalar;
|
|
enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
|
|
typedef Matrix<Scalar, Rows, 1> VectorType;
|
|
typedef Matrix<Scalar, 1, Cols> RowVectorType;
|
|
typedef Matrix<Scalar, Rows, Rows> SquareMatrixType;
|
|
typedef Matrix<Scalar, Dynamic, Dynamic> DynMatrixType;
|
|
typedef DiagonalMatrix<Scalar, Rows> LeftDiagonalMatrix;
|
|
typedef DiagonalMatrix<Scalar, Cols> RightDiagonalMatrix;
|
|
typedef Matrix<Scalar, Rows == Dynamic ? Dynamic : 2 * Rows, Cols == Dynamic ? Dynamic : 2 * Cols> BigMatrix;
|
|
Index rows = m.rows();
|
|
Index cols = m.cols();
|
|
|
|
MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols);
|
|
VectorType v1 = VectorType::Random(rows), v2 = VectorType::Random(rows);
|
|
RowVectorType rv1 = RowVectorType::Random(cols), rv2 = RowVectorType::Random(cols);
|
|
|
|
LeftDiagonalMatrix ldm1(v1), ldm2(v2);
|
|
RightDiagonalMatrix rdm1(rv1), rdm2(rv2);
|
|
|
|
Scalar s1 = internal::random<Scalar>();
|
|
|
|
SquareMatrixType sq_m1(v1.asDiagonal());
|
|
VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
|
|
sq_m1 = v1.asDiagonal();
|
|
VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
|
|
SquareMatrixType sq_m2 = v1.asDiagonal();
|
|
VERIFY_IS_APPROX(sq_m1, sq_m2);
|
|
|
|
ldm1 = v1.asDiagonal();
|
|
LeftDiagonalMatrix ldm3(v1);
|
|
VERIFY_IS_APPROX(ldm1.diagonal(), ldm3.diagonal());
|
|
LeftDiagonalMatrix ldm4 = v1.asDiagonal();
|
|
VERIFY_IS_APPROX(ldm1.diagonal(), ldm4.diagonal());
|
|
|
|
sq_m1.block(0, 0, rows, rows) = ldm1;
|
|
VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
|
|
sq_m1.transpose() = ldm1;
|
|
VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
|
|
|
|
Index i = internal::random<Index>(0, rows - 1);
|
|
Index j = internal::random<Index>(0, cols - 1);
|
|
|
|
internal::set_is_malloc_allowed(false);
|
|
VERIFY_IS_APPROX(((ldm1 * m1)(i, j)), ldm1.diagonal()(i) * m1(i, j));
|
|
VERIFY_IS_APPROX(((ldm1 * (m1 + m2))(i, j)), ldm1.diagonal()(i) * (m1 + m2)(i, j));
|
|
VERIFY_IS_APPROX(((m1 * rdm1)(i, j)), rdm1.diagonal()(j) * m1(i, j));
|
|
VERIFY_IS_APPROX(((v1.asDiagonal() * m1)(i, j)), v1(i) * m1(i, j));
|
|
VERIFY_IS_APPROX(((m1 * rv1.asDiagonal())(i, j)), rv1(j) * m1(i, j));
|
|
VERIFY_IS_APPROX((((v1 + v2).asDiagonal() * m1)(i, j)), (v1 + v2)(i)*m1(i, j));
|
|
VERIFY_IS_APPROX((((v1 + v2).asDiagonal() * (m1 + m2))(i, j)), (v1 + v2)(i) * (m1 + m2)(i, j));
|
|
VERIFY_IS_APPROX(((m1 * (rv1 + rv2).asDiagonal())(i, j)), (rv1 + rv2)(j)*m1(i, j));
|
|
VERIFY_IS_APPROX((((m1 + m2) * (rv1 + rv2).asDiagonal())(i, j)), (rv1 + rv2)(j) * (m1 + m2)(i, j));
|
|
VERIFY_IS_APPROX((ldm1 * ldm1).diagonal()(i), ldm1.diagonal()(i) * ldm1.diagonal()(i));
|
|
VERIFY_IS_APPROX((ldm1 * ldm1 * m1)(i, j), ldm1.diagonal()(i) * ldm1.diagonal()(i) * m1(i, j));
|
|
VERIFY_IS_APPROX(((v1.asDiagonal() * v1.asDiagonal()).diagonal()(i)), v1(i) * v1(i));
|
|
internal::set_is_malloc_allowed(true);
|
|
|
|
if (rows > 1) {
|
|
DynMatrixType tmp = m1.topRows(rows / 2), res;
|
|
VERIFY_IS_APPROX((res = m1.topRows(rows / 2) * rv1.asDiagonal()), tmp * rv1.asDiagonal());
|
|
VERIFY_IS_APPROX((res = v1.head(rows / 2).asDiagonal() * m1.topRows(rows / 2)),
|
|
v1.head(rows / 2).asDiagonal() * tmp);
|
|
}
|
|
|
|
BigMatrix big;
|
|
big.setZero(2 * rows, 2 * cols);
|
|
|
|
big.block(i, j, rows, cols) = m1;
|
|
big.block(i, j, rows, cols) = v1.asDiagonal() * big.block(i, j, rows, cols);
|
|
|
|
VERIFY_IS_APPROX((big.block(i, j, rows, cols)), v1.asDiagonal() * m1);
|
|
|
|
big.block(i, j, rows, cols) = m1;
|
|
big.block(i, j, rows, cols) = big.block(i, j, rows, cols) * rv1.asDiagonal();
|
|
VERIFY_IS_APPROX((big.block(i, j, rows, cols)), m1 * rv1.asDiagonal());
|
|
|
|
// products do not allocate memory
|
|
MatrixType res(rows, cols);
|
|
internal::set_is_malloc_allowed(false);
|
|
res.noalias() = ldm1 * m1;
|
|
res.noalias() = m1 * rdm1;
|
|
res.noalias() = ldm1 * m1 * rdm1;
|
|
res.noalias() = LeftDiagonalMatrix::Identity(rows) * m1 * RightDiagonalMatrix::Zero(cols);
|
|
internal::set_is_malloc_allowed(true);
|
|
|
|
// scalar multiple
|
|
VERIFY_IS_APPROX(LeftDiagonalMatrix(ldm1 * s1).diagonal(), ldm1.diagonal() * s1);
|
|
VERIFY_IS_APPROX(LeftDiagonalMatrix(s1 * ldm1).diagonal(), s1 * ldm1.diagonal());
|
|
|
|
VERIFY_IS_APPROX(m1 * (rdm1 * s1), (m1 * rdm1) * s1);
|
|
VERIFY_IS_APPROX(m1 * (s1 * rdm1), (m1 * rdm1) * s1);
|
|
|
|
// Diagonal to dense
|
|
sq_m1.setRandom();
|
|
sq_m2 = sq_m1;
|
|
VERIFY_IS_APPROX((sq_m1 += (s1 * v1).asDiagonal()), sq_m2 += (s1 * v1).asDiagonal().toDenseMatrix());
|
|
VERIFY_IS_APPROX((sq_m1 -= (s1 * v1).asDiagonal()), sq_m2 -= (s1 * v1).asDiagonal().toDenseMatrix());
|
|
VERIFY_IS_APPROX((sq_m1 = (s1 * v1).asDiagonal()), (s1 * v1).asDiagonal().toDenseMatrix());
|
|
|
|
sq_m1.setRandom();
|
|
sq_m2 = v1.asDiagonal();
|
|
sq_m2 = sq_m1 * sq_m2;
|
|
VERIFY_IS_APPROX((sq_m1 * v1.asDiagonal()).col(i), sq_m2.col(i));
|
|
VERIFY_IS_APPROX((sq_m1 * v1.asDiagonal()).row(i), sq_m2.row(i));
|
|
|
|
sq_m1 = v1.asDiagonal();
|
|
sq_m2 = v2.asDiagonal();
|
|
SquareMatrixType sq_m3 = v1.asDiagonal();
|
|
VERIFY_IS_APPROX(sq_m3 = v1.asDiagonal() + v2.asDiagonal(), sq_m1 + sq_m2);
|
|
VERIFY_IS_APPROX(sq_m3 = v1.asDiagonal() - v2.asDiagonal(), sq_m1 - sq_m2);
|
|
VERIFY_IS_APPROX(sq_m3 = v1.asDiagonal() - 2 * v2.asDiagonal() + v1.asDiagonal(), sq_m1 - 2 * sq_m2 + sq_m1);
|
|
|
|
// Zero and Identity
|
|
LeftDiagonalMatrix zero = LeftDiagonalMatrix::Zero(rows);
|
|
LeftDiagonalMatrix identity = LeftDiagonalMatrix::Identity(rows);
|
|
VERIFY_IS_APPROX(identity.diagonal().sum(), Scalar(rows));
|
|
VERIFY_IS_APPROX(zero.diagonal().sum(), Scalar(0));
|
|
VERIFY_IS_APPROX((zero + 2 * LeftDiagonalMatrix::Identity(rows)).diagonal().sum(), Scalar(2 * rows));
|
|
}
|
|
|
|
template <typename MatrixType>
|
|
void as_scalar_product(const MatrixType& m) {
|
|
typedef typename MatrixType::Scalar Scalar;
|
|
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
|
|
typedef Matrix<Scalar, Dynamic, Dynamic> DynMatrixType;
|
|
typedef Matrix<Scalar, Dynamic, 1> DynVectorType;
|
|
typedef Matrix<Scalar, 1, Dynamic> DynRowVectorType;
|
|
|
|
Index rows = m.rows();
|
|
Index depth = internal::random<Index>(1, EIGEN_TEST_MAX_SIZE);
|
|
|
|
VectorType v1 = VectorType::Random(rows);
|
|
DynVectorType dv1 = DynVectorType::Random(depth);
|
|
DynRowVectorType drv1 = DynRowVectorType::Random(depth);
|
|
DynMatrixType dm1 = dv1;
|
|
DynMatrixType drm1 = drv1;
|
|
|
|
Scalar s = v1(0);
|
|
|
|
VERIFY_IS_APPROX(v1.asDiagonal() * drv1, s * drv1);
|
|
VERIFY_IS_APPROX(dv1 * v1.asDiagonal(), dv1 * s);
|
|
|
|
VERIFY_IS_APPROX(v1.asDiagonal() * drm1, s * drm1);
|
|
VERIFY_IS_APPROX(dm1 * v1.asDiagonal(), dm1 * s);
|
|
}
|
|
|
|
template <int>
|
|
void bug987() {
|
|
Matrix3Xd points = Matrix3Xd::Random(3, 3);
|
|
Vector2d diag = Vector2d::Random();
|
|
Matrix2Xd tmp1 = points.topRows<2>(), res1, res2;
|
|
VERIFY_IS_APPROX(res1 = diag.asDiagonal() * points.topRows<2>(), res2 = diag.asDiagonal() * tmp1);
|
|
Matrix2d tmp2 = points.topLeftCorner<2, 2>();
|
|
VERIFY_IS_APPROX((res1 = points.topLeftCorner<2, 2>() * diag.asDiagonal()), res2 = tmp2 * diag.asDiagonal());
|
|
}
|
|
|
|
template <int>
|
|
void bug2013() {
|
|
Matrix3d m = Matrix3d::Random();
|
|
Vector3d d = Vector3d::Random();
|
|
|
|
Matrix3d ref_unit_lower = m.template triangularView<UnitLower>();
|
|
Matrix3d ref_unit_upper = m.template triangularView<UnitUpper>();
|
|
Matrix3d ref_lower = m.template triangularView<Lower>();
|
|
Matrix3d ref_upper = m.template triangularView<Upper>();
|
|
|
|
VERIFY_IS_APPROX((m.template triangularView<UnitLower>() * d.asDiagonal()).eval(), ref_unit_lower * d.asDiagonal());
|
|
VERIFY_IS_APPROX((d.asDiagonal() * m.template triangularView<UnitLower>()).eval(), d.asDiagonal() * ref_unit_lower);
|
|
|
|
VERIFY_IS_APPROX((m.template triangularView<UnitUpper>() * d.asDiagonal()).eval(), ref_unit_upper * d.asDiagonal());
|
|
VERIFY_IS_APPROX((d.asDiagonal() * m.template triangularView<UnitUpper>()).eval(), d.asDiagonal() * ref_unit_upper);
|
|
|
|
Matrix3d actual = Matrix3d::Random();
|
|
Matrix3d expected = actual;
|
|
actual = m;
|
|
expected = m;
|
|
actual.template triangularView<Upper>() = actual.template triangularView<Upper>() * d.asDiagonal();
|
|
expected.template triangularView<Upper>() = Matrix3d(ref_upper * d.asDiagonal());
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
actual = m;
|
|
expected = m;
|
|
actual.template triangularView<Lower>() = d.asDiagonal() * actual.template triangularView<Lower>();
|
|
expected.template triangularView<Lower>() = Matrix3d(d.asDiagonal() * ref_lower);
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
actual.setRandom();
|
|
expected = actual;
|
|
actual.noalias() += m.template triangularView<UnitLower>() * d.asDiagonal();
|
|
expected.noalias() += ref_unit_lower * d.asDiagonal();
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
actual.setRandom();
|
|
expected = actual;
|
|
actual.noalias() -= d.asDiagonal() * m.template triangularView<UnitUpper>();
|
|
expected.noalias() -= d.asDiagonal() * ref_unit_upper;
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
MatrixXd dynamic_m = MatrixXd::Random(4, 4);
|
|
VectorXd dynamic_d = VectorXd::Random(4);
|
|
MatrixXd dynamic_expected(4, 4);
|
|
MatrixXd no_malloc_result(4, 4);
|
|
|
|
dynamic_expected = MatrixXd(dynamic_m.template triangularView<UnitLower>()) * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.noalias() = dynamic_m.template triangularView<UnitLower>() * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
|
|
dynamic_expected = dynamic_d.asDiagonal() * MatrixXd(dynamic_m.template triangularView<UnitUpper>());
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.noalias() = dynamic_d.asDiagonal() * dynamic_m.template triangularView<UnitUpper>();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
|
|
// Triangular destination assignment from a structured*diagonal product must go through the
|
|
// lazy product_evaluator and avoid materializing a full PlainObject temporary.
|
|
no_malloc_result = dynamic_m;
|
|
dynamic_expected = dynamic_m;
|
|
dynamic_expected.template triangularView<Upper>() =
|
|
MatrixXd(dynamic_m.template triangularView<Upper>()) * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.template triangularView<Upper>() =
|
|
dynamic_m.template triangularView<Upper>() * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
|
|
no_malloc_result = dynamic_m;
|
|
dynamic_expected = dynamic_m;
|
|
dynamic_expected.template triangularView<Lower>() =
|
|
dynamic_d.asDiagonal() * MatrixXd(dynamic_m.template triangularView<Lower>());
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.template triangularView<Lower>() =
|
|
dynamic_d.asDiagonal() * dynamic_m.template triangularView<Lower>();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
}
|
|
|
|
template <int>
|
|
void selfadjoint_diagonal_products() {
|
|
Matrix3cd m = Matrix3cd::Random();
|
|
m.diagonal() = m.diagonal().real();
|
|
Vector3cd d = Vector3cd::Random();
|
|
|
|
Matrix3cd ref_lower = m.template selfadjointView<Lower>();
|
|
Matrix3cd ref_upper = m.template selfadjointView<Upper>();
|
|
|
|
VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * d.asDiagonal()).eval(), ref_lower * d.asDiagonal());
|
|
VERIFY_IS_APPROX((d.asDiagonal() * m.template selfadjointView<Lower>()).eval(), d.asDiagonal() * ref_lower);
|
|
|
|
VERIFY_IS_APPROX((m.template selfadjointView<Upper>() * d.asDiagonal()).eval(), ref_upper * d.asDiagonal());
|
|
VERIFY_IS_APPROX((d.asDiagonal() * m.template selfadjointView<Upper>()).eval(), d.asDiagonal() * ref_upper);
|
|
|
|
Matrix3cd actual = Matrix3cd::Random();
|
|
Matrix3cd expected = actual;
|
|
actual = m;
|
|
expected = m;
|
|
actual.template selfadjointView<Upper>() = actual.template selfadjointView<Upper>() * d.asDiagonal();
|
|
expected.template triangularView<Upper>() = Matrix3cd(ref_upper * d.asDiagonal());
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
actual = m;
|
|
expected = m;
|
|
actual.template selfadjointView<Lower>() = d.asDiagonal() * actual.template selfadjointView<Lower>();
|
|
expected.template triangularView<Lower>() = Matrix3cd(d.asDiagonal() * ref_lower);
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
actual.setRandom();
|
|
expected = actual;
|
|
actual.noalias() += m.template selfadjointView<Lower>() * d.asDiagonal();
|
|
expected.noalias() += ref_lower * d.asDiagonal();
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
actual.setRandom();
|
|
expected = actual;
|
|
actual.noalias() -= d.asDiagonal() * m.template selfadjointView<Upper>();
|
|
expected.noalias() -= d.asDiagonal() * ref_upper;
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
|
|
MatrixXcd dynamic_m = MatrixXcd::Random(4, 4);
|
|
dynamic_m.diagonal() = dynamic_m.diagonal().real();
|
|
VectorXcd dynamic_d = VectorXcd::Random(4);
|
|
MatrixXcd dynamic_expected(4, 4);
|
|
MatrixXcd no_malloc_result(4, 4);
|
|
|
|
dynamic_expected = MatrixXcd(dynamic_m.template selfadjointView<Lower>()) * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.noalias() = dynamic_m.template selfadjointView<Lower>() * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
|
|
dynamic_expected = dynamic_d.asDiagonal() * MatrixXcd(dynamic_m.template selfadjointView<Upper>());
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.noalias() = dynamic_d.asDiagonal() * dynamic_m.template selfadjointView<Upper>();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
|
|
// Triangular destination assignment from selfadjoint*diagonal must use the lazy
|
|
// product_evaluator (with conjugate-mirror coeffs) and not materialize a temporary.
|
|
no_malloc_result = dynamic_m;
|
|
dynamic_expected = dynamic_m;
|
|
dynamic_expected.template triangularView<Lower>() =
|
|
MatrixXcd(dynamic_m.template selfadjointView<Upper>()) * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.template triangularView<Lower>() =
|
|
dynamic_m.template selfadjointView<Upper>() * dynamic_d.asDiagonal();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
|
|
no_malloc_result = dynamic_m;
|
|
dynamic_expected = dynamic_m;
|
|
dynamic_expected.template triangularView<Upper>() =
|
|
dynamic_d.asDiagonal() * MatrixXcd(dynamic_m.template selfadjointView<Lower>());
|
|
internal::set_is_malloc_allowed(false);
|
|
no_malloc_result.template triangularView<Upper>() =
|
|
dynamic_d.asDiagonal() * dynamic_m.template selfadjointView<Lower>();
|
|
internal::set_is_malloc_allowed(true);
|
|
VERIFY_IS_APPROX(no_malloc_result, dynamic_expected);
|
|
}
|
|
|
|
// Exercise the block-tile path of the dense selfadjoint x diagonal kernel
|
|
// (BlockSize = 32 in ProductEvaluators.h). Picks a few sizes that hit
|
|
// full blocks (size = 64), partial blocks (size = 33, 65), and a tiny size
|
|
// that bypasses the block loop entirely (size = 8). Also verifies the
|
|
// overwrite path leaves no stale data when dst is pre-filled.
|
|
template <typename Scalar>
|
|
void selfadjoint_diagonal_products_at(Index n) {
|
|
typedef Matrix<Scalar, Dynamic, Dynamic> MatType;
|
|
typedef Matrix<Scalar, Dynamic, 1> VecType;
|
|
|
|
MatType m = MatType::Random(n, n);
|
|
m.diagonal() = m.diagonal().real().template cast<Scalar>(); // Hermitian diagonal
|
|
VecType d = VecType::Random(n);
|
|
|
|
MatType ref_lower = m.template selfadjointView<Lower>();
|
|
MatType ref_upper = m.template selfadjointView<Upper>();
|
|
|
|
// Plain assignment goes through evalTo (overwrite kernel).
|
|
// Pre-fill dst with garbage to verify no stale entries remain.
|
|
MatType dst = MatType::Constant(n, n, Scalar(42));
|
|
dst.noalias() = m.template selfadjointView<Upper>() * d.asDiagonal();
|
|
VERIFY_IS_APPROX(dst, ref_upper * d.asDiagonal());
|
|
|
|
dst = MatType::Constant(n, n, Scalar(-7));
|
|
dst.noalias() = m.template selfadjointView<Lower>() * d.asDiagonal();
|
|
VERIFY_IS_APPROX(dst, ref_lower * d.asDiagonal());
|
|
|
|
dst = MatType::Constant(n, n, Scalar(13));
|
|
dst.noalias() = d.asDiagonal() * m.template selfadjointView<Upper>();
|
|
VERIFY_IS_APPROX(dst, d.asDiagonal() * ref_upper);
|
|
|
|
dst = MatType::Constant(n, n, Scalar(99));
|
|
dst.noalias() = d.asDiagonal() * m.template selfadjointView<Lower>();
|
|
VERIFY_IS_APPROX(dst, d.asDiagonal() * ref_lower);
|
|
|
|
// Accumulating paths (scaleAndAddTo).
|
|
MatType base = MatType::Random(n, n);
|
|
dst = base;
|
|
dst.noalias() += m.template selfadjointView<Upper>() * d.asDiagonal();
|
|
VERIFY_IS_APPROX(dst, base + ref_upper * d.asDiagonal());
|
|
|
|
dst = base;
|
|
dst.noalias() -= d.asDiagonal() * m.template selfadjointView<Lower>();
|
|
VERIFY_IS_APPROX(dst, base - d.asDiagonal() * ref_lower);
|
|
|
|
// Scalar-scaled products: the "Dense ?= scalar * Product" rewriting rule
|
|
// folds alpha into the SelfAdjointView. For a complex alpha that fold is
|
|
// not Hermitian, so the dispatch must restore alpha rather than apply it
|
|
// straight to the kernel — verify all four orientation x triangle combos.
|
|
Scalar alpha = internal::random<Scalar>();
|
|
dst = base;
|
|
dst.noalias() += alpha * (m.template selfadjointView<Lower>() * d.asDiagonal());
|
|
VERIFY_IS_APPROX(dst, base + alpha * (ref_lower * d.asDiagonal()));
|
|
|
|
dst = base;
|
|
dst.noalias() -= alpha * (m.template selfadjointView<Upper>() * d.asDiagonal());
|
|
VERIFY_IS_APPROX(dst, base - alpha * (ref_upper * d.asDiagonal()));
|
|
|
|
dst = base;
|
|
dst.noalias() += alpha * (d.asDiagonal() * m.template selfadjointView<Upper>());
|
|
VERIFY_IS_APPROX(dst, base + alpha * (d.asDiagonal() * ref_upper));
|
|
|
|
dst = base;
|
|
dst.noalias() -= alpha * (d.asDiagonal() * m.template selfadjointView<Lower>());
|
|
VERIFY_IS_APPROX(dst, base - alpha * (d.asDiagonal() * ref_lower));
|
|
|
|
// Overwrite-with-scalar path: hits evalTo's HasScalarFactor branch.
|
|
dst = MatType::Constant(n, n, Scalar(17));
|
|
dst.noalias() = alpha * (m.template selfadjointView<Lower>() * d.asDiagonal());
|
|
VERIFY_IS_APPROX(dst, alpha * (ref_lower * d.asDiagonal()));
|
|
|
|
dst = MatType::Constant(n, n, Scalar(-3));
|
|
dst.noalias() = alpha * (d.asDiagonal() * m.template selfadjointView<Upper>());
|
|
VERIFY_IS_APPROX(dst, alpha * (d.asDiagonal() * ref_upper));
|
|
|
|
// Conjugated nested expressions go through the same blas_traits extraction
|
|
// path. The extracted matrix must keep NeedToConjugate, otherwise the kernel
|
|
// computes with m instead of m.conjugate().
|
|
MatType conj_ref_lower = m.conjugate().template selfadjointView<Lower>();
|
|
MatType conj_ref_upper = m.conjugate().template selfadjointView<Upper>();
|
|
|
|
dst = MatType::Constant(n, n, Scalar(23));
|
|
dst.noalias() = m.conjugate().template selfadjointView<Upper>() * d.asDiagonal();
|
|
VERIFY_IS_APPROX(dst, conj_ref_upper * d.asDiagonal());
|
|
|
|
dst = MatType::Constant(n, n, Scalar(-29));
|
|
dst.noalias() = d.asDiagonal() * m.conjugate().template selfadjointView<Lower>();
|
|
VERIFY_IS_APPROX(dst, d.asDiagonal() * conj_ref_lower);
|
|
|
|
dst = base;
|
|
dst.noalias() += alpha * (m.conjugate().template selfadjointView<Lower>() * d.asDiagonal());
|
|
VERIFY_IS_APPROX(dst, base + alpha * (conj_ref_lower * d.asDiagonal()));
|
|
|
|
dst = base;
|
|
dst.noalias() -= alpha * (d.asDiagonal() * m.conjugate().template selfadjointView<Upper>());
|
|
VERIFY_IS_APPROX(dst, base - alpha * (d.asDiagonal() * conj_ref_upper));
|
|
}
|
|
|
|
template <int>
|
|
void selfadjoint_diagonal_products_block_path() {
|
|
selfadjoint_diagonal_products_at<double>(8);
|
|
selfadjoint_diagonal_products_at<double>(33); // partial off-diagonal block
|
|
selfadjoint_diagonal_products_at<double>(64); // exact multiple of BlockSize
|
|
selfadjoint_diagonal_products_at<double>(65); // off-by-one
|
|
selfadjoint_diagonal_products_at<std::complex<double>>(33);
|
|
selfadjoint_diagonal_products_at<std::complex<double>>(65);
|
|
}
|
|
|
|
// In-place patterns where the diagonal operand shares storage with the destination view.
|
|
// The fast path in triangular_product_assignment_dispatcher detects the run-time overlap
|
|
// and materializes a temporary, so the result must match a reference computed against a
|
|
// materialized source. The structured (triangular/selfadjoint) operand can safely share
|
|
// storage with dst because the kernel reads each cell before writing it.
|
|
template <unsigned int Mode, typename Mat>
|
|
void verify_triangular_in_place_with_aliased_diagonal(const Mat& m) {
|
|
// diagonal * tri_view
|
|
{
|
|
Mat actual = m, expected = m;
|
|
Mat ref_diag = actual.diagonal().asDiagonal();
|
|
Mat ref_tri = actual.template triangularView<Mode>();
|
|
expected.template triangularView<Mode>() = (ref_diag * ref_tri).eval();
|
|
actual.template triangularView<Mode>() = actual.diagonal().asDiagonal() * actual.template triangularView<Mode>();
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
}
|
|
// tri_view * diagonal
|
|
{
|
|
Mat actual = m, expected = m;
|
|
Mat ref_diag = actual.diagonal().asDiagonal();
|
|
Mat ref_tri = actual.template triangularView<Mode>();
|
|
expected.template triangularView<Mode>() = (ref_tri * ref_diag).eval();
|
|
actual.template triangularView<Mode>() = actual.template triangularView<Mode>() * actual.diagonal().asDiagonal();
|
|
VERIFY_IS_APPROX(actual, expected);
|
|
}
|
|
}
|
|
|
|
template <unsigned int Mode, typename Mat>
|
|
void verify_selfadjoint_in_place_with_aliased_diagonal(const Mat& m) {
|
|
// diagonal * sa_view
|
|
{
|
|
Mat actual = m, expected = m;
|
|
Mat ref_diag = actual.diagonal().asDiagonal();
|
|
Mat ref_sa = actual.template selfadjointView<Mode>();
|
|
expected.template triangularView<Mode>() = (ref_diag * ref_sa).eval();
|
|
actual.template selfadjointView<Mode>() = actual.diagonal().asDiagonal() * actual.template selfadjointView<Mode>();
|
|
VERIFY_IS_APPROX(actual.template triangularView<Mode>().toDenseMatrix(),
|
|
expected.template triangularView<Mode>().toDenseMatrix());
|
|
}
|
|
// sa_view * diagonal
|
|
{
|
|
Mat actual = m, expected = m;
|
|
Mat ref_diag = actual.diagonal().asDiagonal();
|
|
Mat ref_sa = actual.template selfadjointView<Mode>();
|
|
expected.template triangularView<Mode>() = (ref_sa * ref_diag).eval();
|
|
actual.template selfadjointView<Mode>() = actual.template selfadjointView<Mode>() * actual.diagonal().asDiagonal();
|
|
VERIFY_IS_APPROX(actual.template triangularView<Mode>().toDenseMatrix(),
|
|
expected.template triangularView<Mode>().toDenseMatrix());
|
|
}
|
|
}
|
|
|
|
template <int>
|
|
void structured_diagonal_aliasing() {
|
|
for (int n : {3, 5, 8, 17, 32, 33, 64, 65}) {
|
|
MatrixXcd m = MatrixXcd::Random(n, n);
|
|
m.diagonal() = m.diagonal().real(); // Hermitian-friendly diagonal
|
|
|
|
verify_triangular_in_place_with_aliased_diagonal<Upper>(m);
|
|
verify_triangular_in_place_with_aliased_diagonal<Lower>(m);
|
|
verify_selfadjoint_in_place_with_aliased_diagonal<Upper>(m);
|
|
verify_selfadjoint_in_place_with_aliased_diagonal<Lower>(m);
|
|
}
|
|
}
|
|
|
|
EIGEN_DECLARE_TEST(diagonalmatrices) {
|
|
for (int i = 0; i < g_repeat; i++) {
|
|
CALL_SUBTEST_1(diagonalmatrices(Matrix<float, 1, 1>()));
|
|
CALL_SUBTEST_1(as_scalar_product(Matrix<float, 1, 1>()));
|
|
|
|
CALL_SUBTEST_2(diagonalmatrices(Matrix3f()));
|
|
CALL_SUBTEST_3(diagonalmatrices(Matrix<double, 3, 3, RowMajor>()));
|
|
CALL_SUBTEST_4(diagonalmatrices(Matrix4d()));
|
|
CALL_SUBTEST_5(diagonalmatrices(Matrix<float, 4, 4, RowMajor>()));
|
|
CALL_SUBTEST_6(diagonalmatrices(
|
|
MatrixXcf(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
|
|
CALL_SUBTEST_6(as_scalar_product(MatrixXcf(1, 1)));
|
|
CALL_SUBTEST_7(diagonalmatrices(
|
|
MatrixXi(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
|
|
CALL_SUBTEST_8(diagonalmatrices(Matrix<double, Dynamic, Dynamic, RowMajor>(
|
|
internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
|
|
CALL_SUBTEST_9(diagonalmatrices(
|
|
MatrixXf(internal::random<int>(1, EIGEN_TEST_MAX_SIZE), internal::random<int>(1, EIGEN_TEST_MAX_SIZE))));
|
|
CALL_SUBTEST_9(diagonalmatrices(MatrixXf(1, 1)));
|
|
CALL_SUBTEST_9(as_scalar_product(MatrixXf(1, 1)));
|
|
}
|
|
CALL_SUBTEST_10(bug987<0>());
|
|
CALL_SUBTEST_10(bug2013<0>());
|
|
CALL_SUBTEST_10(selfadjoint_diagonal_products<0>());
|
|
CALL_SUBTEST_10(selfadjoint_diagonal_products_block_path<0>());
|
|
CALL_SUBTEST_10(structured_diagonal_aliasing<0>());
|
|
}
|