eigen_base.hpp

Go to the documentation of this file.

// OpenBEM - Copyright (C) 2026 Shashwat Sharma
 
// This file is part of OpenBEM.
 
// OpenBEM is free software: 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 3
// of the License, or (at your option) any later version.
 
// You should have received a copy of the GNU General Public License along with OpenBEM.
// If not, see <https://www.gnu.org/licenses/>.
 
 
#ifndef EIGEN_MATRIX_BASE_HPP
#define EIGEN_MATRIX_BASE_HPP
 
#include <memory>
 
#include <external/Eigen/Dense>
#include <external/Eigen/IterativeLinearSolvers>
#include <external/EigenUnsupported/Eigen/IterativeSolvers>
 
#include "types.hpp"
#include "matrix/base.hpp"
 
 
namespace bem
{
 
template <typename T, typename MatrixType>
class EigenMatrixBase: public MatrixBase<T>
{
public:
 
    EigenMatrixBase(Index rows, Index cols): matrix_(std::make_shared<MatrixType> ())
    {
        resize(rows, cols);
        set_zero();
        return;
    };
 
 
    EigenMatrixBase(): EigenMatrixBase(0, 0) {};
 
 
    virtual void set_block(
        const MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        const T& a = T(1)
        ) = 0;
 
 
    virtual void add_block(
        const MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        const T& a = T(1)
        ) = 0;
 
 
    virtual void get_block(
        MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        Index b_rows,
        Index b_cols
        ) const = 0;
 
 
    MatrixType& raw_matrix()
    {
        return *matrix_;
    };
 
 
    const MatrixType& raw_matrix() const
    {
        return *matrix_;
    };
 
 
    void set_raw_matrix(const MatrixType& matrix)
    {
        matrix_ = std::make_shared<MatrixType> (matrix);
        return;
    };
 
 
    void bind_raw_matrix(std::shared_ptr<MatrixType> matrix)
    {
        matrix_ = std::move(matrix);
        return;
    };
 
 
    Index num_rows() const override
    {
        return matrix_->rows();
    };
 
 
    Index num_cols() const override
    {
        return matrix_->cols();
    };
 
 
    Index size() const override
    {
        return matrix_->size();
    };
 
 
    void resize(Index rows, Index cols) override
    {
        matrix_->resize(rows, cols);
        set_zero();
        return;
    };
 
 
    void clear() override
    {
        resize(0, 0);
        return;
    };
 
 
    void set_zero() override
    {
        matrix_->setZero();
        return;
    };
 
 
    void set_identity() override
    {
        matrix_->setIdentity();
        return;
    };
 
 
 
    void set_transpose(const MatrixBase<T>& x)
    {
        const EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (x);
        (*matrix_) = xd.raw_matrix().transpose();
        return;
    };
 
 
    void get_diagonal(MatrixBase<T>& x) const
    {
        EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<EigenMatrixBase<T, MatrixType>&> (x);
        xd.raw_matrix() = matrix_->diagonal().asDiagonal();
        return;
    };
 
 
    void add_ax(const MatrixBase<T>& x, const T& a = T(1)) override
    {
        const EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (x);
        raw_matrix() += a * xd.raw_matrix();
        return;
    };
 
 
    void set_axpby(const MatrixBase<T>& x, const MatrixBase<T>& y, const T& a = T(1), const T& b = T(1))
    {
        const EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (x);
        const EigenMatrixBase<T, MatrixType>& yd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (y);
        raw_matrix() = a * xd.raw_matrix() + b * yd.raw_matrix();
        return;
    };
 
 
    void set_mat_mul(const MatrixBase<T>& x, const MatrixBase<T>& y, const T& a = T(1))
    {
        const EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (x);
        const EigenMatrixBase<T, MatrixType>& yd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (y);
        raw_matrix() = xd.raw_matrix() * yd.raw_matrix() * a;
        return;
    };
 
 
    void add_mat_mul(const MatrixBase<T>& x, const MatrixBase<T>& y, const T& a = T(1))
    {
        const EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (x);
        const EigenMatrixBase<T, MatrixType>& yd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (y);
        raw_matrix() += xd.raw_matrix() * yd.raw_matrix() * a;
        return;
    };
 
 
 
