istlsparsematrixadapter.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM 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 2 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_ISTL_SPARSE_MATRIX_ADAPTER_HH
#define EWOMS_ISTL_SPARSE_MATRIX_ADAPTER_HH

#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>

namespace Opm {
namespace Linear {

template <class MatrixBlockType, class AllocatorType=std::allocator<MatrixBlockType> >
class IstlSparseMatrixAdapter
{
public:
    using IstlMatrix = Dune::BCRSMatrix<MatrixBlockType, AllocatorType>;

    using MatrixBlock = typename IstlMatrix::block_type;
    static_assert(std::is_same<MatrixBlock, MatrixBlockType>::value,
                  "IstlMatrix::block_type and MatrixBlockType must be identical");

    using Scalar = typename MatrixBlock::field_type;

    IstlSparseMatrixAdapter(const size_t rows, const size_t columns)
        : rows_(rows)
        , columns_(columns)
        , istlMatrix_()
    {}

    template <class Simulator>
    explicit IstlSparseMatrixAdapter(const Simulator& simulator)
        : IstlSparseMatrixAdapter(simulator.model().numTotalDof(), simulator.model().numTotalDof())
    {}

    template <class Set>
    void reserve(const std::vector<Set>& sparsityPattern)
    {
        // allocate raw matrix
        istlMatrix_.reset(new IstlMatrix(rows_, columns_, IstlMatrix::random));

        // make sure sparsityPattern is consistent with number of rows
        assert(rows_ == sparsityPattern.size());

        // allocate space for the rows of the matrix
        for (size_t dofIdx = 0; dofIdx < rows_; ++ dofIdx)
            istlMatrix_->setrowsize(dofIdx, sparsityPattern[dofIdx].size());

        istlMatrix_->endrowsizes();

        // fill the rows with indices. each degree of freedom talks to
        // all of its neighbors. (it also talks to itself since
        // degrees of freedom are sometimes quite egocentric.)
        for (size_t dofIdx = 0; dofIdx < rows_; ++ dofIdx) {
            auto nIt    = sparsityPattern[dofIdx].begin();
            auto nEndIt = sparsityPattern[dofIdx].end();
            for (; nIt != nEndIt; ++nIt)
                istlMatrix_->addindex(dofIdx, *nIt);
        }
        istlMatrix_->endindices();
    }

    IstlMatrix& istlMatrix()
    { return *istlMatrix_; }
    const IstlMatrix& istlMatrix() const
    { return *istlMatrix_; }

    size_t rows() const
    { return rows_; }

    size_t cols() const
    { return columns_; }

    void clear()
    { (*istlMatrix_) = Scalar(0.0); }

    void clearRow(const size_t row, const Scalar diag = 1.0)
    {
        MatrixBlock diagBlock(Scalar(0));
        for (int i = 0; i < diagBlock.rows; ++i)
            diagBlock[i][i] = diag;

        auto& matRow = (*istlMatrix_)[row];
        auto colIt = matRow.begin();
        const auto& colEndIt = matRow.end();
        for (; colIt != colEndIt; ++colIt) {
            if (colIt.index() == row)
                *colIt = diagBlock;
            else
                *colIt = Scalar(0.0);
        }
    }

    void block(const size_t rowIdx, const size_t colIdx, MatrixBlock& value) const
    { value = (*istlMatrix_)[rowIdx][colIdx]; }

    MatrixBlockType* blockAddress(const size_t rowIdx, const size_t colIdx) const
    { return &(*istlMatrix_)[rowIdx][colIdx]; }


    void setBlock(const size_t rowIdx, const size_t colIdx, const MatrixBlock& value)
    { (*istlMatrix_)[rowIdx][colIdx] = value; }

    void addToBlock(const size_t rowIdx, const size_t colIdx, const MatrixBlock& value)
    { (*istlMatrix_)[rowIdx][colIdx] += value; }

    void commit()
    { }

    void finalize()
    { }

protected:
    size_t rows_;
    size_t columns_;

    std::unique_ptr<IstlMatrix> istlMatrix_;
};

}} // namespace Linear, Opm

#endif