test_EigenvalueSolver.cpp

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

#include <algebra/CRSMatrix.hpp>
#include <algebra/CRSMatrixDescriptor.hpp>
#include <algebra/EigenvalueSolver.hpp>
#include <algebra/SmallMatrix.hpp>

#include <utils/pugs_config.hpp>

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

TEST_CASE("EigenvalueSolver", "[algebra]")
{
  SECTION("Sparse Matrices")
  {
    SECTION("symmetric system")
    {
      Array<int> non_zeros(3);
      non_zeros.fill(3);
      CRSMatrixDescriptor<double> S{3, 3, non_zeros};

      S(0, 0) = 3;
      S(0, 1) = 2;
      S(0, 2) = 4;

      S(1, 0) = 2;
      S(1, 1) = 0;
      S(1, 2) = 2;

      S(2, 0) = 4;
      S(2, 1) = 2;
      S(2, 2) = 3;

      CRSMatrix A{S.getCRSMatrix()};

      SECTION("eigenvalues")
      {
        SmallArray<double> eigenvalues;

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues);
        REQUIRE(eigenvalues[0] == Catch::Approx(-1));
        REQUIRE(eigenvalues[1] == Catch::Approx(-1));
        REQUIRE(eigenvalues[2] == Catch::Approx(8));
#else    //  PUGS_HAS_SLEPC
        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }

      SECTION("eigenvalues and eigenvectors")
      {
        SmallArray<double> eigenvalue_list;
        std::vector<SmallVector<double>> eigenvector_list;

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);

        REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
        REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
        REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

        SmallMatrix<double> P{3};
        SmallMatrix<double> PT{3};
        SmallMatrix<double> D{3};
        D = zero;

        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = 0; j < 3; ++j) {
            P(i, j)  = eigenvector_list[j][i];
            PT(i, j) = eigenvector_list[i][j];
          }
          D(i, i) = eigenvalue_list[i];
        }

        SmallMatrix PDPT = P * D * PT;
        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = 0; j < 3; ++j) {
            REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
          }
        }
#else    //  PUGS_HAS_SLEPC
        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }

      SECTION("eigenvalues and passage matrix")
      {
        SmallArray<double> eigenvalue_list;
        SmallMatrix<double> P{};

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);

        REQUIRE(eigenvalue_list[0] == Catch::Approx(-1));
        REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
        REQUIRE(eigenvalue_list[2] == Catch::Approx(8));

        SmallMatrix<double> D{3};
        D = zero;
        for (size_t i = 0; i < 3; ++i) {
          D(i, i) = eigenvalue_list[i];
        }
        SmallMatrix PT = transpose(P);

        SmallMatrix PDPT = P * D * PT;
        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = 0; j < 3; ++j) {
            REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
          }
        }
#else    //  PUGS_HAS_SLEPC
        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }
      SECTION("exponential of a matrix")
      {
        SmallMatrix<double> expA{};
        SmallMatrix<double> expS{3};

        expS(0, 0) = 1325.074593930307812;
        expS(0, 1) = 662.353357244568285;
        expS(0, 2) = 1324.70671448913637;
        expS(1, 0) = expS(0, 1);
        expS(1, 1) = 331.544558063455535;
        expS(1, 2) = 662.353357244568185;
        expS(2, 0) = expS(0, 2);
        expS(2, 1) = expS(1, 2);
        expS(2, 2) = expS(0, 0);

#ifdef PUGS_HAS_SLEPC
        EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expA);

        for (size_t i = 0; i < 3; ++i) {
          for (size_t j = 0; j < 3; ++j) {
            REQUIRE(expA(i, j) - expS(i, j) == Catch::Approx(0).margin(1E-12));
          }
        }
#else    //  PUGS_HAS_SLEPC
        REQUIRE_THROWS_WITH(EigenvalueSolver{}.computeExpForSymmetricMatrix(A, expP),
                            "not implemented yet: SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
      }
    }
  }
}