3#ifndef DUNE_ISTL_MATRIXMARKET_HH
4#define DUNE_ISTL_MATRIXMARKET_HH
47 namespace MatrixMarketImpl
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()
159 template<
typename T,
typename A,
int i,
int j>
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>
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;
210 template<
typename T,
typename A,
int i>
215 static void print(std::ostream& os,
const M&)
217 os<<
"% ISTL_STRUCT blocked ";
218 os<<i<<
" "<<1<<std::endl;
222 template<
typename T,
typename A,
int i,
int j>
227 static void print(std::ostream& os,
const 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 };
265 : type(coordinate_type), ctype(double_type), structure(general)
269 MM_STRUCTURE structure;
272 inline bool lineFeed(std::istream& file)
296 inline void skipComments(std::istream& file)
304 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
310 inline 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 template<std::
size_t brows, std::
size_t bcols>
502 Dune::tuple<std::size_t, std::size_t, std::size_t>
503 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries,
const MMHeader& header)
505 std::size_t blockrows=rows/brows;
506 std::size_t blockcols=cols/bcols;
507 std::size_t blocksize=brows*bcols;
508 std::size_t blockentries=0;
510 switch(header.structure)
513 blockentries = entries/blocksize;
break;
514 case skew_symmetric :
515 blockentries = 2*entries/blocksize;
break;
517 blockentries = (2*entries-rows)/blocksize;
break;
519 blockentries = (2*entries-rows)/blocksize;
break;
523 return Dune::make_tuple(blockrows, blockcols, blockentries);
533 struct IndexData :
public T
570 std::istream&
operator>>(std::istream& is, NumericWrapper<T>& num)
572 return is>>num.number;
575 inline std::istream&
operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
587 bool operator<(
const IndexData<T>& i1,
const IndexData<T>& i2)
589 return i1.index<i2.index;
598 std::istream&
operator>>(std::istream& is, IndexData<T>& data)
603 return is>>data.number;
612 template<
typename D,
int brows,
int bcols>
621 void operator()(
const std::vector<std::set<IndexData<D> > >& rows,
624 for(
typename M::RowIterator iter=matrix.begin();
625 iter!= matrix.end(); ++iter)
627 for(
typename M::size_type brow=iter.index()*brows,
628 browend=iter.index()*brows+brows;
629 brow<browend; ++brow)
631 typedef typename std::set<IndexData<D> >::const_iterator Siter;
632 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
633 siter != send; ++siter)
634 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
640 template<
int brows,
int bcols>
641 struct MatrixValuesSetter<PatternDummy,brows,bcols>
644 void operator()(
const std::vector<std::set<IndexData<PatternDummy> > >& rows,
649 template<
typename T,
typename A,
int brows,
int bcols,
typename D>
651 std::istream& file, std::size_t entries,
652 const MMHeader& mmHeader,
const D&)
658 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
660 for(; entries>0; --entries) {
666 assert(row/bcols<matrix.N());
668 assert(data.index/bcols<matrix.M());
669 rows[row].insert(data);
673 if(mmHeader.structure!= general)
678 for(
typename Matrix::CreateIterator iter=matrix.createbegin();
679 iter!= matrix.createend(); ++iter)
681 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
682 brow<browend; ++brow)
684 typedef typename std::set<IndexData<D> >::const_iterator Siter;
685 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
686 siter != send; ++siter, ++nnz)
687 iter.insert(siter->index/bcols);
694 MatrixValuesSetter<D,brows,bcols> Setter;
696 Setter(rows, matrix);
704 inline void mm_read_header(std::size_t& rows, std::size_t& cols,
705 MatrixMarketImpl::MMHeader& header, std::istream& istr,
708 using namespace MatrixMarketImpl;
710 if(!readMatrixMarketBanner(istr, header)) {
711 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
712 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
715 istr.seekg(0, std::ios::beg);
717 header.type=array_type;
723 throw MatrixMarketFormatError();
728 throw MatrixMarketFormatError();
732 template<
typename T,
typename A,
int entries>
737 for(
int i=0; size>0; ++i, --size) {
740 vector[i/entries][i%entries]=val;
751 template<
typename T,
typename A,
int entries>
755 using namespace MatrixMarketImpl;
758 std::size_t rows, cols;
759 mm_read_header(rows,cols,header,istr,
true);
761 DUNE_THROW(MatrixMarketFormatError,
"cols!=1, therefore this is no vector!");
763 if(header.type!=array_type)
764 DUNE_THROW(MatrixMarketFormatError,
"Vectors have to be stored in array format!");
766 std::size_t size=rows/entries;
767 if(size*entries!=rows)
768 DUNE_THROW(MatrixMarketFormatError,
"Block size of vector is not correct1");
772 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
774 mm_read_vector_entries(vector, rows, istr);
784 template<
typename T,
typename A,
int brows,
int bcols>
788 using namespace MatrixMarketImpl;
791 if(!readMatrixMarketBanner(istr, header)) {
792 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
793 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
796 istr.seekg(0, std::ios::beg);
800 std::size_t rows, cols, entries;
803 throw MatrixMarketFormatError();
808 throw MatrixMarketFormatError();
812 throw MatrixMarketFormatError();
816 std::size_t nnz, blockrows, blockcols;
818 Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
820 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
823 matrix.setSize(blockrows, blockcols);
826 if(header.type==array_type)
829 readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
833 struct mm_multipliers
836 template<
typename B,
int i,
int j,
typename A>
837 struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
845 template<
typename B,
int i,
int j>
846 void mm_print_entry(
const FieldMatrix<B,i,j>& entry,
847 typename FieldMatrix<B,i,j>::size_type rowidx,
848 typename FieldMatrix<B,i,j>::size_type colidx,
853 for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
855 for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
856 ostr<< rowidx<<
" "<<coli<<
" "<<*col<<std::endl;
862 void mm_print_vector_entry(
const V& entry, std::ostream& ostr,
863 const integral_constant<int,1>&)
865 ostr<<entry<<std::endl;
870 void mm_print_vector_entry(
const V& vector, std::ostream& ostr,
871 const integral_constant<int,0>&)
873 using namespace MatrixMarketImpl;
876 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
877 typedef typename V::const_iterator VIter;
879 for(VIter i=vector.begin(); i != vector.end(); ++i)
881 mm_print_vector_entry(*i, ostr,
882 integral_constant<int,isnumeric>());
885 template<
typename T,
typename A,
int i>
886 std::size_t countEntries(
const BlockVector<FieldVector<T,i>,A>& vector)
888 return vector.size()*i;
893 void writeMatrixMarket(
const V& vector, std::ostream& ostr,
894 const integral_constant<int,0>&)
896 using namespace MatrixMarketImpl;
898 ostr<<countEntries(vector)<<
" "<<1<<std::endl;
899 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
900 mm_print_vector_entry(vector,ostr, integral_constant<int,isnumeric>());
905 void writeMatrixMarket(
const M& matrix,
907 const integral_constant<int,1>&)
909 ostr<<matrix.N()*mm_multipliers<M>::rows<<
" "
910 <<matrix.M()*mm_multipliers<M>::cols<<
" "
913 typedef typename M::const_iterator riterator;
914 typedef typename M::ConstColIterator citerator;
915 for(riterator row=matrix.begin(); row != matrix.end(); ++row)
916 for(citerator col = row->begin(); col != row->end(); ++col)
918 mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
919 col.index()*mm_multipliers<M>::cols+1, ostr);
927 void writeMatrixMarket(
const M& matrix,
930 using namespace MatrixMarketImpl;
933 mm_header_printer<M>::print(ostr);
934 mm_block_structure_header<M>::print(ostr,matrix);
952 std::string filename)
954 std::ofstream file(filename.c_str());
955 file.setf(std::ios::scientific,std::ios::floatfield);
956 writeMatrixMarket(matrix, file);
974 template<
typename M,
typename G,
typename L>
976 std::string filename,
978 bool storeIndices=
true)
981 int rank = comm.communicator().
rank();
983 std::ostringstream rfilename;
984 rfilename<<filename <<
"_"<<rank<<
".mm";
985 std::cout<< rfilename.str()<<std::endl;
986 std::ofstream file(rfilename.str().c_str());
987 file.setf(std::ios::scientific,std::ios::floatfield);
988 writeMatrixMarket(matrix, file);
996 rfilename<<filename<<
"_"<<rank<<
".idx";
997 file.open(rfilename.str().c_str());
998 file.setf(std::ios::scientific,std::ios::floatfield);
1000 typedef typename IndexSet::const_iterator Iterator;
1003 file << iter->global()<<
" "<<(std::size_t)iter->local()<<
" "
1004 <<(int)iter->local().attribute()<<
" "<<(int)iter->local().isPublic()<<std::endl;
1007 file<<
"neighbours:";
1008 const std::set<int>& neighbours=comm.
remoteIndices().getNeighbours();
1009 typedef std::set<int>::const_iterator SIter;
1010 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1011 file<<
" "<< *neighbour;
1030 template<
typename M,
typename G,
typename L>
1032 const std::string& filename,
1034 bool readIndices=
true)
1036 using namespace MatrixMarketImpl;
1039 typedef typename LocalIndex::Attribute Attribute;
1041 int rank = comm.communicator().
rank();
1043 std::ostringstream rfilename;
1044 rfilename<<filename <<
"_"<<rank<<
".mm";
1046 file.open(rfilename.str().c_str(), std::ios::in);
1061 rfilename<<filename<<
"_"<<rank<<
".idx";
1062 file.open(rfilename.str().c_str());
1067 while(!file.eof() && file.peek()!=
'n') {
1084 if(s!=
"neighbours:")
1085 DUNE_THROW(MatrixMarketFormatError,
"was exspecting the string: \"neighbours:\"");
1087 while(!file.eof()) {
1093 comm.ri.setNeighbours(nb);
1095 comm.ri.template rebuild<false>();
1110 template<
typename M>
1112 const std::string& filename)
1115 file.open(filename.c_str(), std::ios::in);
Implementation of the BCRSMatrix class.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
A vector of blocks with memory management.
Definition: bvector.hh:253
int rank() const
Return rank, is between 0 and size()-1.
Definition: mpicollectivecommunication.hh:166
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:283
remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:287
Base class for Dune-Exceptions.
Definition: exceptions.hh:91
A dense n x m matrix.
Definition: fmatrix.hh:67
vector space out of a tensor product of fields.
Definition: fvector.hh:94
Default exception class for I/O errors.
Definition: exceptions.hh:256
Index Set Interface base class.
Definition: indexidset.hh:76
Exception indicating that the index set is not in the expected state.
Definition: indexset.hh:204
An index present on the local process.
Definition: localindex.hh:33
Default exception for dummy implementations.
Definition: exceptions.hh:288
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:173
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:464
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:473
Implements a matrix constructed from a given type representing a field and compile-time given number ...
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:238
std::istream & operator>>(std::istream &is, tuple< T1 > &t)
Read a tuple.
Definition: tuples.hh:203
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
void readMatrixMarket(Dune::BlockVector< Dune::FieldVector< T, entries >, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:752
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:1031
void storeMatrixMarket(const M &matrix, std::string filename)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:951
int countNonZeros(const M &matrix)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:155
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:626
Some handy generic functions for ISTL matrices.
Dune namespace.
Definition: alignment.hh:10
Classes providing communication interfaces for overlapping Schwarz methods.
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:480
Functor to the data values of the matrix.
Definition: matrixmarket.hh:614
void operator()(const std::vector< std::set< IndexData< D > > > &rows, M &matrix)
Sets the matrixvalues.
Definition: matrixmarket.hh:621
a wrapper class of numeric values.
Definition: matrixmarket.hh:551
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:563
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:59
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:64
Fallback implementation of the std::tuple class.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18