eigen_dense.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_DENSE_MATRIX_HPP
#define EIGEN_DENSE_MATRIX_HPP
 
#include <external/Eigen/Dense>
#include <external/Eigen/IterativeLinearSolvers>
#include <external/EigenUnsupported/Eigen/IterativeSolvers>
 
#include "types.hpp"
#include "matrix/base.hpp"
#include "matrix/eigen_base.hpp"
 
 
namespace bem
{
 
template <typename T, int StorageOrder = Eigen::ColMajor>
class EigenDenseMatrix: public EigenMatrixBase<T, Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, StorageOrder>>
{
 
    using MatrixType = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, StorageOrder>;
    using base = EigenMatrixBase<T, MatrixType>;
    using base::base;
    using base::matrix_;
 
public:
 
 
    T value(Index row, Index col) const override
    {
        return (*matrix_)(row, col);
    };
 
 
    void set_value(Index row, Index col, const T& a) override
    {
        (*matrix_)(row, col) = a;
        return;
    };
 
 
    void add_value(Index row, Index col, const T& a) override
    {
        (*matrix_)(row, col) += a;
        return;
    };
 
 
    void scale(const T& a) override
    {
        (*matrix_) *= a;
        return;
    };
 
 
    void set_constant(const T& a)
    {
        matrix_->setConstant(a);
        return;
    };
 
 
    void mat_solve(MatrixBase<T>& x, const MatrixBase<T>& b) override
    {
        const EigenDenseMatrix<T>& bd = dynamic_cast<const EigenDenseMatrix<T>&> (b);
        EigenDenseMatrix<T>& xd = dynamic_cast<EigenDenseMatrix<T>&> (x);
        // xd.raw_matrix() = matrix_->householderQr().solve(bd.raw_matrix());
        // xd.raw_matrix() = matrix_->colPivHouseholderQr().solve(bd.raw_matrix());
        // xd.raw_matrix() = matrix_->fullPivHouseholderQr().solve(bd.raw_matrix());
        xd.raw_matrix() = matrix_->partialPivLu().solve(bd.raw_matrix());
        // xd.raw_matrix() = matrix_->fullPivLu().solve(bd.raw_matrix());
        return;
    };
 
 
    void set_block(
        const MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        const T& a = T(1)
        ) override
    {
        const EigenDenseMatrix<T>& xd = dynamic_cast<const EigenDenseMatrix<T>&> (x);
        matrix_->block(row_start, col_start, xd.num_rows(), xd.num_cols()) = xd.raw_matrix() * a;
        return;
    };
 
 
    void add_block(
        const MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        const T& a = T(1)
        ) override
    {
        const EigenDenseMatrix<T>& xd = dynamic_cast<const EigenDenseMatrix<T>&> (x);
        matrix_->block(row_start, col_start, xd.num_rows(), xd.num_cols()) += a * xd.raw_matrix();
        return;
    };
 
 
    void get_block(
        MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        Index b_rows,
        Index b_cols
        ) const override
    {
        EigenDenseMatrix<T>& xd = dynamic_cast<EigenDenseMatrix<T>&> (x);
        xd.resize(b_rows, b_cols);
        xd.raw_matrix() = matrix_->block(row_start, col_start, b_rows, b_cols);
        return;
    };
 
 
    virtual const T* data() const override
    { return matrix_->data(); };
 
 
    virtual T* data() override
    { return matrix_->data(); };
 
 
    T frobenius_norm() const
    {
        return matrix_->norm();
    };
 
 
    T one_norm() const
    {
        return matrix_->colwise().sum().array().abs().maxCoeff();
    };
 
 
    T infinity_norm() const
    {
        return matrix_->rowwise().sum().array().abs().maxCoeff();
    };
 
 
    Index rank() const
    {
        return matrix_->fullPivLu().rank();
    };
 
};
 
}
 
#endif