    void mat_solve_iterative(MatrixBase<T>& x, const MatrixBase<T>& b, const Float tol = 1e-3)
    {
        mat_solve_gmres(x, b, tol);
        // mat_solve_bicgstab(x, b, tol);
        return;
    };
 
 
    void mat_solve_gmres(
        MatrixBase<T>& x, const MatrixBase<T>& b, const Float tol = 1e-3, const Index restart = 100
        )
    {
        Eigen::GMRES<MatrixType, Eigen::IdentityPreconditioner> solver;
        // Eigen::GMRES<MatrixType, Eigen::DiagonalPreconditioner<T>> solver;
 
        solver.setTolerance(tol);
        solver.setMaxIterations(1000);
        solver.set_restart(restart);
 
        solver.compute(raw_matrix());
        if (solver.info() != Eigen::Success)
            throw std::runtime_error("Matrix solver initialization failed.");
 
        EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<EigenMatrixBase<T, MatrixType>&> (x);
        const EigenMatrixBase<T, MatrixType>& bd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (b);
 
        xd.raw_matrix() = solver.solve(bd.raw_matrix());
 
        std::cout << "GMRES iterations: " << solver.iterations() << " | tolerance: " << solver.tolerance() << " | residual: " << solver.error() << std::endl;
 
        // xd.resize(num_cols(), bd.num_cols());
        // for (Index ii = 0; ii < bd.num_cols(); ++ii)
        //     xd.raw_matrix().col(ii) = solver.solve(bd.raw_matrix().col(ii));
 
        if (solver.info() != Eigen::Success)
            std::cout << "Warning: GMRES solve failed." << std::endl;
 
        return;
    };
 
 
    void mat_solve_bicgstab(MatrixBase<T>& x, const MatrixBase<T>& b, const Float tol = 1e-3)
    {
        Eigen::BiCGSTAB<MatrixType, Eigen::IdentityPreconditioner> solver;
        // Eigen::BiCGSTAB<MatrixType, Eigen::DiagonalPreconditioner<T>> solver;
 
        solver.setTolerance(tol);
        solver.setMaxIterations(1000);
 
        solver.compute(raw_matrix());
        if (solver.info() != Eigen::Success)
            throw std::runtime_error("Matrix solver initialization failed.");
 
        EigenMatrixBase<T, MatrixType>& xd = dynamic_cast<EigenMatrixBase<T, MatrixType>&> (x);
        const EigenMatrixBase<T, MatrixType>& bd = dynamic_cast<const EigenMatrixBase<T, MatrixType>&> (b);
 
        xd.raw_matrix() = solver.solve(bd.raw_matrix());
 
        std::cout << "BiCGStab iterations: " << solver.iterations() << " | tolerance: " << solver.tolerance() << " | residual: " << solver.error() << std::endl;
 
        // xd.resize(num_cols(), bd.num_cols());
        // for (Index ii = 0; ii < bd.num_cols(); ++ii)
        //     xd.raw_matrix().col(ii) = solver.solve(bd.raw_matrix().col(ii));
 
        if (solver.info() != Eigen::Success)
            std::cout << "Warning: BiCGStab solve failed." << std::endl;
 
        return;
    };
 
 
    Float cond() const
    {
        // Eigen::JacobiSVD<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> svd (
        //     (Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>) raw_matrix()
        //     );
        Eigen::BDCSVD<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> svd (
            (Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>) raw_matrix()
            );
        EigColVec<Float> s = svd.singularValues();
        return s(0) / s(s.size() - 1);
    };
 
 
    void print(const std::string name = "matrix") const override
    {
        std::cout << "\n--- " << name << " (" << num_rows() << " x " << num_cols() << ") ---" << std::endl;
        std::cout << raw_matrix() << std::endl;
        std::cout << "---" << std::endl;
        return;
    };
 
 
protected:
 
    std::shared_ptr<MatrixType> matrix_ = nullptr;
 
};
 
}
 
#endif