EigenMatrix.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <Eigen/Sparse>
#include <cassert>
#include <iosfwd>
#include <string>
 
#include "MathLib/LinAlg/RowColumnIndices.h"
#include "MathLib/LinAlg/SetMatrixSparsity.h"
 
namespace MathLib
{
class EigenMatrix final
{
public:
    using RawMatrixType = Eigen::SparseMatrix<double, Eigen::RowMajor>;
    using IndexType = RawMatrixType::Index;
 
    // TODO The matrix constructor should take num_rows and num_cols as
    // arguments
    //      that is left for a later refactoring.

    explicit EigenMatrix(IndexType const n,
                         IndexType const n_nonzero_columns = 0)
        : mat_(n, n)
    {
        if (n_nonzero_columns > 0)
        {
            mat_.reserve(Eigen::Matrix<IndexType, Eigen::Dynamic, 1>::Constant(
                n, n_nonzero_columns));
        }
    }

    IndexType getNumberOfRows() const { return mat_.rows(); }

    IndexType getNumberOfColumns() const { return mat_.cols(); }

    static constexpr IndexType getRangeBegin() { return 0; }

    IndexType getRangeEnd() const { return getNumberOfRows(); }

    void setZero()
    {
        auto const N = mat_.nonZeros();
        for (auto i = decltype(N){0}; i < N; i++)
        {
            mat_.valuePtr()[i] = 0;
        }
        // don't use mat_.setZero(). it makes a matrix uncompressed
    }

    int setValue(IndexType row, IndexType col, double val)
    {
        assert(row < (IndexType)getNumberOfRows() &&
               col < (IndexType)getNumberOfColumns());
        mat_.coeffRef(row, col) = val;
        return 0;
    }

    int add(IndexType row, IndexType col, double val)
    {
        mat_.coeffRef(row, col) += val;
        return 0;
    }

    void addToDiagonal(double const value);

    template <class T_DENSE_MATRIX>
    void add(std::vector<IndexType> const& row_pos,
             const T_DENSE_MATRIX& sub_matrix,
             double fkt = 1.0)
    {
        this->add(row_pos, row_pos, sub_matrix, fkt);
    }

    template <class T_DENSE_MATRIX>
    void add(RowColumnIndices<IndexType> const& indices,
             const T_DENSE_MATRIX& sub_matrix,
             double fkt = 1.0)
    {
        this->add(indices.rows, indices.columns, sub_matrix, fkt);
    }

    template <class T_DENSE_MATRIX>
    void add(std::vector<IndexType> const& row_pos,
             std::vector<IndexType> const& col_pos,
             const T_DENSE_MATRIX& sub_matrix, double fkt = 1.0);

    double get(IndexType row, IndexType col) const
    {
        return mat_.coeff(row, col);
    }

    static constexpr bool isAssembled() { return true; }

    void write(const std::string& filename) const;

    void write(std::ostream& os) const;
 
    RawMatrixType& getRawMatrix() { return mat_; }
    const RawMatrixType& getRawMatrix() const { return mat_; }
 
protected:
    RawMatrixType mat_;
};
 
template <class T_DENSE_MATRIX>
void EigenMatrix::add(std::vector<IndexType> const& row_pos,
                      std::vector<IndexType> const& col_pos,
                      const T_DENSE_MATRIX& sub_matrix, double fkt)
{
    auto const n_rows = row_pos.size();
    auto const n_cols = col_pos.size();
    for (auto i = decltype(n_rows){0}; i < n_rows; i++)
    {
        auto const row = row_pos[i];
        for (auto j = decltype(n_cols){0}; j < n_cols; j++)
        {
            auto const col = col_pos[j];
            add(row, col, fkt * sub_matrix(i, j));
        }
    }
};

template <typename SPARSITY_PATTERN>
struct SetMatrixSparsity<EigenMatrix, SPARSITY_PATTERN>
{
    void operator()(EigenMatrix& matrix,
                    SPARSITY_PATTERN const& sparsity_pattern) const
    {
        static_assert(EigenMatrix::RawMatrixType::IsRowMajor,
                      "Set matrix sparsity relies on the EigenMatrix to be in "
                      "row-major storage order.");
 
        assert(matrix.getNumberOfRows() ==
               static_cast<EigenMatrix::IndexType>(
                   sparsity_pattern.numberOfRows()));
 
        matrix.getRawMatrix().reserve(
            sparsity_pattern.number_non_zeros_per_row);
    }
};
 
}  // end namespace MathLib