eigen/test/triangular.cpp

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

#if defined(EIGEN_TEST_PART_100) || defined(EIGEN_TEST_PART_ALL)
#define EIGEN_NO_DEPRECATED_WARNING
#endif

#include "main.h"

template <typename MatrixType>
void triangular_deprecated(const MatrixType& m) {
  Index rows = m.rows();
  Index cols = m.cols();
  MatrixType m1, m2, m3, m4;
  m1.setRandom(rows, cols);
  m2.setRandom(rows, cols);
  m3 = m1;
  m4 = m2;
  // deprecated method:
  m1.template triangularView<Eigen::Upper>().swap(m2);
  // use this method instead:
  m3.template triangularView<Eigen::Upper>().swap(m4.template triangularView<Eigen::Upper>());
  VERIFY_IS_APPROX(m1, m3);
  VERIFY_IS_APPROX(m2, m4);
  // deprecated method:
  m1.template triangularView<Eigen::Lower>().swap(m4);
  // use this method instead:
  m3.template triangularView<Eigen::Lower>().swap(m2.template triangularView<Eigen::Lower>());
  VERIFY_IS_APPROX(m1, m3);
  VERIFY_IS_APPROX(m2, m4);
}

template <typename MatrixType>
void triangular_square(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>();

  Index rows = m.rows();
  Index cols = m.cols();

  MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), m3(rows, cols), m4(rows, cols),
             r1(rows, cols), r2(rows, cols);
  VectorType v2 = VectorType::Random(rows);
  VectorType v3 = VectorType::Zero(rows);

  MatrixType m1up = m1.template triangularView<Upper>();
  MatrixType m2up = m2.template triangularView<Upper>();

  if (rows * cols > 1) {
    VERIFY(m1up.isUpperTriangular());
    VERIFY(m2up.transpose().isLowerTriangular());
    VERIFY(!m2.isLowerTriangular());
  }

  //   VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2);

  // test overloaded operator+=
  r1.setZero();
  r2.setZero();
  r1.template triangularView<Upper>() += m1;
  r2 += m1up;
  VERIFY_IS_APPROX(r1, r2);

  // test overloaded operator=
  m1.setZero();
  m1.template triangularView<Upper>() = m2.transpose() + m2;
  m3 = m2.transpose() + m2;
  VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1);

  // test overloaded operator=
  m1.setZero();
  m1.template triangularView<Lower>() = m2.transpose() + m2;
  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);

  VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(),
                   m3.conjugate().template triangularView<Lower>().toDenseMatrix());

  m1 = MatrixType::Random(rows, cols);
  for (int i = 0; i < rows; ++i)
    if (numext::abs2(m1(i, i)) < RealScalar(1e-1)) m1(i, i) = Scalar(1);

  Transpose<MatrixType> trm4(m4);
  // test back and forward substitution with a vector as the rhs
  m3 = m1.template triangularView<Upper>();
  v3 = m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2));
  VERIFY(v2.isApprox(v3, largerEps));
  m3 = m1.template triangularView<Lower>();
  v3 = m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2));
  VERIFY(v2.isApprox(v3, largerEps));
  m3 = m1.template triangularView<Upper>();
  v3 = m3 * (m1.template triangularView<Upper>().solve(v2));
  VERIFY(v2.isApprox(v3, largerEps));
  m3 = m1.template triangularView<Lower>();
  v3 = m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2));
  VERIFY(v2.isApprox(v3, largerEps));

  // test back and forward substitution with a matrix as the rhs
  m3 = m1.template triangularView<Upper>();
  m4 = m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2));
  VERIFY(m2.isApprox(m4, largerEps));
  m3 = m1.template triangularView<Lower>();
  m4 = m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2));
  VERIFY(m2.isApprox(m4, largerEps));
  m3 = m1.template triangularView<Upper>();
  m4 = m3 * (m1.template triangularView<Upper>().solve(m2));
  VERIFY(m2.isApprox(m4, largerEps));
  m3 = m1.template triangularView<Lower>();
  m4 = m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2));
  VERIFY(m2.isApprox(m4, largerEps));

  // check M * inv(L) using in place API
  m4 = m3;
  m1.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4);
  VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Lower>(), m3);

  // check M * inv(U) using in place API
  m3 = m1.template triangularView<Upper>();
  m4 = m3;
  m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4);
  VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Upper>(), m3);

  // check solve with unit diagonal
  m3 = m1.template triangularView<UnitUpper>();
  VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps));

  //   VERIFY((  m1.template triangularView<Upper>()
  //           * m2.template triangularView<Upper>()).isUpperTriangular());

  // test swap
  m1.setOnes();
  m2.setZero();
  m2.template triangularView<Upper>().swap(m1.template triangularView<Eigen::Upper>());
  m3.setZero();
  m3.template triangularView<Upper>().setOnes();
  VERIFY_IS_APPROX(m2, m3);

  m1.setRandom();
  m3 = m1.template triangularView<Upper>();
  Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> m5(cols, internal::random<int>(1, 20));
  m5.setRandom();
  Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> m6(internal::random<int>(1, 20), rows);
  m6.setRandom();
  VERIFY_IS_APPROX(m1.template triangularView<Upper>() * m5, m3 * m5);
  VERIFY_IS_APPROX(m6 * m1.template triangularView<Upper>(), m6 * m3);

  m1up = m1.template triangularView<Upper>();
  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
  VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(),
                   m1up.adjoint());
  VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(),
                   m1up.adjoint());

  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal());

  m3.setRandom();
  const MatrixType& m3c(m3);
  VERIFY(is_same_type(m3c.template triangularView<Lower>(),
                      m3.template triangularView<Lower>().template conjugateIf<false>()));
  VERIFY(is_same_type(m3c.template triangularView<Lower>().conjugate(),
                      m3.template triangularView<Lower>().template conjugateIf<true>()));
  VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<true>().toDenseMatrix(),
                   m3.conjugate().template triangularView<Lower>().toDenseMatrix());
  VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<false>().toDenseMatrix(),
                   m3.template triangularView<Lower>().toDenseMatrix());

  VERIFY(is_same_type(m3c.template selfadjointView<Lower>(),
                      m3.template selfadjointView<Lower>().template conjugateIf<false>()));
  VERIFY(is_same_type(m3c.template selfadjointView<Lower>().conjugate(),
                      m3.template selfadjointView<Lower>().template conjugateIf<true>()));
  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<true>().toDenseMatrix(),
                   m3.conjugate().template selfadjointView<Lower>().toDenseMatrix());
  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<false>().toDenseMatrix(),
                   m3.template selfadjointView<Lower>().toDenseMatrix());
}

