3 #ifndef DUNE_MATRIXMARKET_HH
4 #define DUNE_MATRIXMARKET_HH
15 #include <dune/common/fmatrix.hh>
16 #include <dune/common/tuples.hh>
17 #include <dune/common/unused.hh>
59 struct mm_numeric_type {
69 struct mm_numeric_type<int>
78 static std::string str()
85 struct mm_numeric_type<double>
94 static std::string str()
101 struct mm_numeric_type<float>
110 static std::string str()
117 struct mm_numeric_type<std::complex<double> >
126 static std::string str()
133 struct mm_numeric_type<std::complex<float> >
142 static std::string str()
157 struct mm_header_printer;
159 template<
typename T,
typename A,
int i,
int j>
160 struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,
A> >
162 static void print(std::ostream& os)
164 os<<
"%%MatrixMarket matrix coordinate ";
165 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
169 template<
typename B,
typename A>
170 struct mm_header_printer<BlockVector<B,
A> >
172 static void print(std::ostream& os)
174 os<<
"%%MatrixMarket matrix array ";
175 os<<mm_numeric_type<typename B::field_type>::str()<<
" general"<<std::endl;
179 template<
typename T,
int j>
180 struct mm_header_printer<FieldVector<T,j> >
182 static void print(std::ostream& os)
184 os<<
"%%MatrixMarket matrix array ";
185 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
189 template<
typename T,
int i,
int j>
190 struct mm_header_printer<FieldMatrix<T,i,j> >
192 static void print(std::ostream& os)
194 os<<
"%%MatrixMarket matrix array ";
195 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
208 struct mm_block_structure_header;
210 template<
typename T,
typename A,
int i>
211 struct mm_block_structure_header<BlockVector<FieldVector<T,i>,
A> >
213 typedef BlockVector<FieldVector<T,i>,
A> M;
215 static void print(std::ostream& os,
const M& m)
217 os<<
"% ISTL_STRUCT blocked ";
218 os<<i<<
" "<<1<<std::endl;
222 template<
typename T,
typename A,
int i,
int j>
223 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,
A> >
225 typedef BCRSMatrix<FieldMatrix<T,i,j>,
A> M;
227 static void print(std::ostream& os,
const M& m)
229 os<<
"% ISTL_STRUCT blocked ";
230 os<<i<<
" "<<j<<std::endl;
235 template<
typename T,
int i,
int j>
236 struct mm_block_structure_header<FieldMatrix<T,i,j> >
238 typedef FieldMatrix<T,i,j> M;
240 static void print(std::ostream& os,
const M& m)
244 template<
typename T,
int i>
245 struct mm_block_structure_header<FieldVector<T,i> >
247 typedef FieldVector<T,i> M;
249 static void print(std::ostream& os,
const M& m)
253 enum LineType { MM_HEADER, MM_ISTLSTRUCT, DATA };
254 enum { MM_MAX_LINE_LENGTH=1025 };
256 enum MM_TYPE { coordinate_type, array_type, unknown_type };
258 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
260 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
272 bool lineFeed(std::istream& file)
296 void skipComments(std::istream& file)
304 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
310 bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
319 std::cout<<buffer<<std::endl;
321 if(buffer!=
"%%MatrixMarket") {
323 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
334 if(buffer !=
"matrix")
337 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
351 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
358 if(buffer !=
"array")
360 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
363 mmHeader.type=array_type;
367 if(buffer !=
"coordinate")
369 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
372 mmHeader.type=coordinate_type;
375 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
389 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
395 if(buffer !=
"integer")
397 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
400 mmHeader.ctype=integer_type;
406 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
409 mmHeader.ctype=double_type;
413 if(buffer !=
"complex")
415 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
418 mmHeader.ctype=complex_type;
422 if(buffer !=
"pattern")
424 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
427 mmHeader.ctype=pattern;
430 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
439 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
445 if(buffer !=
"general")
447 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
450 mmHeader.structure=general;
454 if(buffer !=
"hermitian")
456 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
459 mmHeader.structure=hermitian;
462 if(buffer.size()==1) {
463 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
471 if(buffer !=
"symmetric")
473 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
476 mmHeader.structure=symmetric;
480 if(buffer !=
"skew-symmetric")
482 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
485 mmHeader.structure=skew_symmetric;
488 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
492 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
495 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
501 void readNextLine(std::istream& file, std::ostringstream&, LineType&
type)
503 DUNE_UNUSED_PARAMETER(type);
508 while(index==0&&!file.eof())
511 while(!file.eof() && (c=file.get())==
' ') ;
514 while(!file.eof() && (c=file.get())==
'\n') {
519 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
527 template<std::
size_t brows, std::
size_t bcols>
528 Dune::tuple<std::size_t, std::size_t, std::size_t>
529 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries,
const MMHeader& header)
531 std::size_t blockrows=rows/brows;
532 std::size_t blockcols=cols/bcols;
533 std::size_t blocksize=brows*bcols;
534 std::size_t blockentries=0;
536 switch(header.structure)
539 blockentries = entries/blocksize;
break;
540 case skew_symmetric :
541 blockentries = 2*entries/blocksize;
break;
543 blockentries = (2*entries-rows)/blocksize;
break;
545 blockentries = (2*entries-rows)/blocksize;
break;
547 throw Dune::NotImplemented();
549 return Dune::make_tuple(blockrows, blockcols, blockentries);
559 struct IndexData :
public T
576 struct NumericWrapper
592 struct NumericWrapper<PatternDummy>
596 std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
598 return is>>num.number;
601 std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
603 DUNE_UNUSED_PARAMETER(num);
613 bool operator<(const IndexData<T>& i1,
const IndexData<T>& i2)
615 return i1.index<i2.index;
624 std::istream& operator>>(std::istream& is, IndexData<T>& data)
629 return is>>data.number;
638 template<
typename D,
int brows,
int bcols>
639 struct MatrixValuesSetter
647 void operator()(
const std::vector<std::set<IndexData<D> > >& rows,
650 for(
typename M::RowIterator iter=matrix.begin();
651 iter!= matrix.end(); ++iter)
653 for(
typename M::size_type brow=iter.index()*brows,
654 browend=iter.index()*brows+brows;
655 brow<browend; ++brow)
657 typedef typename std::set<IndexData<D> >::const_iterator Siter;
658 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
659 siter != send; ++siter)
660 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
666 template<
int brows,
int bcols>
667 struct MatrixValuesSetter<PatternDummy,brows,bcols>
670 void operator()(
const std::vector<std::set<IndexData<PatternDummy> > >& rows,
675 template<
typename T,
typename A,
int brows,
int bcols,
typename D>
677 std::istream& file, std::size_t entries,
678 const MMHeader& mmHeader,
const D&)
684 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
686 for(; entries>0; --entries) {
692 assert(row/bcols<matrix.N());
694 assert(data.index/bcols<matrix.M());
695 rows[
row].insert(data);
699 if(mmHeader.structure!= general)
700 DUNE_THROW(Dune::NotImplemented,
"Only general is supported right now!");
704 for(
typename Matrix::CreateIterator iter=matrix.createbegin();
705 iter!= matrix.createend(); ++iter)
707 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
708 brow<browend; ++brow)
710 typedef typename std::set<IndexData<D> >::const_iterator Siter;
711 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
712 siter != send; ++siter, ++nnz)
713 iter.insert(siter->index/bcols);
720 MatrixValuesSetter<D,brows,bcols> Setter;
722 Setter(rows, matrix);
730 void mm_read_header(std::size_t& rows, std::size_t& cols, MMHeader& header, std::istream& istr,
733 if(!readMatrixMarketBanner(istr, header)) {
734 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
735 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
738 istr.seekg(0, std::ios::beg);
740 header.type=array_type;
755 template<
typename T,
typename A,
int entries>
760 for(
int i=0; size>0; ++i, --size) {
763 vector[i/entries][i%entries]=val;
774 template<
typename T,
typename A,
int entries>
779 std::size_t rows, cols;
784 if(header.type!=array_type)
787 std::size_t size=rows/entries;
788 if(size*entries!=rows)
793 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
805 template<
typename T,
typename A,
int brows,
int bcols>
810 if(!readMatrixMarketBanner(istr, header)) {
811 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
812 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
815 istr.seekg(0, std::ios::beg);
819 std::size_t rows, cols, entries;
835 std::size_t nnz, blockrows, blockcols;
837 Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
839 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
842 matrix.setSize(blockrows, blockcols);
845 if(header.type==array_type)
846 DUNE_THROW(Dune::NotImplemented,
"Array format currently not supported for matrices!");
848 readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
855 template<
typename B,
int i,
int j,
typename A>
864 template<
typename B,
int i,
int j>
872 for(riterator row=entry.begin(); row != entry.end(); ++
row, ++rowidx) {
874 for(citerator
col = row->begin();
col != row->end(); ++
col, ++coli)
875 ostr<< rowidx<<
" "<<coli<<
" "<<*
col<<std::endl;
882 const integral_constant<int,1>&)
884 ostr<<entry<<std::endl;
890 const integral_constant<int,0>&)
893 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
894 typedef typename V::const_iterator VIter;
896 for(VIter i=vector.begin(); i != vector.end(); ++i)
899 integral_constant<int,isnumeric>());
902 template<
typename T,
typename A,
int i>
905 return vector.size()*i;
911 const integral_constant<int,0>&)
914 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
922 const integral_constant<int,1>&)
928 typedef typename M::const_iterator riterator;
929 typedef typename M::ConstColIterator citerator;
930 for(riterator row=matrix.begin(); row != matrix.end(); ++
row)
931 for(citerator
col = row->begin();
col != row->end(); ++
col)
946 mm_header_printer<M>::print(ostr);
947 mm_block_structure_header<M>::print(ostr,matrix);
965 std::string filename)
967 std::ofstream file(filename.c_str());
968 file.setf(std::ios::scientific,std::ios::floatfield);
987 template<
typename M,
typename G,
typename L>
989 std::string filename,
991 bool storeIndices=
true)
996 std::ostringstream rfilename;
997 rfilename<<filename <<
"_"<<rank<<
".mm";
998 std::cout<< rfilename.str()<<std::endl;
999 std::ofstream file(rfilename.str().c_str());
1000 file.setf(std::ios::scientific,std::ios::floatfield);
1009 rfilename<<filename<<
"_"<<rank<<
".idx";
1010 file.open(rfilename.str().c_str());
1011 file.setf(std::ios::scientific,std::ios::floatfield);
1013 typedef typename IndexSet::const_iterator Iterator;
1014 for(Iterator iter = comm.
indexSet().begin();
1015 iter != comm.
indexSet().end(); ++iter) {
1016 file << iter->global()<<
" "<<(std::size_t)iter->local()<<
" "
1017 <<(int)iter->local().attribute()<<
" "<<(int)iter->local().isPublic()<<std::endl;
1020 file<<
"neighbours:";
1021 const std::set<int>& neighbours=comm.
remoteIndices().getNeighbours();
1022 typedef std::set<int>::const_iterator SIter;
1023 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1024 file<<
" "<< *neighbour;
1043 template<
typename M,
typename G,
typename L>
1045 const std::string& filename,
1047 bool readIndices=
true)
1050 typedef typename LocalIndex::Attribute Attribute;
1054 std::ostringstream rfilename;
1055 rfilename<<filename <<
"_"<<rank<<
".mm";
1057 file.open(rfilename.str().c_str(), std::ios::in);
1059 DUNE_THROW(IOError,
"Could not open file" << rfilename.str().c_str());
1070 IndexSet& pis=comm.pis;
1072 rfilename<<filename<<
"_"<<rank<<
".idx";
1073 file.open(rfilename.str().c_str());
1075 DUNE_THROW(InvalidIndexSetState,
"Index set is not empty!");
1078 while(!file.eof() && file.peek()!=
'n') {
1087 pis.add(g,LocalIndex(l,Attribute(c),b));
1095 if(s!=
"neighbours:")
1098 while(!file.eof()) {
1104 comm.ri.setNeighbours(nb);
1106 comm.ri.template rebuild<false>();
1121 template<
typename M>
1123 const std::string& filename)
1126 file.open(filename.c_str(), std::ios::in);
1128 DUNE_THROW(IOError,
"Could not open file" << filename);
int countNonZeros(const M &matrix)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:155
Classes providing communication interfaces for overlapping Schwarz methods.
std::size_t countEntries(const BlockVector< FieldVector< T, i >, A > &vector)
Definition: matrixmarket.hh:903
void mm_read_header(std::size_t &rows, std::size_t &cols, MMHeader &header, std::istream &istr, bool isVector)
Definition: matrixmarket.hh:730
void storeMatrixMarket(const M &matrix, std::string filename)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:964
MM_STRUCTURE structure
Definition: matrixmarket.hh:269
void mm_print_vector_entry(const V &entry, std::ostream &ostr, const integral_constant< int, 1 > &)
Definition: matrixmarket.hh:881
Some handy generic functions for ISTL matrices.
A vector of blocks with memory management.
Definition: bvector.hh:253
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
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, std::ostream &ostr)
Definition: matrixmarket.hh:865
Col col
Definition: matrixmatrix.hh:347
void readMatrixMarket(Dune::BlockVector< Dune::FieldVector< T, entries >, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:775
const CollectiveCommunication< MPI_Comm > & communicator() const
Definition: owneroverlapcopy.hh:306
Matrix & A
Definition: matrixmatrix.hh:216
Implementation of the BCRSMatrix class.
MM_CTYPE ctype
Definition: matrixmarket.hh:268
MM_TYPE type
Definition: matrixmarket.hh:267
void writeMatrixMarket(const V &vector, std::ostream &ostr, const integral_constant< int, 0 > &)
Definition: matrixmarket.hh:910
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:452
void mm_read_vector_entries(Dune::BlockVector< Dune::FieldVector< T, entries >, A > &vector, std::size_t size, std::istream &istr)
Definition: matrixmarket.hh:756
Row row
Definition: matrixmatrix.hh:345
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:453
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:475
std::size_t index
Definition: matrixmarket.hh:561
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:466
T number
Definition: matrixmarket.hh:578
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:1044
Definition: matrixmarket.hh:852
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
Definition: matrixmarket.hh:726
Definition: matrixutils.hh:25