From 389a9d44a4157c88cc527fb54e3327287d7fe8ad Mon Sep 17 00:00:00 2001
From: Jonathan Youett <youett@math.fu-berlin.de>
Date: Tue, 8 Nov 2016 13:37:36 +0100
Subject: [PATCH] Add a blockmatrix view that computes global indices from
block-indices
This class also provides a static method to combine submatrices into
one block matrix by copying the entries.
---
dune/matrix-vector/CMakeLists.txt | 1 +
dune/matrix-vector/blockmatrixview.hh | 109 ++++++++++++++++++++++++++
2 files changed, 110 insertions(+)
create mode 100644 dune/matrix-vector/blockmatrixview.hh
diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt
index 9acfb4a..9112065 100644
--- a/dune/matrix-vector/CMakeLists.txt
+++ b/dune/matrix-vector/CMakeLists.txt
@@ -5,6 +5,7 @@ install(FILES
addtodiagonal.hh
axpy.hh
axy.hh
+ blockmatrixview.hh
componentwisematrixmap.hh
crossproduct.hh
genericvectortools.hh
diff --git a/dune/matrix-vector/blockmatrixview.hh b/dune/matrix-vector/blockmatrixview.hh
new file mode 100644
index 0000000..54af5ae
--- /dev/null
+++ b/dune/matrix-vector/blockmatrixview.hh
@@ -0,0 +1,109 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+#ifndef DUNE_MATRIX_VECTOR_BLOCK_MATRIX_VIEW_HH
+#define DUNE_MATRIX_VECTOR_BLOCK_MATRIX_VIEW_HH
+
+#include <array>
+#include <vector>
+
+#include <dune/istl/matrixindexset.hh>
+
+namespace Dune {
+namespace MatrixVector {
+
+/** \brief Class that provides global indices for a block matrix given local indices of the blocks. Furthermore, it provides a method to construct a block matrix from given sub matrices.
+ *
+ * This class provides a method to combine matrices \f A_1,\ldots,A_n\f into one single block matrix
+ * \f B = \begin{pmatrix} A_1 & 0 & \ldots & 0 \\
+ * 0 & A_2 & \ldots & 0 \\
+ * \ldots \\
+ * 0 & \ldots & 0 & A_n \end{pmatrix} \f,
+ * by copying the corresponding entries.
+ * A BlockMatrixView object provides methods that return global indices for \f B\f given a local index of \f A_i\f and the block index \f i\f.
+ *
+ * \tparam MatrixType The matrix type of the matrix blocks.
+ */
+template <class MatrixType>
+class BlockMatrixView
+{
+public:
+ //! Construct block view from vector of matrices corr. to the blocks
+ BlockMatrixView(const std::vector<const MatrixType*>& submat) {
+
+ offsets_.resize(submat.size());
+ offsets_[0] = {0, 0};
+
+ for (size_t i=0; i<submat.size()-1; i++)
+ offsets_[i+1] = {offsets_[i][0] + submat[i]->N(),
+ offsets_[i][1] + submat[i]->M()};
+
+ nRows_ = offsets_.back()[0] + submat.back()->N();
+ nCols_ = offsets_.back()[1] + submat.back()->M();
+ }
+
+ //! Get global row index for given block and a local row index
+ size_t row(size_t block, size_t index) const {
+ return offsets_[block][0] + index;
+ }
+
+ //! Get global column index given a block and a local column index
+ size_t col(size_t block, size_t index) const {
+ return offsets_[block][1] + index;
+ }
+
+ //! Return the total number of rows
+ size_t nRows() const {
+ return nRows_;
+ }
+
+ //! Return the total number of columns
+ size_t nCols() const {
+ return nCols_;
+ }
+
+ /** \brief Combine submatrices in a big matrix. */
+ static void setupBlockMatrix(const std::vector<MatrixType>& submat, MatrixType& mat) {
+ std::vector<const MatrixType*> dummy(submat.size());
+ for (size_t i=0; i<dummy.size(); i++)
+ dummy[i] = &submat[i];
+ setupBlockMatrix(dummy,mat);
+ }
+
+ /** \brief Combine submatrices in a big matrix. */
+ static void setupBlockMatrix(const std::vector<const MatrixType*>& submat, MatrixType& mat)
+ {
+ BlockMatrixView<MatrixType> blockView(submat);
+ MatrixIndexSet indexSet(blockView.nRows(), blockView.nCols());
+
+ for (size_t i=0; i<submat.size(); i++)
+ indexSet.import(*submat[i], blockView.row(i,0), blockView.col(i,0));
+
+ indexSet.exportIdx(mat);
+
+ for (size_t i=0; i<submat.size(); i++) {
+
+ auto row = submat[i]->begin();
+ auto rowEnd = submat[i]->end();
+
+ for (; row != rowEnd; row++) {
+
+ auto col = row->begin();
+ auto colEnd = row->end();
+
+ for (; col != colEnd; col++)
+ mat[blockView.row(i,row.index())][blockView.col(i,col.index())] = *col;
+ }
+ }
+ }
+
+protected:
+ //! Offsets for the different blocks
+ std::vector<std::array<size_t, 2> > offsets_;
+ //! Total numbers of rows and columns
+ size_t nRows_, nCols_;
+};
+
+} // end namespace MatrixVector
+} // end namespace Dune
+
+#endif
--
GitLab