template <typename MatrixType>
void triangular_rect(const MatrixType& m) {
  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };

  Index rows = m.rows();
  Index cols = m.cols();

  MatrixType m1 = MatrixType::Random(rows, cols), m2 = MatrixType::Random(rows, cols), m3(rows, cols), m4(rows, cols),
             r1(rows, cols), r2(rows, cols);

  MatrixType m1up = m1.template triangularView<Upper>();
  MatrixType m2up = m2.template triangularView<Upper>();

  if (rows > 1 && cols > 1) {
    VERIFY(m1up.isUpperTriangular());
    VERIFY(m2up.transpose().isLowerTriangular());
    VERIFY(!m2.isLowerTriangular());
  }

  // test overloaded operator+=
  r1.setZero();
  r2.setZero();
  r1.template triangularView<Upper>() += m1;
  r2 += m1up;
  VERIFY_IS_APPROX(r1, r2);

  // test overloaded operator=
  m1.setZero();
  m1.template triangularView<Upper>() = 3 * m2;
  m3 = 3 * m2;
  VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);

  m1.setZero();
  m1.template triangularView<Lower>() = 3 * m2;
  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);

  m1.setZero();
  m1.template triangularView<StrictlyUpper>() = 3 * m2;
  VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);

  m1.setZero();
  m1.template triangularView<StrictlyLower>() = 3 * m2;
  VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
  m1.setRandom();
  m2 = m1.template triangularView<Upper>();
  VERIFY(m2.isUpperTriangular());
  VERIFY(!m2.isLowerTriangular());
  m2 = m1.template triangularView<StrictlyUpper>();
  VERIFY(m2.isUpperTriangular());
  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
  m2 = m1.template triangularView<UnitUpper>();
  VERIFY(m2.isUpperTriangular());
  m2.diagonal().array() -= Scalar(1);
  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
  m2 = m1.template triangularView<Lower>();
  VERIFY(m2.isLowerTriangular());
  VERIFY(!m2.isUpperTriangular());
  m2 = m1.template triangularView<StrictlyLower>();
  VERIFY(m2.isLowerTriangular());
  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
  m2 = m1.template triangularView<UnitLower>();
  VERIFY(m2.isLowerTriangular());
  m2.diagonal().array() -= Scalar(1);
  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
  // test swap
  m1.setOnes();
  m2.setZero();
  m2.template triangularView<Upper>().swap(m1.template triangularView<Eigen::Upper>());
  m3.setZero();
  m3.template triangularView<Upper>().setOnes();
  VERIFY_IS_APPROX(m2, m3);
}

