From f542b0d1abb6092fda06548e28bb9c4fa700a1f2 Mon Sep 17 00:00:00 2001
From: Markus Blatt <mblatt@dune-project.org>
Date: Wed, 3 Nov 2010 16:01:52 +0000
Subject: [PATCH] Store row size in redistribute information. Make
redistribution more modular
[[Imported from SVN: r1345]]
---
dune/istl/matrixredistribute.hh | 173 ++++++++++++++++++++------------
1 file changed, 110 insertions(+), 63 deletions(-)
diff --git a/dune/istl/matrixredistribute.hh b/dune/istl/matrixredistribute.hh
index 5b10aac13..a5b1ee505 100644
--- a/dune/istl/matrixredistribute.hh
+++ b/dune/istl/matrixredistribute.hh
@@ -31,6 +31,14 @@ namespace Dune
void resetSetup()
{}
+ void setNoRows(std::size_t size)
+ {}
+
+ std::size_t getRowSize(std::size_t index) const
+ {
+ return -1;
+ }
+
};
#if HAVE_MPI
@@ -123,7 +131,26 @@ namespace Dune
return setup_;
}
+ void reserve(std::size_t size)
+ {}
+
+ std::size_t& getRowSize(std::size_t index)
+ {
+ return rowSize[index];
+ }
+
+ std::size_t getRowSize(std::size_t index) const
+ {
+ return rowSize[index];
+ }
+ void setNoRows(std::size_t rows)
+ {
+ rowSize.resize(rows, 0);
+ }
+
private:
+ std::vector<std::size_t> rowSize;
+
RedistributeInterface interface;
bool setup_;
};
@@ -136,7 +163,7 @@ namespace Dune
* is communicated of.
* @tparam I The type of the index set.
*/
- template<class M>
+ template<class M, class RI>
struct CommMatrixRowSize
{
// Make the default communication policy work.
@@ -148,11 +175,11 @@ namespace Dune
* @param m_ The matrix whose sparsity pattern is communicated.
* @param[out] rowsize_ The vector containing the row sizes
*/
- CommMatrixRowSize(const M& m_, std::vector<size_type>& rowsize_)
+ CommMatrixRowSize(const M& m_, RI& rowsize_)
: matrix(m_), rowsize(rowsize_)
{}
const M& matrix;
- std::vector<size_type>& rowsize;
+ RI& rowsize;
};
@@ -165,7 +192,7 @@ namespace Dune
* is communicated of.
* @tparam I The type of the index set.
*/
- template<class M, class I>
+ template<class M, class I, class RI>
struct CommMatrixSparsityPattern
{
typedef typename M::size_type size_type;
@@ -176,7 +203,7 @@ namespace Dune
* @param idxset_ The index set corresponding to the local matrix.
* @param aggidxset_ The index set corresponding to the redistributed matrix.
*/
- CommMatrixSparsityPattern(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
+ CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_)
: matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize()
{}
@@ -187,8 +214,8 @@ namespace Dune
* @param aggidxset_ The index set corresponding to the redistributed matrix.
* @param rowsize_ The row size for the redistributed owner rows.
*/
- CommMatrixSparsityPattern(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
- const std::vector<typename M::size_type>& rowsize_)
+ CommMatrixSparsityPattern(const M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
+ const RI& rowsize_)
: matrix(m_), idxset(idxset_), aggidxset(aggidxset_), sparsity(aggidxset_.size()), rowsize(&rowsize_)
{}
@@ -255,18 +282,18 @@ namespace Dune
}
}
- M& matrix;
+ const M& matrix;
typedef Dune::GlobalLookupIndexSet<I> LookupIndexSet;
const Dune::GlobalLookupIndexSet<I>& idxset;
const I& aggidxset;
std::vector<std::set<size_type> > sparsity;
- const std::vector<size_type>* rowsize;
+ const RI* rowsize;
};
- template<class M, class I>
- struct CommPolicy<CommMatrixSparsityPattern<M,I> >
+ template<class M, class I, class RI>
+ struct CommPolicy<CommMatrixSparsityPattern<M,I,RI> >
{
- typedef CommMatrixSparsityPattern<M,I> Type;
+ typedef CommMatrixSparsityPattern<M,I,RI> Type;
/**
* @brief The indexed type we send.
@@ -283,8 +310,8 @@ namespace Dune
return t.matrix[i].size();
else
{
- assert((*t.rowsize)[i]>0);
- return (*t.rowsize)[i];
+ assert(t.rowsize->getRowSize(i)>0);
+ return t.rowsize->getRowSize(i);
}
}
};
@@ -295,7 +322,7 @@ namespace Dune
* @tparam M The type of the matrix.
* @tparam I The type of the ParallelIndexSet.
*/
- template<class M, class I>
+ template<class M, class I, class RI>
struct CommMatrixRow
{
/**
@@ -314,7 +341,7 @@ namespace Dune
* @brief Constructor.
*/
CommMatrixRow(M& m_, const Dune::GlobalLookupIndexSet<I>& idxset_, const I& aggidxset_,
- std::vector<size_t>& rowsize_)
+ RI& rowsize_)
: matrix(m_), idxset(idxset_), aggidxset(aggidxset_), rowsize(&rowsize_)
{}
/**
@@ -351,13 +378,13 @@ namespace Dune
/** @brief Index set for the redistributed matrix. */
const I& aggidxset;
/** @brief row size information for the receiving side. */
- std::vector<size_t>* rowsize; // row sizes differ from sender side in ovelap!
+ RI* rowsize; // row sizes differ from sender side in ovelap!
};
- template<class M, class I>
- struct CommPolicy<CommMatrixRow<M,I> >
+ template<class M, class I, class RI>
+ struct CommPolicy<CommMatrixRow<M,I,RI> >
{
- typedef CommMatrixRow<M,I> Type;
+ typedef CommMatrixRow<M,I,RI> Type;
/**
* @brief The indexed type we send.
@@ -374,16 +401,16 @@ namespace Dune
return t.matrix[i].size();
else
{
- assert((*t.rowsize)[i]>0);
- return (*t.rowsize)[i];
+ assert(t.rowsize->getRowSize(i)>0);
+ return t.rowsize->getRowSize(i);
}
}
};
- template<class M, class I>
+ template<class M, class I, class RI>
struct MatrixRowSizeGatherScatter
{
- typedef CommMatrixRowSize<M> Container;
+ typedef CommMatrixRowSize<M,RI> Container;
static const typename M::size_type gather(const Container& cont, std::size_t i)
{
@@ -392,15 +419,15 @@ namespace Dune
static void scatter(Container& cont, const typename M::size_type& rowsize, std::size_t i)
{
assert(rowsize);
- cont.rowsize[i]=rowsize;
+ cont.rowsize.getRowSize(i)=rowsize;
}
};
- template<class M, class I>
+ template<class M, class I, class RI>
struct MatrixSparsityPatternGatherScatter
{
typedef typename I::GlobalIndex GlobalIndex;
- typedef CommMatrixSparsityPattern<M,I> Container;
+ typedef CommMatrixSparsityPattern<M,I,RI> Container;
typedef typename M::ConstColIterator ColIter;
static ColIter col;
@@ -451,14 +478,14 @@ namespace Dune
}
};
- template<class M, class I>
- typename MatrixSparsityPatternGatherScatter<M,I>::ColIter MatrixSparsityPatternGatherScatter<M,I>::col;
+ template<class M, class I, class RI>
+ typename MatrixSparsityPatternGatherScatter<M,I,RI>::ColIter MatrixSparsityPatternGatherScatter<M,I,RI>::col;
- template<class M, class I>
+ template<class M, class I, class RI>
struct MatrixRowGatherScatter
{
typedef typename I::GlobalIndex GlobalIndex;
- typedef CommMatrixRow<M,I> Container;
+ typedef CommMatrixRow<M,I,RI> Container;
typedef typename M::ConstColIterator ColIter;
typedef typename std::pair<GlobalIndex,typename M::block_type> Data;
static ColIter col;
@@ -490,44 +517,31 @@ namespace Dune
}
};
- template<class M, class I>
- typename MatrixRowGatherScatter<M,I>::ColIter MatrixRowGatherScatter<M,I>::col;
+ template<class M, class I, class RI>
+ typename MatrixRowGatherScatter<M,I,RI>::ColIter MatrixRowGatherScatter<M,I,RI>::col;
+
+ template<class M, class I, class RI>
+ typename MatrixRowGatherScatter<M,I,RI>::Data MatrixRowGatherScatter<M,I,RI>::datastore;
+
- template<class M, class I>
- typename MatrixRowGatherScatter<M,I>::Data MatrixRowGatherScatter<M,I>::datastore;
- /**
- * @brief Redistribute a matrix according to given domain decompositions.
- *
- * All the parameters for this function can be obtained by calling
- * graphRepartition with the graph of the original matrix.
- *
- * @param origMatrix The matrix on the original partitioning.
- * @param newMatrix An empty matrix to store the new redistributed matrix in.
- * @param origComm The parallel information of the original partitioning.
- * @param newComm The parallel information of the new partitioning.
- * @param ri The remote index information between the original and the new partitioning.
- * Upon exit of this method it will be prepared for copying from owner to owner vertices
- * for data redistribution.
- * @tparam M The matrix type. It is assumed to be sparse. E.g. BCRSMatrix.
- * @tparam C The type of the parallel information, see OwnerOverlapCopyCommunication.
- */
template<typename M, typename C>
- void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
- RedistributeInformation<C>& ri)
+ void redistributeSparsityPattern(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
+ RedistributeInformation<C>& ri)
{
- std::vector<typename M::size_type> rowsize(newComm.indexSet().size(), 0);
-
typedef typename C::ParallelIndexSet IndexSet;
- CommMatrixRowSize<M> commRowSize(origMatrix, rowsize);
- ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet> >(commRowSize,commRowSize);
+ typedef RedistributeInformation<C> RI;
+ CommMatrixRowSize<M,RI> commRowSize(origMatrix, ri);
+ ri.template redistribute<MatrixRowSizeGatherScatter<M,IndexSet,RI> >(commRowSize,commRowSize);
origComm.buildGlobalLookup();
- CommMatrixSparsityPattern<M,IndexSet> origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
- CommMatrixSparsityPattern<M,IndexSet> newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), rowsize);
+ CommMatrixSparsityPattern<M,IndexSet,RedistributeInformation<C> >
+ origsp(origMatrix, origComm.globalLookup(), newComm.indexSet());
+ CommMatrixSparsityPattern<M,IndexSet,RedistributeInformation<C> >
+ newsp(origMatrix, origComm.globalLookup(), newComm.indexSet(), ri);
- ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet> >(origsp,newsp);
+ ri.template redistribute<MatrixSparsityPatternGatherScatter<M,IndexSet,RI> >(origsp,newsp);
newsp.storeSparsityPattern(newMatrix);
@@ -554,16 +568,49 @@ namespace Dune
if(ret)
DUNE_THROW(ISTLError, "Matrix not symmetric!");
#endif
+ }
- CommMatrixRow<M,IndexSet> origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
- CommMatrixRow<M,IndexSet> newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),rowsize);
+ template<typename M, typename C>
+ void redistributeMatrixEntries(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
+ RedistributeInformation<C>& ri)
+ {
+ typedef typename C::ParallelIndexSet IndexSet;
+ typedef RedistributeInformation<C> RI;
+ CommMatrixRow<M,IndexSet,RedistributeInformation<C> >
+ origrow(origMatrix, origComm.globalLookup(), newComm.indexSet());
+ CommMatrixRow<M,IndexSet,RedistributeInformation<C> >
+ newrow(newMatrix, origComm.globalLookup(), newComm.indexSet(),ri);
- ri.template redistribute<MatrixRowGatherScatter<M,IndexSet> >(origrow,newrow);
+ ri.template redistribute<MatrixRowGatherScatter<M,IndexSet,RI> >(origrow,newrow);
newrow.setOverlapRowsToDirichlet();
if(newMatrix.N()>0&&newMatrix.N()<20)
printmatrix(std::cout, newMatrix, "redist", "row");
}
+ /**
+ * @brief Redistribute a matrix according to given domain decompositions.
+ *
+ * All the parameters for this function can be obtained by calling
+ * graphRepartition with the graph of the original matrix.
+ *
+ * @param origMatrix The matrix on the original partitioning.
+ * @param newMatrix An empty matrix to store the new redistributed matrix in.
+ * @param origComm The parallel information of the original partitioning.
+ * @param newComm The parallel information of the new partitioning.
+ * @param ri The remote index information between the original and the new partitioning.
+ * Upon exit of this method it will be prepared for copying from owner to owner vertices
+ * for data redistribution.
+ * @tparam M The matrix type. It is assumed to be sparse. E.g. BCRSMatrix.
+ * @tparam C The type of the parallel information, see OwnerOverlapCopyCommunication.
+ */
+ template<typename M, typename C>
+ void redistributeMatrix(M& origMatrix, M& newMatrix, C& origComm, C& newComm,
+ RedistributeInformation<C>& ri)
+ {
+ ri.setNoRows(newComm.indexSet().size());
+ redistributeSparsityPattern(origMatrix, newMatrix, origComm, newComm, ri);
+ redistributeMatrixEntries(origMatrix, newMatrix, origComm, newComm, ri);
+ }
#endif
}
#endif
--
GitLab