From f31ff8893856c6087397941154c11a7f0ee5ab9d Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 11 Dec 2018 14:44:08 +0100
Subject: [PATCH] Implement MatrixMarket support for BCRSMatrix<double> etc
---
dune/istl/matrixmarket.hh | 196 +++++++++++++++++++++--------
dune/istl/test/matrixmarkettest.cc | 107 +++++++++-------
2 files changed, 206 insertions(+), 97 deletions(-)
diff --git a/dune/istl/matrixmarket.hh b/dune/istl/matrixmarket.hh
index d9ba94229..a8f4faaaf 100644
--- a/dune/istl/matrixmarket.hh
+++ b/dune/istl/matrixmarket.hh
@@ -169,13 +169,13 @@ namespace Dune
template<class M>
struct mm_header_printer;
- template<typename T, typename A, int i, int j>
- struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> >
+ template<typename T, typename A>
+ struct mm_header_printer<BCRSMatrix<T,A> >
{
static void print(std::ostream& os)
{
os<<"%%MatrixMarket matrix coordinate ";
- os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
+ os<<mm_numeric_type<typename Imp::BlockTraits<T>::field_type>::str()<<" general"<<std::endl;
}
};
@@ -185,7 +185,7 @@ namespace Dune
static void print(std::ostream& os)
{
os<<"%%MatrixMarket matrix array ";
- os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl;
+ os<<mm_numeric_type<typename Imp::BlockTraits<B>::field_type>::str()<<" general"<<std::endl;
}
};
@@ -220,6 +220,19 @@ namespace Dune
template<class M>
struct mm_block_structure_header;
+ template<typename T, typename A>
+ struct mm_block_structure_header<BlockVector<T,A> >
+ {
+ typedef BlockVector<T,A> M;
+ static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
+
+ static void print(std::ostream& os, const M&)
+ {
+ os<<"% ISTL_STRUCT blocked ";
+ os<<"1 1"<<std::endl;
+ }
+ };
+
template<typename T, typename A, int i>
struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
{
@@ -232,6 +245,19 @@ namespace Dune
}
};
+ template<typename T, typename A>
+ struct mm_block_structure_header<BCRSMatrix<T,A> >
+ {
+ typedef BCRSMatrix<T,A> M;
+ static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
+
+ static void print(std::ostream& os, const M&)
+ {
+ os<<"% ISTL_STRUCT blocked ";
+ os<<"1 1"<<std::endl;
+ }
+ };
+
template<typename T, typename A, int i, int j>
struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
{
@@ -651,19 +677,35 @@ namespace Dune
* @param row The row data as read from file.
* @param matrix The matrix whose data we set.
*/
- template<typename M>
+ template<typename T>
void operator()(const std::vector<std::set<IndexData<D> > >& rows,
- M& matrix)
+ BCRSMatrix<T>& matrix)
{
- for(typename M::RowIterator iter=matrix.begin();
- iter!= matrix.end(); ++iter)
+ static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
+ for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
{
- for(typename M::size_type brow=iter.index()*brows,
+ auto brow=iter.index();
+ for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
+ (*iter)[siter->index] = siter->number;
+ }
+ }
+
+ /**
+ * @brief Sets the matrix values.
+ * @param row The row data as read from file.
+ * @param matrix The matrix whose data we set.
+ */
+ template<typename T>
+ void operator()(const std::vector<std::set<IndexData<D> > >& rows,
+ BCRSMatrix<FieldMatrix<T,brows,bcols> >& matrix)
+ {
+ for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
+ {
+ for (auto brow=iter.index()*brows,
browend=iter.index()*brows+brows;
brow<browend; ++brow)
{
- typedef typename std::set<IndexData<D> >::const_iterator Siter;
- for(Siter siter=rows[brow].begin(), send=rows[brow].end();
+ for (auto siter=rows[brow].begin(), send=rows[brow].end();
siter != send; ++siter)
(*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
}
@@ -694,12 +736,39 @@ namespace Dune
return std::conj(r);
}
- template<typename T, typename A, int brows, int bcols, typename D>
- void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
+ template<typename M>
+ struct mm_multipliers
+ {};
+
+ template<typename B, typename A>
+ struct mm_multipliers<BCRSMatrix<B,A> >
+ {
+ enum {
+ rows = 1,
+ cols = 1
+ };
+ };
+
+ template<typename B, int i, int j, typename A>
+ struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
+ {
+ enum {
+ rows = i,
+ cols = j
+ };
+ };
+
+ template<typename T, typename A, typename D>
+ void readSparseEntries(Dune::BCRSMatrix<T,A>& matrix,
std::istream& file, std::size_t entries,
const MMHeader& mmHeader, const D&)
{
- typedef Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A> Matrix;
+ typedef Dune::BCRSMatrix<T,A> Matrix;
+
+ // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
+ constexpr int brows = mm_multipliers<Matrix>::rows;
+ constexpr int bcols = mm_multipliers<Matrix>::cols;
+
// First path
// store entries together with column index in a separate
// data structure
@@ -810,6 +879,15 @@ namespace Dune
istr >> cols;
}
+ template<typename T, typename A>
+ void mm_read_vector_entries(Dune::BlockVector<T,A>& vector,
+ std::size_t size,
+ std::istream& istr)
+ {
+ for (int i=0; size>0; ++i, --size)
+ istr>>vector[i];
+ }
+
template<typename T, typename A, int entries>
void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
std::size_t size,
@@ -829,8 +907,8 @@ namespace Dune
* @param istr The input stream to read the data from.
* @warning Not all formats are supported!
*/
- template<typename T, typename A, int entries>
- void readMatrixMarket(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
+ template<typename T, typename A>
+ void readMatrixMarket(Dune::BlockVector<T,A>& vector,
std::istream& istr)
{
using namespace MatrixMarketImpl;
@@ -844,29 +922,37 @@ namespace Dune
if(header.type!=array_type)
DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
- std::size_t size=rows/entries;
- if(size*entries!=rows)
- DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
- vector.resize(size);
+ Dune::Hybrid::ifElse(Dune::IsNumber<T>(),
+ [&](auto id){
+ vector.resize(rows);
+ },
+ [&](auto id){ // Assume that T is a FieldVector
+ T dummy;
+ auto blocksize = id(dummy).size();
+ std::size_t size=rows/blocksize;
+ if(size*blocksize!=rows)
+ DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
- istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
+ vector.resize(size);
+ });
+ istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
mm_read_vector_entries(vector, rows, istr);
}
-
/**
* @brief Reads a sparse matrix from a matrix market file.
* @param matrix The matrix to store the data in.
* @param istr The input stream to read the data from.
* @warning Not all formats are supported!
*/
- template<typename T, typename A, int brows, int bcols>
- void readMatrixMarket(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
+ template<typename T, typename A>
+ void readMatrixMarket(Dune::BCRSMatrix<T,A>& matrix,
std::istream& istr)
{
using namespace MatrixMarketImpl;
+ using Matrix = Dune::BCRSMatrix<T,A>;
MMHeader header;
if(!readMatrixMarketBanner(istr, header)) {
@@ -896,46 +982,42 @@ namespace Dune
std::size_t nnz, blockrows, blockcols;
+ // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
+ constexpr int brows = mm_multipliers<Matrix>::rows;
+ constexpr int bcols = mm_multipliers<Matrix>::cols;
+
std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
matrix.setSize(blockrows, blockcols);
- matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise);
+ matrix.setBuildMode(Dune::BCRSMatrix<T,A>::row_wise);
if(header.type==array_type)
DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
- readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
+ readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
}
- template<typename M>
- struct mm_multipliers
- {};
-
- template<typename B, int i, int j, typename A>
- struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
- {
- enum {
- rows = i,
- cols = j
- };
- };
-
- template<typename B, int i, int j>
- void mm_print_entry(const FieldMatrix<B,i,j>& entry,
- typename FieldMatrix<B,i,j>::size_type rowidx,
- typename FieldMatrix<B,i,j>::size_type colidx,
+ // Print a scalar entry
+ template<typename B>
+ void mm_print_entry(const B& entry,
+ std::size_t rowidx,
+ std::size_t colidx,
std::ostream& ostr)
{
- typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
- typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
- for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
- int coli=colidx;
- for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
- ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
- }
+ Hybrid::ifElse(IsNumber<B>(),
+ [&](auto id) {
+ ostr << rowidx << " " << colidx << " " << entry << std::endl;
+ },
+ [&](auto id) {
+ for (auto row=id(entry).begin(); row != id(entry).end(); ++row, ++rowidx) {
+ int coli=colidx;
+ for (auto col = row->begin(); col != row->end(); ++col, ++coli)
+ ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
+ }
+ });
}
// Write a vector entry
@@ -963,6 +1045,12 @@ namespace Dune
std::integral_constant<int,isnumeric>());
}
+ template<typename T, typename A>
+ std::size_t countEntries(const BlockVector<T,A>& vector)
+ {
+ return vector.size();
+ }
+
template<typename T, typename A, int i>
std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
{
@@ -987,8 +1075,8 @@ namespace Dune
std::ostream& ostr,
const std::integral_constant<int,1>&)
{
- ostr<<matrix.N()*mm_multipliers<M>::rows<<" "
- <<matrix.M()*mm_multipliers<M>::cols<<" "
+ ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
+ <<matrix.M()*MatrixMarketImpl::mm_multipliers<M>::cols<<" "
<<countNonZeros(matrix)<<std::endl;
typedef typename M::const_iterator riterator;
@@ -996,8 +1084,8 @@ namespace Dune
for(riterator row=matrix.begin(); row != matrix.end(); ++row)
for(citerator col = row->begin(); col != row->end(); ++col)
// Matrix Market indexing start with 1!
- mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
- col.index()*mm_multipliers<M>::cols+1, ostr);
+ mm_print_entry(*col, row.index()*MatrixMarketImpl::mm_multipliers<M>::rows+1,
+ col.index()*MatrixMarketImpl::mm_multipliers<M>::cols+1, ostr);
}
diff --git a/dune/istl/test/matrixmarkettest.cc b/dune/istl/test/matrixmarkettest.cc
index 7d9b20864..a2483cd0a 100644
--- a/dune/istl/test/matrixmarkettest.cc
+++ b/dune/istl/test/matrixmarkettest.cc
@@ -21,54 +21,40 @@
#include "laplacian.hh"
#endif
-int main(int argc, char** argv)
+template <class Matrix, class Vector>
+int testMatrixMarket(int N)
{
-#if HAVE_MPI
- MPI_Init(&argc, &argv);
- int size;
- MPI_Comm_size(MPI_COMM_WORLD, &size);
-#endif
- const int BS=1;
- int N=100;
-
- if(argc>1)
- N = atoi(argv[1]);
- std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
-
-
- typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
- typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
- typedef Dune::FieldVector<double,BS> VectorBlock;
- typedef Dune::BlockVector<VectorBlock> BVector;
-
#if HAVE_MPI
typedef int GlobalId;
typedef Dune::OwnerOverlapCopyCommunication<GlobalId> Communication;
Communication comm(MPI_COMM_WORLD);
- std::cout<<comm.communicator().rank()<<" "<<comm.communicator().size()<<
- " "<<size<<std::endl;
+ std::cout << comm.communicator().rank() << " " << comm.communicator().size() << std::endl;
int n;
- BCRSMat mat = setupAnisotropic2d<BCRSMat>(N, comm.indexSet(), comm.communicator(), &n, .011);
+ Matrix mat = setupAnisotropic2d<typename Matrix::block_type>(N, comm.indexSet(), comm.communicator(), &n, .011);
#else
- BCRSMat mat;
+ Matrix mat;
setupLaplacian(mat, N);
#endif
- BVector bv(mat.N()), cv(mat.N());
- typedef BVector::iterator VIter;
+ Vector bv(mat.N()), cv(mat.N());
int i=0;
- for(VIter entry=bv.begin(); bv.end() != entry; ++entry) {
- typedef BVector::block_type::iterator SIter;
- for(SIter sentry=entry->begin(); sentry != entry->end(); ++sentry,++i)
- *sentry=i;
- }
+ Dune::Hybrid::ifElse(Dune::IsNumber<typename Vector::block_type>(),
+ [&](auto id) {
+ for(auto entry=bv.begin(); bv.end() != entry; ++entry, ++i)
+ *entry=i;
+ },
+ [&](auto id) {
+ for(auto entry=bv.begin(); bv.end() != entry; ++entry)
+ for (auto sentry=id(entry)->begin(); sentry != id(entry)->end(); ++sentry,++i)
+ *sentry=i;
+ });
#if HAVE_MPI
comm.remoteIndices().rebuild<false>();
comm.copyOwnerToAll(bv,bv);
- Dune::OverlappingSchwarzOperator<BCRSMat,BVector,BVector,Communication> op(mat, comm);
+ Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,Communication> op(mat, comm);
op.apply(bv, cv);
storeMatrixMarket(mat, std::string("testmat"), comm);
storeMatrixMarket(bv, std::string("testvec"), comm, false);
@@ -81,8 +67,8 @@ int main(int argc, char** argv)
storeMatrixMarket(bv, std::string("testvec"));
#endif
- BCRSMat mat1;
- BVector bv1,cv1;
+ Matrix mat1;
+ Vector bv1,cv1;
#if HAVE_MPI
Communication comm1(MPI_COMM_WORLD);
@@ -100,11 +86,9 @@ int main(int argc, char** argv)
++ret;
std::cerr<<"matrix sizes do not match"<<std::endl;
}
- typedef BCRSMat::const_iterator RowIterator;
- typedef BCRSMat::ConstColIterator ColIterator;
- for(RowIterator row=mat.begin(), row1=mat1.begin(); row!=mat.end(); ++row, ++row1)
- for(ColIterator col=row->begin(), col1=row1->begin(); col!= row->end(); ++col, ++col1)
+ for (auto row=mat.begin(), row1=mat1.begin(); row!=mat.end(); ++row, ++row1)
+ for (auto col=row->begin(), col1=row1->begin(); col!= row->end(); ++col, ++col1)
{
if(col.index()!=col1.index()) {
std::cerr <<"Column indices do not match"<<std::endl;
@@ -116,8 +100,8 @@ int main(int argc, char** argv)
}
}
- for(VIter entry=bv.begin(), entry1=bv1.begin(); bv.end() != entry; ++entry, ++entry1)
- if(*entry!=*entry1)
+ for (auto entry=bv.begin(), entry1=bv1.begin(); bv.end() != entry; ++entry, ++entry1)
+ if (*entry!=*entry1)
{
std::cerr<<"written and read vector do not match"<<std::endl;
++ret;
@@ -126,7 +110,7 @@ int main(int argc, char** argv)
cv1.resize(mat1.M());
#if HAVE_MPI
- Dune::OverlappingSchwarzOperator<BCRSMat,BVector,BVector,Communication> op1(mat1, comm1);
+ Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,Communication> op1(mat1, comm1);
op1.apply(bv1, cv1);
if(comm1.indexSet()!=comm.indexSet())
@@ -135,21 +119,58 @@ int main(int argc, char** argv)
++ret;
}
#else
- typedef Dune::MatrixAdapter<BCRSMat,BVector,BVector> Operator;
+ typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
Operator op1(mat1);
op1.apply(bv1, cv1);
#endif
- for(VIter entry=cv.begin(), entry1=cv1.begin(); cv.end() != entry; ++entry, ++entry1)
- if(*entry!=*entry1)
+ for (auto entry=cv.begin(), entry1=cv1.begin(); cv.end() != entry; ++entry, ++entry1)
+ if (*entry!=*entry1)
{
std::cerr<<"computed vectors do not match"<<std::endl;
++ret;
}
+ return ret;
+}
+
+int main(int argc, char** argv)
+{
+#if HAVE_MPI
+ MPI_Init(&argc, &argv);
+ int size;
+ MPI_Comm_size(MPI_COMM_WORLD, &size);
+#endif
+ const int BS=1;
+ int N=2;
+
+ if(argc>1)
+ N = atoi(argv[1]);
+ std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
+
+ // Test scalar matrices and vectors
+ int ret = testMatrixMarket<Dune::BCRSMatrix<double>, Dune::BlockVector<double> >(N);
+
#if HAVE_MPI
if(ret!=0)
MPI_Abort(MPI_COMM_WORLD, ret);
+#endif
+
+ // Test block matrices and vectors with trivial blocks
+ typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
+ typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+ typedef Dune::FieldVector<double,BS> VectorBlock;
+ typedef Dune::BlockVector<VectorBlock> BVector;
+
+ ret = testMatrixMarket<BCRSMat, BVector>(N);
+
+#if HAVE_MPI
+ if(ret!=0)
+ MPI_Abort(MPI_COMM_WORLD, ret);
+#endif
+
+
+#if HAVE_MPI
MPI_Finalize();
#endif
}
--
GitLab