// Test triangular solve and product at sizes that exercise GEBP blocking.
// The standard test caps at maxsize=20, which never triggers the blocked code paths
// in TriangularSolverMatrix.h (requires size >= 48 with EIGEN_DEBUG_SMALL_PRODUCT_BLOCKS).
template <int>
void triangular_at_blocking_boundaries() {
  typedef double Scalar;
  typedef Matrix<Scalar, Dynamic, Dynamic> Mat;
  typedef Matrix<Scalar, Dynamic, 1> Vec;

  const int sizes[] = {47, 48, 49, 64, 96, 128};
  for (int si = 0; si < 6; ++si) {
    int n = sizes[si];
    Mat m1 = Mat::Random(n, n);
    // Make well-conditioned: dominant diagonal
    for (int i = 0; i < n; ++i) m1(i, i) += Scalar(n);

    Vec v = Vec::Random(n);
    Mat rhs = Mat::Random(n, 5);

    // Upper triangular solve with vector
    Mat U = m1.triangularView<Upper>();
    Vec x = m1.triangularView<Upper>().solve(v);
    VERIFY_IS_APPROX(U * x, v);

    // Lower triangular solve with vector
    Mat L = m1.triangularView<Lower>();
    x = m1.triangularView<Lower>().solve(v);
    VERIFY_IS_APPROX(L * x, v);

    // Upper triangular solve with matrix rhs
    Mat X = m1.triangularView<Upper>().solve(rhs);
    VERIFY_IS_APPROX(U * X, rhs);

    // Lower triangular solve with matrix rhs
    X = m1.triangularView<Lower>().solve(rhs);
    VERIFY_IS_APPROX(L * X, rhs);

    // Triangular product
    Mat prod = m1.triangularView<Upper>() * rhs;
    VERIFY_IS_APPROX(prod, U * rhs);
    prod = rhs.transpose() * m1.triangularView<Upper>();
    VERIFY_IS_APPROX(prod, rhs.transpose() * U);
  }

  // Also test with float and RowMajor
  {
    typedef Matrix<float, Dynamic, Dynamic, RowMajor> RMat;
    typedef Matrix<float, Dynamic, 1> FVec;
    for (int si = 0; si < 6; ++si) {
      int n = sizes[si];
      RMat m1 = RMat::Random(n, n);
      for (int i = 0; i < n; ++i) m1(i, i) += float(n);

      FVec v = FVec::Random(n);
      RMat U = m1.triangularView<Upper>();
      FVec x = m1.triangularView<Upper>().solve(v);
      VERIFY_IS_APPROX(U * x, v);

      RMat L = m1.triangularView<Lower>();
      x = m1.triangularView<Lower>().solve(v);
      VERIFY_IS_APPROX(L * x, v);
    }
  }
}

void bug_159() {
  Matrix3d m = Matrix3d::Random().triangularView<Lower>();
  EIGEN_UNUSED_VARIABLE(m);
}

EIGEN_DECLARE_TEST(triangular) {
  int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE, 20);
  for (int i = 0; i < g_repeat; i++) {
    int r = internal::random<int>(2, maxsize);
    TEST_SET_BUT_UNUSED_VARIABLE(r);
    int c = internal::random<int>(2, maxsize);
    TEST_SET_BUT_UNUSED_VARIABLE(c);

    CALL_SUBTEST_1(triangular_square(Matrix<float, 1, 1>()));
    CALL_SUBTEST_2(triangular_square(Matrix<float, 2, 2>()));
    CALL_SUBTEST_3(triangular_square(Matrix3d()));
    CALL_SUBTEST_4(triangular_square(Matrix<std::complex<float>, 8, 8>()));
    CALL_SUBTEST_5(triangular_square(MatrixXcd(r, r)));
    CALL_SUBTEST_6(triangular_square(Matrix<float, Dynamic, Dynamic, RowMajor>(r, r)));

    CALL_SUBTEST_7(triangular_rect(Matrix<float, 4, 5>()));
    CALL_SUBTEST_8(triangular_rect(Matrix<double, 6, 2>()));
    CALL_SUBTEST_9(triangular_rect(MatrixXcf(r, c)));
    CALL_SUBTEST_5(triangular_rect(MatrixXcd(r, c)));
    CALL_SUBTEST_6(triangular_rect(Matrix<float, Dynamic, Dynamic, RowMajor>(r, c)));

    CALL_SUBTEST_100(triangular_deprecated(Matrix<float, 5, 7>()));
    CALL_SUBTEST_100(triangular_deprecated(MatrixXd(r, c)));
  }

  CALL_SUBTEST_1(bug_159());

  // Triangular solve/product at blocking boundaries (deterministic, outside g_repeat).
  CALL_SUBTEST_11(triangular_at_blocking_boundaries<0>());
}