3#ifndef DUNE_ISTL_MATRIXMARKET_HH
4#define DUNE_ISTL_MATRIXMARKET_HH
25#include <dune/common/hybridutilities.hh>
62 namespace MatrixMarketImpl
84 struct mm_numeric_type<int>
93 static std::string str()
100 struct mm_numeric_type<double>
109 static std::string str()
116 struct mm_numeric_type<float>
125 static std::string str()
132 struct mm_numeric_type<
std::complex<double> >
141 static std::string str()
148 struct mm_numeric_type<
std::complex<float> >
157 static std::string str()
174 template<
typename T,
typename A>
177 static void print(std::ostream& os)
179 os<<
"%%MatrixMarket matrix coordinate ";
180 os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<T>::field_type>>::str()<<
" general"<<std::endl;
184 template<
typename B,
typename A>
187 static void print(std::ostream& os)
189 os<<
"%%MatrixMarket matrix array ";
190 os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<B>::field_type>>::str()<<
" general"<<std::endl;
194 template<
typename T,
int j>
195 struct mm_header_printer<FieldVector<T,j> >
197 static void print(std::ostream& os)
199 os<<
"%%MatrixMarket matrix array ";
200 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
204 template<
typename T,
int i,
int j>
205 struct mm_header_printer<FieldMatrix<T,i,j> >
207 static void print(std::ostream& os)
209 os<<
"%%MatrixMarket matrix array ";
210 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
225 template<
typename T,
typename A>
231 static void print(std::ostream& os,
const M&)
233 os<<
"% ISTL_STRUCT blocked ";
234 os<<
"1 1"<<std::endl;
238 template<
typename T,
typename A,
int i>
243 static void print(std::ostream& os,
const M&)
245 os<<
"% ISTL_STRUCT blocked ";
246 os<<i<<
" "<<1<<std::endl;
250 template<
typename T,
typename A>
251 struct mm_block_structure_header<BCRSMatrix<T,A> >
253 typedef BCRSMatrix<T,A> M;
254 static_assert(IsNumber<T>::value,
"Only scalar entries are expected here!");
256 static void print(std::ostream& os,
const M&)
258 os<<
"% ISTL_STRUCT blocked ";
259 os<<
"1 1"<<std::endl;
263 template<
typename T,
typename A,
int i,
int j>
264 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
266 typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M;
268 static void print(std::ostream& os,
const M&)
270 os<<
"% ISTL_STRUCT blocked ";
271 os<<i<<
" "<<j<<std::endl;
276 template<
typename T,
int i,
int j>
277 struct mm_block_structure_header<FieldMatrix<T,i,j> >
279 typedef FieldMatrix<T,i,j> M;
281 static void print(std::ostream& os,
const M& m)
285 template<
typename T,
int i>
286 struct mm_block_structure_header<FieldVector<T,i> >
288 typedef FieldVector<T,i> M;
290 static void print(std::ostream& os,
const M& m)
294 enum LineType { MM_HEADER, MM_ISTLSTRUCT, DATA };
295 enum { MM_MAX_LINE_LENGTH=1025 };
297 enum MM_TYPE { coordinate_type, array_type, unknown_type };
299 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
301 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
306 : type(coordinate_type), ctype(double_type), structure(general)
310 MM_STRUCTURE structure;
313 inline bool lineFeed(std::istream& file)
337 inline void skipComments(std::istream& file)
351 inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
360 dverb<<buffer<<std::endl;
362 if(buffer!=
"%%MatrixMarket") {
375 if(buffer !=
"matrix")
392 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
399 if(buffer !=
"array")
404 mmHeader.type=array_type;
408 if(buffer !=
"coordinate")
413 mmHeader.type=coordinate_type;
430 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
436 if(buffer !=
"integer")
441 mmHeader.ctype=integer_type;
450 mmHeader.ctype=double_type;
454 if(buffer !=
"complex")
459 mmHeader.ctype=complex_type;
463 if(buffer !=
"pattern")
468 mmHeader.ctype=pattern;
480 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
486 if(buffer !=
"general")
491 mmHeader.structure=general;
495 if(buffer !=
"hermitian")
500 mmHeader.structure=hermitian;
503 if(buffer.size()==1) {
512 if(buffer !=
"symmetric")
517 mmHeader.structure=symmetric;
521 if(buffer !=
"skew-symmetric")
526 mmHeader.structure=skew_symmetric;
543 template<std::
size_t brows, std::
size_t bcols>
544 std::tuple<std::size_t, std::size_t, std::size_t>
545 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries,
const MMHeader& header)
547 std::size_t blockrows=rows/brows;
548 std::size_t blockcols=cols/bcols;
549 std::size_t blocksize=brows*bcols;
550 std::size_t blockentries=0;
552 switch(header.structure)
555 blockentries = entries/blocksize;
break;
556 case skew_symmetric :
557 blockentries = 2*entries/blocksize;
break;
559 blockentries = (2*entries-rows)/blocksize;
break;
561 blockentries = (2*entries-rows)/blocksize;
break;
565 return std::make_tuple(blockrows, blockcols, blockentries);
575 struct IndexData :
public T
612 std::istream&
operator>>(std::istream& is, NumericWrapper<T>& num)
614 return is>>num.number;
617 inline std::istream&
operator>>(std::istream& is, [[maybe_unused]] NumericWrapper<PatternDummy>& num)
628 bool operator<(
const IndexData<T>& i1,
const IndexData<T>& i2)
630 return i1.index<i2.index;
639 std::istream&
operator>>(std::istream& is, IndexData<T>& data)
644 return is>>data.number;
663 data.number = {real.number, imag.number};
673 template<
typename D,
int brows,
int bcols>
682 void operator()(
const std::vector<std::set<IndexData<D> > >& rows,
685 static_assert(
IsNumber<T>::value && brows==1 && bcols==1,
"Only scalar entries are expected here!");
686 for (
auto iter=matrix.
begin(); iter!= matrix.
end(); ++iter)
688 auto brow=iter.index();
689 for (
auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
690 (*iter)[siter->index] = siter->number;
700 void operator()(
const std::vector<std::set<IndexData<D> > >& rows,
703 for (
auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
705 for (
auto brow=iter.index()*brows,
706 browend=iter.index()*brows+brows;
707 brow<browend; ++brow)
709 for (
auto siter=rows[brow].begin(), send=rows[brow].end();
710 siter != send; ++siter)
711 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
717 template<
int brows,
int bcols>
718 struct MatrixValuesSetter<PatternDummy,brows,bcols>
721 void operator()(
const std::vector<std::set<IndexData<PatternDummy> > >& rows,
726 template<
class T>
struct is_complex : std::false_type {};
727 template<
class T>
struct is_complex<
std::complex<T>> : std::true_type {};
731 std::enable_if_t<!is_complex<T>::value, T> conj(
const T& r){
736 std::enable_if_t<is_complex<T>::value, T> conj(
const T& r){
741 struct mm_multipliers
744 template<
typename B,
typename A>
745 struct mm_multipliers<BCRSMatrix<B,A> >
753 template<
typename B,
int i,
int j,
typename A>
754 struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
762 template<
typename T,
typename A,
typename D>
764 std::istream& file, std::size_t entries,
765 const MMHeader& mmHeader,
const D&)
770 constexpr int brows = mm_multipliers<Matrix>::rows;
771 constexpr int bcols = mm_multipliers<Matrix>::cols;
776 std::vector<std::set<IndexData<D> > > rows(matrix.
N()*brows);
778 auto readloop = [&] (
auto symmetryFixup) {
779 for(std::size_t i = 0; i < entries; ++i) {
785 assert(row/bcols<matrix.
N());
787 assert(data.index/bcols<matrix.
M());
788 rows[row].insert(data);
790 symmetryFixup(row, data);
794 switch(mmHeader.structure)
797 readloop([](
auto...){});
800 readloop([&](
auto row,
auto data) {
801 IndexData<D> data_sym(data);
802 data_sym.index = row;
803 rows[data.index].insert(data_sym);
806 case skew_symmetric :
807 readloop([&](
auto row,
auto data) {
808 IndexData<D> data_sym;
809 data_sym.number = -data.number;
810 data_sym.index = row;
811 rows[data.index].insert(data_sym);
815 readloop([&](
auto row,
auto data) {
816 IndexData<D> data_sym;
817 data_sym.number = conj(data.number);
818 data_sym.index = row;
819 rows[data.index].insert(data_sym);
824 "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
829 for(
typename Matrix::CreateIterator iter=matrix.
createbegin();
832 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
833 brow<browend; ++brow)
835 typedef typename std::set<IndexData<D> >::const_iterator Siter;
836 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
837 siter != send; ++siter, ++nnz)
838 iter.insert(siter->index/bcols);
845 MatrixValuesSetter<D,brows,bcols> Setter;
847 Setter(rows, matrix);
855 inline void mm_read_header(std::size_t& rows, std::size_t& cols,
856 MatrixMarketImpl::MMHeader& header, std::istream& istr,
859 using namespace MatrixMarketImpl;
861 if(!readMatrixMarketBanner(istr, header)) {
862 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
863 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
866 istr.seekg(0, std::ios::beg);
868 header.type=array_type;
874 throw MatrixMarketFormatError();
879 throw MatrixMarketFormatError();
883 template<
typename T,
typename A>
889 for (
int i=0; size>0; ++i, --size)
893 template<
typename T,
typename A,
int entries>
899 for(
int i=0; size>0; ++i, --size) {
913 template<
typename T,
typename A>
918 using namespace MatrixMarketImpl;
921 std::size_t rows, cols;
922 mm_read_header(rows,cols,header,istr,
true);
923 if(cols!=Simd::lanes<field_type>()) {
924 if(Simd::lanes<field_type>() == 1)
925 DUNE_THROW(MatrixMarketFormatError,
"cols!=1, therefore this is no vector!");
927 DUNE_THROW(MatrixMarketFormatError,
"cols does not match the number of lanes in the field_type!");
930 if(header.type!=array_type)
931 DUNE_THROW(MatrixMarketFormatError,
"Vectors have to be stored in array format!");
939 auto blocksize = dummy.size();
940 std::size_t size=rows/blocksize;
941 if(size*blocksize!=rows)
942 DUNE_THROW(MatrixMarketFormatError,
"Block size of vector is not correct!");
948 for(
size_t l=0;l<Simd::lanes<field_type>();++l){
949 mm_read_vector_entries(vector, rows, istr, l);
959 template<
typename T,
typename A>
963 using namespace MatrixMarketImpl;
967 if(!readMatrixMarketBanner(istr, header)) {
968 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
969 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
972 istr.seekg(0, std::ios::beg);
976 std::size_t rows, cols, entries;
979 throw MatrixMarketFormatError();
984 throw MatrixMarketFormatError();
988 throw MatrixMarketFormatError();
992 std::size_t nnz, blockrows, blockcols;
995 constexpr int brows = mm_multipliers<Matrix>::rows;
996 constexpr int bcols = mm_multipliers<Matrix>::cols;
998 std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
1003 matrix.
setSize(blockrows, blockcols);
1006 if(header.type==array_type)
1009 readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1013 template<
typename B>
1014 void mm_print_entry(
const B& entry,
1019 if constexpr (IsNumber<B>())
1020 ostr << rowidx <<
" " << colidx <<
" " << entry << std::endl;
1023 for (
auto row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
1025 for (
auto col = row->begin(); col != row->end(); ++col, ++coli)
1026 ostr<< rowidx<<
" "<<coli<<
" "<<*col<<std::endl;
1032 template<
typename V>
1033 void mm_print_vector_entry(
const V& entry, std::ostream& ostr,
1034 const std::integral_constant<int,1>&,
1041 template<
typename V>
1042 void mm_print_vector_entry(
const V& vector, std::ostream& ostr,
1043 const std::integral_constant<int,0>&,
1046 using namespace MatrixMarketImpl;
1049 const int isnumeric = mm_numeric_type<Simd::Scalar<typename V::block_type>>::is_numeric;
1050 typedef typename V::const_iterator VIter;
1052 for(VIter i=vector.begin(); i != vector.end(); ++i)
1054 mm_print_vector_entry(*i, ostr,
1055 std::integral_constant<int,isnumeric>(),
1059 template<
typename T,
typename A>
1060 std::size_t countEntries(
const BlockVector<T,A>& vector)
1062 return vector.size();
1065 template<
typename T,
typename A,
int i>
1066 std::size_t countEntries(
const BlockVector<FieldVector<T,i>,A>& vector)
1068 return vector.size()*i;
1072 template<
typename V>
1073 void writeMatrixMarket(
const V& vector, std::ostream& ostr,
1074 const std::integral_constant<int,0>&)
1076 using namespace MatrixMarketImpl;
1077 typedef typename V::field_type field_type;
1079 ostr<<countEntries(vector)<<
" "<<Simd::lanes<field_type>()<<std::endl;
1080 const int isnumeric = mm_numeric_type<Simd::Scalar<V>>::is_numeric;
1081 for(
size_t l=0;l<Simd::lanes<field_type>(); ++l){
1082 mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>(), l);
1087 template<
typename M>
1088 void writeMatrixMarket(
const M& matrix,
1090 const std::integral_constant<int,1>&)
1092 ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<
" "
1093 <<matrix.M()*MatrixMarketImpl::mm_multipliers<M>::cols<<
" "
1096 typedef typename M::const_iterator riterator;
1097 typedef typename M::ConstColIterator citerator;
1098 for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1099 for(citerator col = row->begin(); col != row->end(); ++col)
1101 mm_print_entry(*col, row.index()*MatrixMarketImpl::mm_multipliers<M>::rows+1,
1102 col.index()*MatrixMarketImpl::mm_multipliers<M>::cols+1, ostr);
1109 template<
typename M>
1110 void writeMatrixMarket(
const M& matrix,
1113 using namespace MatrixMarketImpl;
1116 mm_header_printer<M>::print(ostr);
1117 mm_block_structure_header<M>::print(ostr,matrix);
1122 static const int default_precision = -1;
1134 template<
typename M>
1136 std::string filename,
1137 int prec=default_precision)
1139 std::ofstream file(filename.c_str());
1140 file.setf(std::ios::scientific,std::ios::floatfield);
1142 file.precision(prec);
1143 writeMatrixMarket(matrix, file);
1162 template<
typename M,
typename G,
typename L>
1164 std::string filename,
1166 bool storeIndices=
true,
1167 int prec=default_precision)
1170 int rank = comm.communicator().rank();
1172 std::ostringstream rfilename;
1173 rfilename<<filename <<
"_"<<rank<<
".mm";
1174 dverb<< rfilename.str()<<std::endl;
1175 std::ofstream file(rfilename.str().c_str());
1176 file.setf(std::ios::scientific,std::ios::floatfield);
1178 file.precision(prec);
1179 writeMatrixMarket(matrix, file);
1187 rfilename<<filename<<
"_"<<rank<<
".idx";
1188 file.open(rfilename.str().c_str());
1189 file.setf(std::ios::scientific,std::ios::floatfield);
1191 typedef typename IndexSet::const_iterator Iterator;
1194 file << iter->global()<<
" "<<(std::size_t)iter->local()<<
" "
1195 <<(int)iter->local().attribute()<<
" "<<(int)iter->local().isPublic()<<std::endl;
1198 file<<
"neighbours:";
1199 const std::set<int>& neighbours=comm.
remoteIndices().getNeighbours();
1200 typedef std::set<int>::const_iterator SIter;
1201 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1202 file<<
" "<< *neighbour;
1221 template<
typename M,
typename G,
typename L>
1223 const std::string& filename,
1225 bool readIndices=
true)
1227 using namespace MatrixMarketImpl;
1230 typedef typename LocalIndexT::Attribute Attribute;
1232 int rank = comm.communicator().rank();
1234 std::ostringstream rfilename;
1235 rfilename<<filename <<
"_"<<rank<<
".mm";
1237 file.open(rfilename.str().c_str(), std::ios::in);
1252 rfilename<<filename<<
"_"<<rank<<
".idx";
1253 file.open(rfilename.str().c_str());
1258 while(!file.eof() && file.peek()!=
'n') {
1267 pis.add(g,LocalIndexT(l,Attribute(c),b));
1275 if(s!=
"neighbours:")
1276 DUNE_THROW(MatrixMarketFormatError,
"was expecting the string: \"neighbours:\"");
1278 while(!file.eof()) {
1284 comm.ri.setNeighbours(nb);
1286 comm.ri.template rebuild<false>();
1301 template<
typename M>
1303 const std::string& filename)
1306 file.open(filename.c_str(), std::ios::in);
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:673
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:679
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1101
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1976
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1095
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:831
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:859
A vector of blocks with memory management.
Definition: bvector.hh:393
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:501
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:399
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Default exception class for I/O errors.
Definition: exceptions.hh:229
Index Set Interface base class.
Definition: indexidset.hh:76
auto size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:221
Exception indicating that the index set is not in the expected state.
Definition: indexset.hh:203
A generic dynamic dense matrix.
Definition: matrix.hh:559
Default exception for dummy implementations.
Definition: exceptions.hh:261
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:460
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:469
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
iterator begin()
Get an iterator over the indices positioned at the first index.
iterator end()
Get an iterator over the indices positioned after the last index.
LI LocalIndex
The type of the local index, e.g. ParallelLocalIndex.
Definition: indexset.hh:237
Stream & operator>>(Stream &stream, std::tuple< Ts... > &t)
Read a std::tuple.
Definition: streamoperators.hh:41
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:914
void storeMatrixMarket(const M &matrix, std::string filename, int prec=default_precision)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:1135
void loadMatrixMarket(M &matrix, const std::string &filename, OwnerOverlapCopyCommunication< G, L > &comm, bool readIndices=true)
Load a parallel matrix/vector stored in matrix market format.
Definition: matrixmarket.hh:1222
auto countNonZeros(const M &, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:117
EnableIfInterOperable< T1, T2, bool >::type operator<(const RandomAccessIteratorFacade< T1, V1, R1, D > &lhs, const RandomAccessIteratorFacade< T2, V2, R2, D > &rhs)
Comparison operator.
Definition: iteratorfacades.hh:635
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition: interface.hh:322
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:114
Some handy generic functions for ISTL matrices.
Dune namespace.
Definition: alignedallocator.hh:11
T lane(std::size_t l, const T &v)
access a lane of a simd vector (scalar version)
Definition: simd.hh:364
Classes providing communication interfaces for overlapping Schwarz methods.
Include file for users of the SIMD abstraction layer.
Standard Dune debug streams.
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Functor to the data values of the matrix.
Definition: matrixmarket.hh:675
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:682
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:700
a wrapper class of numeric values.
Definition: matrixmarket.hh:593
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:605
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:74
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:79