eigen_sparse.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_SPARSE_MATRIX_HPP
#define EIGEN_SPARSE_MATRIX_HPP
 
#include <vector>
#include <numeric>
#include <stdexcept>
 
#include <external/Eigen/Sparse>
#include <external/Eigen/SparseLU>
#include <external/Eigen/IterativeLinearSolvers>
#include <external/Eigen/SVD>
#include <external/EigenUnsupported/Eigen/IterativeSolvers>
 
#include "types.hpp"
#include "matrix/base.hpp"
#include "matrix/eigen_base.hpp"
#include "matrix/eigen_dense.hpp"
 
 
namespace bem
{
 
template <typename T, int StorageOrder = Eigen::ColMajor>
class EigenSparseMatrix: public EigenMatrixBase<T, Eigen::SparseMatrix<T, StorageOrder>>
{
 
    using MatrixType = Eigen::SparseMatrix<T, StorageOrder>;
    using base = EigenMatrixBase<T, MatrixType>;
    using base::base;
    using base::matrix_;
 
public:
 
    using base::raw_matrix;
 
    void preallocate(const std::vector<Index>& nnz) override
    {
        triplets_.reserve(std::accumulate(nnz.begin(), nnz.end(), 0));
        assembled_ = false;
        return;
    };
 
 
    void preallocate(Index num_entries) override
    {
        triplets_.reserve(num_entries);
        assembled_ = false;
        return;
    };
 
 
    void preallocate_directly(const std::vector<Index>& nnz)
    {
        matrix_->reserve(nnz);
        assembled_ = false;
        return;
    };
 
 
    void assemble() override
    {
        if (assembled_)
            return;
 
        triplets_.shrink_to_fit();
        if (insert_mode_on_)
            matrix_->setFromTriplets(
                triplets_.begin(),
                triplets_.end(),
                [] (const T&, const T &b) { return b; }
                );
        else
            matrix_->setFromTriplets(triplets_.begin(), triplets_.end());
        triplets_.clear();
        assembled_ = true;
        return;
    };
 
 
    T value(Index row, Index col) const override
    {
        return matrix_->coeff(row, col);
    };
 
 
    void set_value(Index row, Index col, const T& a) override
    {
        triplets_.push_back(Eigen::Triplet<T> (row, col, a));
        insert_mode_on_ = true;
        assembled_ = false;
        return;
    };
 
 
    void set_value_directly(Index row, Index col, const T& a)
    {
        matrix_->coeffRef(row, col) = a;
        matrix_->makeCompressed();
        return;
    };
 
 
    void add_value(Index row, Index col, const T& a) override
    {
        triplets_.push_back(Eigen::Triplet<T> (row, col, a));
        insert_mode_on_ = false;
        assembled_ = false;
        return;
    };
 
 
    void add_value_directly(Index row, Index col, const T& a)
    {
        matrix_->coeffRef(row, col) += a;
        return;
    };
 
 
    void scale(const T& a) override
    {
        for (Index kk = 0; kk < matrix_->outerSize(); ++kk)
            for (typename MatrixType::InnerIterator it((*matrix_), kk); it; ++it)
                it.valueRef() = a * it.value();
        matrix_->makeCompressed();
        return;
    };
 
 
    void factorize()
    {
        solver_.compute((*matrix_));
        if (solver_.info() != Eigen::Success)
            throw std::runtime_error("Matrix decomposition failed.");
        factorized_ = true;
        return;
    };
 
 
    void mat_solve(MatrixBase<T>& x, const MatrixBase<T>& b) override
    {
        if (!factorized_)
            factorize();
 
        const EigenDenseMatrix<T>& bd = dynamic_cast<const EigenDenseMatrix<T>&> (b);
        EigenDenseMatrix<T>& xd = dynamic_cast<EigenDenseMatrix<T>&> (x);
 
        xd.raw_matrix() = solver_.solve(bd.raw_matrix());
 
        if (solver_.info() != Eigen::Success)
            throw std::runtime_error("Matrix solve failed.");
 
        return;
    };
 
 
    void set_block(
        const MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        const T& a = T(1)
        ) override
    {
        const EigenSparseMatrix<T>& xd = dynamic_cast<const EigenSparseMatrix<T>&> (x);
 
        std::vector<Eigen::Triplet<T>> x_triplets;
        x_triplets.reserve(xd.raw_matrix().nonZeros());
 
        for (Index kk = 0; kk < xd.raw_matrix().outerSize(); ++kk)
            for (typename MatrixType::InnerIterator it (xd.raw_matrix(), kk); it; ++it)
                x_triplets.push_back(Eigen::Triplet<T> (row_start + it.row(), col_start + it.col(), it.value() * a));
 
        matrix_->insertFromTriplets(
            x_triplets.begin(),
            x_triplets.end(),
            [] (const T&, const T &b) { return b; }
            );
 
        return;
    };
 
 
    void add_block(
        const MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        const T& a = T(1)
        ) override
    {
        const EigenSparseMatrix<T>& xd = dynamic_cast<const EigenSparseMatrix<T>&> (x);
 
        std::vector<Eigen::Triplet<T>> x_triplets;
        x_triplets.reserve(xd.raw_matrix().nonZeros());
 
        for (Index kk = 0; kk < xd.raw_matrix().outerSize(); ++kk)
            for (typename MatrixType::InnerIterator it (xd.raw_matrix(), kk); it; ++it)
                x_triplets.push_back(Eigen::Triplet<T> (row_start + it.row(), col_start + it.col(), it.value() * a));
 
        matrix_->insertFromTriplets(x_triplets.begin(), x_triplets.end());
 
        return;
    };
 
 
    void get_block(
        MatrixBase<T>& x,
        Index row_start,
        Index col_start,
        Index b_rows,
        Index b_cols
        ) const override
    {
        EigenSparseMatrix<T>& xd = dynamic_cast<EigenSparseMatrix<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_->valuePtr(); };
 
 
    virtual T* data() override
    { return matrix_->valuePtr(); };
 
 
private:
 
    std::vector<Eigen::Triplet<T>> triplets_;
    bool insert_mode_on_ = true;
    bool assembled_ = false;
 
    Eigen::SparseLU<MatrixType> solver_;
    bool factorized_ = false;
 
};
 
}
 
#endif