test_TinyMatrix.cpp

#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <Kokkos_Core.hpp>

#include <utils/PugsAssert.hpp>
#include <utils/Types.hpp>

#include <algebra/TinyMatrix.hpp>

#include <sstream>

// Instantiate to ensure full coverage is performed
template class TinyMatrix<1, 1, int>;
template class TinyMatrix<2, 2, int>;
template class TinyMatrix<3, 4, int>;
template class TinyMatrix<4, 4, double>;

// clazy:excludeall=non-pod-global-static

TEST_CASE("TinyMatrix", "[algebra]")
{
  REQUIRE(TinyMatrix<1, 1, int>::Dimension == 1);
  REQUIRE(TinyMatrix<1, 1, int>::NumberOfRows == 1);
  REQUIRE(TinyMatrix<1, 1, int>::NumberOfColumns == 1);
  REQUIRE(TinyMatrix<2, 3, int>::Dimension == 6);
  REQUIRE(TinyMatrix<2, 3, int>::NumberOfRows == 2);
  REQUIRE(TinyMatrix<2, 3, int>::NumberOfColumns == 3);
  REQUIRE(TinyMatrix<5, 4, int>::Dimension == 20);
  REQUIRE(TinyMatrix<5, 4, int>::NumberOfRows == 5);
  REQUIRE(TinyMatrix<5, 4, int>::NumberOfColumns == 4);

  TinyMatrix<3, 4, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12);
  TinyVector<3, int> x1(1, 5, 9);
  TinyVector<3, int> x2(2, 6, 10);
  TinyVector<3, int> x3(3, 7, 11);
  TinyVector<3, int> x4(4, 8, 12);
  REQUIRE(((A(0, 0) == 1) and (A(0, 1) == 2) and (A(0, 2) == 3) and (A(0, 3) == 4) and   //
           (A(1, 0) == 5) and (A(1, 1) == 6) and (A(1, 2) == 7) and (A(1, 3) == 8) and   //
           (A(2, 0) == 9) and (A(2, 1) == 10) and (A(2, 2) == 11) and (A(2, 3) == 12)));
  REQUIRE(A == TinyMatrix<3, 4, int>(x1, x2, x3, x4));
  TinyMatrix<3, 4, int> B(6, 5, 3, 8, 34, 6, 35, 6, 7, 1, 3, 6);

  SECTION("checking for opposed matrix")
  {
    const TinyMatrix minus_A = -A;
    REQUIRE(
      ((minus_A(0, 0) == -1) and (minus_A(0, 1) == -2) and (minus_A(0, 2) == -3) and (minus_A(0, 3) == -4) and   //
       (minus_A(1, 0) == -5) and (minus_A(1, 1) == -6) and (minus_A(1, 2) == -7) and (minus_A(1, 3) == -8) and   //
       (minus_A(2, 0) == -9) and (minus_A(2, 1) == -10) and (minus_A(2, 2) == -11) and (minus_A(2, 3) == -12)));
  }

  SECTION("checking for equality and difference tests")
  {
    const TinyMatrix copy_A = A;
    REQUIRE(((copy_A(0, 0) == 1) and (copy_A(0, 1) == 2) and (copy_A(0, 2) == 3) and (copy_A(0, 3) == 4) and   //
             (copy_A(1, 0) == 5) and (copy_A(1, 1) == 6) and (copy_A(1, 2) == 7) and (copy_A(1, 3) == 8) and   //
             (copy_A(2, 0) == 9) and (copy_A(2, 1) == 10) and (copy_A(2, 2) == 11) and (copy_A(2, 3) == 12)));
    REQUIRE(copy_A == A);
    REQUIRE_FALSE(copy_A != A);

    TinyMatrix<3, 4, int> affected_A;
    affected_A = A;
    REQUIRE(affected_A == A);
    REQUIRE_FALSE(affected_A != A);

    REQUIRE(A != B);
    REQUIRE_FALSE(A == B);
  }

  SECTION("checking for scalar left product")
  {
    const int a                    = 2;
    const TinyMatrix<3, 4, int> aA = a * A;

    REQUIRE(aA == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
  }

  SECTION("checking for scalar seft product")
  {
    const int a = 2;

    TinyMatrix copy_A = A;

    REQUIRE((copy_A *= a) == TinyMatrix<3, 4, int>(2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24));
  }

  SECTION("checking for null matrix management")
  {
    TinyMatrix<4, 3, int> Z = zero;
    REQUIRE(Z == TinyMatrix<4, 3, int>(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));

    TinyMatrix<4, 3, int> affected_Z;
    affected_Z = zero;
    REQUIRE(affected_Z == Z);
  }

  SECTION("checking for identity management")
  {
    TinyMatrix<3, 3, int> I = identity;
    REQUIRE(I == TinyMatrix<3, 3, int>(1, 0, 0, 0, 1, 0, 0, 0, 1));

    TinyMatrix<3, 3, int> affected_I;
    affected_I = identity;
    REQUIRE(affected_I == I);
  }

  SECTION("checking for matrices sum")
  {
    REQUIRE(A + B == TinyMatrix<3, 4, int>(7, 7, 6, 12, 39, 12, 42, 14, 16, 11, 14, 18));

    TinyMatrix<3, 4, int> ApB = A;
    ApB += B;
    REQUIRE(ApB == A + B);

    TinyMatrix<3, 4, int> Ap2B = A + 2 * B;
    REQUIRE(Ap2B == ApB + B);
  }

  SECTION("checking for matrices difference ")
  {
    REQUIRE(A - B == TinyMatrix<3, 4, int>(-5, -3, 0, -4, -29, 0, -28, 2, 2, 9, 8, 6));

    TinyMatrix<3, 4, int> AmB = A;
    AmB -= B;
    REQUIRE(AmB == A - B);

    TinyMatrix<3, 4, int> Am2B = A - 2 * B;
    REQUIRE(Am2B == AmB - B);
  }

  SECTION("checking for matrices product")
  {
    TinyMatrix<4, 2, int> C{3, -2, 2, 6, -2, 5, 7, 2};
    REQUIRE(A * C == TinyMatrix<3, 2, int>(29, 33, 69, 77, 109, 121));

    TinyMatrix<2, 3, int> D{-3, 2, 3, 2, 5, -7};
    REQUIRE(D * A == TinyMatrix<2, 4, int>(34, 36, 38, 40, -36, -36, -36, -36));
  }

  SECTION("checking for matrix-vector product")
  {
    REQUIRE(A * TinyVector<4, int>(2, -3, 5, 2) == TinyVector<3, int>(19, 43, 67));
  }

  SECTION("checking for tensor product")
  {
    const TinyVector<4, int> u(1, 3, 7, 5);
    const TinyVector<3, int> v(6, 2, -3);

    REQUIRE(tensorProduct(u, v) == TinyMatrix<4, 3, int>(6, 2, -3, 18, 6, -9, 42, 14, -21, 30, 10, -15));
  }

  SECTION("checking for minor calculation")
  {
    REQUIRE(getMinor(A, 0, 0) == TinyMatrix<2, 3, int>(6, 7, 8, 10, 11, 12));
    REQUIRE(getMinor(A, 1, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 10, 11, 12));
    REQUIRE(getMinor(A, 2, 0) == TinyMatrix<2, 3, int>(2, 3, 4, 6, 7, 8));

    REQUIRE(getMinor(A, 0, 1) == TinyMatrix<2, 3, int>(5, 7, 8, 9, 11, 12));
    REQUIRE(getMinor(A, 1, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 9, 11, 12));
    REQUIRE(getMinor(A, 2, 1) == TinyMatrix<2, 3, int>(1, 3, 4, 5, 7, 8));

    REQUIRE(getMinor(A, 0, 2) == TinyMatrix<2, 3, int>(5, 6, 8, 9, 10, 12));
    REQUIRE(getMinor(A, 1, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 9, 10, 12));
    REQUIRE(getMinor(A, 2, 2) == TinyMatrix<2, 3, int>(1, 2, 4, 5, 6, 8));

    REQUIRE(getMinor(A, 0, 3) == TinyMatrix<2, 3, int>(5, 6, 7, 9, 10, 11));
    REQUIRE(getMinor(A, 1, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 9, 10, 11));
    REQUIRE(getMinor(A, 2, 3) == TinyMatrix<2, 3, int>(1, 2, 3, 5, 6, 7));
  }

  SECTION("checking for cofactors")
  {
    TinyMatrix<3, 3, int> A(1, 2, 3, 4, 5, 6, 7, 8, 9);
    REQUIRE(cofactor(A, 0, 0) == (5 * 9 - 8 * 6));
    REQUIRE(cofactor(A, 1, 0) == -(2 * 9 - 8 * 3));
    REQUIRE(cofactor(A, 2, 0) == (2 * 6 - 5 * 3));

    REQUIRE(cofactor(A, 0, 1) == -(4 * 9 - 7 * 6));
    REQUIRE(cofactor(A, 1, 1) == (1 * 9 - 7 * 3));
    REQUIRE(cofactor(A, 2, 1) == -(1 * 6 - 4 * 3));

    REQUIRE(cofactor(A, 0, 2) == (4 * 8 - 5 * 7));
    REQUIRE(cofactor(A, 1, 2) == -(1 * 8 - 7 * 2));
    REQUIRE(cofactor(A, 2, 2) == (1 * 5 - 4 * 2));
  }

  SECTION("checking for determinant calculations")
  {
    REQUIRE(det(TinyMatrix<1, 1, int>(6)) == 6);
    REQUIRE(det(TinyMatrix<2, 2, int>(3, 1, -3, 6)) == 21);
    REQUIRE(det(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1);
    REQUIRE(det(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == -1444);
    REQUIRE(det(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
            Catch::Approx(6661.455).epsilon(1E-14));
    REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
  }

  SECTION("checking for trace calculations")
  {
    REQUIRE(trace(TinyMatrix<1, 1, int>(6)) == 6);
    REQUIRE(trace(TinyMatrix<2, 2, int>(5, 1, -3, 6)) == 5 + 6);
    REQUIRE(trace(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1 + 2 + 3);
    REQUIRE(trace(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == 6 + 34 + 7);
    REQUIRE(trace(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
            Catch::Approx(1 + 4 + 2 + 17.5).epsilon(1E-14));
    REQUIRE(trace(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 1 + 0 + 1 + 2);
  }

  SECTION("checking for inverse calculations")
  {
    {
      const TinyMatrix<1> A1(2);
      REQUIRE(inverse(A1)(0, 0) == Catch::Approx(0.5).epsilon(1E-14));
    }

    {
      const TinyMatrix<2> A2(2, 3, 4, 1);
      const TinyMatrix<2> I = inverse(A2) * A2;

      REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
      REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(1, 0) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(1, 1) == Catch::Approx(1).epsilon(1E-14));
    }

    {
      const TinyMatrix<3> A3(2, 3, 1, 4, -1, 5, -2, 3, 4);
      const TinyMatrix<3> I = inverse(A3) * A3;

      REQUIRE(I(0, 0) == Catch::Approx(1).epsilon(1E-14));
      REQUIRE(I(0, 1) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(0, 2) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(1, 0) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(1, 1) == Catch::Approx(1).epsilon(1E-14));
      REQUIRE(I(1, 2) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(2, 0) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(2, 1) == Catch::Approx(0).margin(1E-14));
      REQUIRE(I(2, 2) == Catch::Approx(1).epsilon(1E-14));
    }
  }

  SECTION("checking for sizes")
  {
    REQUIRE(TinyMatrix<1>{}.numberOfRows() == 1);
    REQUIRE(TinyMatrix<1>{}.numberOfColumns() == 1);
    REQUIRE(TinyMatrix<1>{}.dimension() == 1);
    REQUIRE(TinyMatrix<1>{}.numberOfValues() == 1);

    REQUIRE(TinyMatrix<2>{}.numberOfRows() == 2);
    REQUIRE(TinyMatrix<2>{}.numberOfColumns() == 2);
    REQUIRE(TinyMatrix<2>{}.numberOfValues() == 4);

    REQUIRE(TinyMatrix<3>{}.numberOfRows() == 3);
    REQUIRE(TinyMatrix<3>{}.numberOfColumns() == 3);
    REQUIRE(TinyMatrix<3>{}.dimension() == 9);
    REQUIRE(TinyMatrix<3>{}.numberOfValues() == 9);

    REQUIRE(TinyMatrix<3, 4>{}.numberOfRows() == 3);
    REQUIRE(TinyMatrix<3, 4>{}.numberOfColumns() == 4);
    REQUIRE(TinyMatrix<3, 4>{}.dimension() == 12);
    REQUIRE(TinyMatrix<3, 4>{}.numberOfValues() == 12);
  }

  SECTION("is square matrix")
  {
    REQUIRE(TinyMatrix<1>{}.isSquare());
    REQUIRE(TinyMatrix<2>{}.isSquare());
    REQUIRE(TinyMatrix<3>{}.isSquare());

    REQUIRE_FALSE(TinyMatrix<3, 4>{}.isSquare());
  }

  SECTION("transpose")
  {
    TinyMatrix tA = transpose(A);

    REQUIRE(((tA(0, 0) == 1) and (tA(1, 0) == 2) and (tA(2, 0) == 3) and (tA(3, 0) == 4) and   //
             (tA(0, 1) == 5) and (tA(1, 1) == 6) and (tA(2, 1) == 7) and (tA(3, 1) == 8) and   //
             (tA(0, 2) == 9) and (tA(1, 2) == 10) and (tA(2, 2) == 11) and (tA(3, 2) == 12)));

    TinyMatrix ttA = transpose(tA);
    REQUIRE(ttA == A);
  }

  SECTION("checking for matrices output")
  {
    REQUIRE(Catch::Detail::stringify(A) == "[[1,2,3,4],[5,6,7,8],[9,10,11,12]]");
    REQUIRE(Catch::Detail::stringify(TinyMatrix<1, 1, int>(7)) == "[[7]]");
  }

#ifndef NDEBUG
  SECTION("output with signaling NaN")
  {
    TinyMatrix<2, 3> A;
    A(0, 0) = 1;
    A(0, 2) = 3;
    A(1, 0) = 2;
    std::ostringstream A_ost;
    A_ost << A;

    std::ostringstream ref_ost;
    ref_ost << "[[1,nan,3],[2,nan,nan]]";

    REQUIRE(A_ost.str() == ref_ost.str());
  }

  SECTION("checking for bounds violation")
  {
    REQUIRE_THROWS_AS(A(3, 0), AssertError);
    REQUIRE_THROWS_AS(A(0, 4), AssertError);

    REQUIRE_THROWS_AS(getMinor(A, 3, 0), AssertError);
    REQUIRE_THROWS_AS(getMinor(A, 0, 4), AssertError);

    const TinyMatrix<3, 4, int>& constA = A;
    REQUIRE_THROWS_AS(constA(3, 0), AssertError);
    REQUIRE_THROWS_AS(constA(0, 4), AssertError);
  }

  SECTION("checking for nan initialization")
  {
    TinyMatrix<3, 4, double> B;

    for (size_t i = 0; i < B.numberOfRows(); ++i) {
      for (size_t j = 0; j < B.numberOfColumns(); ++j) {
        REQUIRE(std::isnan(B(i, j)));
      }
    }
  }

  SECTION("checking for bad initialization")
  {
    TinyMatrix<3, 4, int> B;

    for (size_t i = 0; i < B.numberOfRows(); ++i) {
      for (size_t j = 0; j < B.numberOfColumns(); ++j) {
        REQUIRE(B(i, j) == std::numeric_limits<int>::max() / 2);
      }
    }
  }
#endif   // NDEBUG
}