LisMatrix.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 <string>
#include <vector>
 
#include "LisCheck.h"
#include "LisWrapper.h"
#include "MathLib/LinAlg/RowColumnIndices.h"
#include "MathLib/LinAlg/SetMatrixSparsity.h"
 
namespace MathLib
{
// Forward declarations.
class LisVector;
class LisMatrix;
 
template <typename SPARSITY_PATTERN>
struct SetMatrixSparsity<LisMatrix, SPARSITY_PATTERN>;

class LisMatrix
{
public:
    enum class MatrixType : int
    {
        CRS = 1,
        CCS = 2,
        MSR = 3,
        DIA = 4,
        ELL = 5,
        JDS = 6,
        BSR = 7,
        BSC = 8,
        VBR = 9,
        COO = 10,
        DNS = 11
    };
 
    using IndexType = LIS_INT;
 
public:
    LisMatrix(std::size_t n_rows, MatrixType mat_type = MatrixType::CRS);

    LisMatrix(std::size_t n_rows, int nonzero, IndexType* row_ptr,
              IndexType* col_idx, double* data);

    virtual ~LisMatrix();

    std::size_t getNumberOfRows() const { return n_rows_; }

    std::size_t getNumberOfColumns() const { return getNumberOfRows(); }

    std::size_t getRangeBegin() const { return is_; }

    std::size_t getRangeEnd() const { return ie_; }

    void setZero();

    int setValue(IndexType rowId, IndexType colId, double v);

    int add(IndexType rowId, IndexType colId, double v);

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

    LIS_MATRIX& getRawMatrix() { return AA_; }

    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);

    MatrixType getMatrixType() const { return mat_type_; }

    bool isAssembled() const { return is_assembled_; }
 
private:
    std::size_t const n_rows_;
    MatrixType const mat_type_;
    LIS_MATRIX AA_;
    LIS_VECTOR diag_;
    bool is_assembled_;
    IndexType is_;  
    IndexType
        ie_;  
    bool use_external_arrays_;
 
    // friend function
    friend bool finalizeMatrixAssembly(LisMatrix& mat);
 
    template <typename MATRIX, typename SPARSITY_PATTERN>
    friend struct SetMatrixSparsity;
};
 
template <class T_DENSE_MATRIX>
void LisMatrix::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));
        }
    }
};

bool finalizeMatrixAssembly(LisMatrix& mat);

template <typename SPARSITY_PATTERN>
struct SetMatrixSparsity<LisMatrix, SPARSITY_PATTERN>
{
    void operator()(LisMatrix& matrix,
                    SPARSITY_PATTERN const& sparsity_pattern) const
    {
        auto const n_rows = matrix.getNumberOfRows();
        std::vector<LisMatrix::IndexType> row_sizes;
        row_sizes.reserve(n_rows);
 
        // LIS needs 1 more entry, otherwise it starts reallocating arrays.
        transform(cbegin(sparsity_pattern), cend(sparsity_pattern),
                  back_inserter(row_sizes), [](auto const i) { return i + 1; });
 
        int ierr = lis_matrix_malloc(matrix.AA_, 0, row_sizes.data());
        checkLisError(ierr);
    }
};
 
}  // namespace MathLib