1 #ifndef DUNE_MATRIXMARKET_HH
2 #define DUNE_MATRIXMARKET_HH
13 #include<dune/common/fmatrix.hh>
14 #include<dune/common/tuples.hh>
15 #include<dune/common/misc.hh>
57 struct mm_numeric_type{
67 struct mm_numeric_type<int>
76 static std::string str()
83 struct mm_numeric_type<double>
92 static std::string str()
99 struct mm_numeric_type<float>
108 static std::string str()
115 struct mm_numeric_type<std::complex<double> >
124 static std::string str()
131 struct mm_numeric_type<std::complex<float> >
140 static std::string str()
155 struct mm_header_printer;
157 template<
typename T,
typename A,
int i,
int j>
158 struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,
A> >
160 static void print(std::ostream& os)
162 os<<
"%%MatrixMarket matrix coordinate ";
163 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
167 template<
typename B,
typename A>
168 struct mm_header_printer<BlockVector<B,
A> >
170 static void print(std::ostream& os)
172 os<<
"%%MatrixMarket matrix array ";
173 os<<mm_numeric_type<typename B::field_type>::str()<<
" general"<<std::endl;
177 template<
typename T,
int j>
178 struct mm_header_printer<FieldVector<T,j> >
180 static void print(std::ostream& os)
182 os<<
"%%MatrixMarket matrix array ";
183 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
187 template<
typename T,
int i,
int j>
188 struct mm_header_printer<FieldMatrix<T,i,j> >
190 static void print(std::ostream& os)
192 os<<
"%%MatrixMarket matrix array ";
193 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
206 struct mm_block_structure_header;
208 template<
typename T,
typename A,
int i>
209 struct mm_block_structure_header<BlockVector<FieldVector<T,i>,
A> >
211 typedef BlockVector<FieldVector<T,i>,
A> M;
213 static void print(std::ostream& os,
const M& m)
215 os<<
"% ISTL_STRUCT blocked ";
216 os<<i<<
" "<<1<<std::endl;
220 template<
typename T,
typename A,
int i,
int j>
221 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,
A> >
223 typedef BCRSMatrix<FieldMatrix<T,i,j>,
A> M;
225 static void print(std::ostream& os,
const M& m)
227 os<<
"% ISTL_STRUCT blocked ";
228 os<<i<<
" "<<j<<std::endl;
233 template<
typename T,
int i,
int j>
234 struct mm_block_structure_header<FieldMatrix<T,i,j> >
236 typedef FieldMatrix<T,i,j> M;
238 static void print(std::ostream& os,
const M& m)
242 template<
typename T,
int i>
243 struct mm_block_structure_header<FieldVector<T,i> >
245 typedef FieldVector<T,i> M;
247 static void print(std::ostream& os,
const M& m)
251 enum LineType{ MM_HEADER, MM_ISTLSTRUCT, DATA };
252 enum{ MM_MAX_LINE_LENGTH=1025 };
254 enum MM_TYPE { coordinate_type, array_type, unknown_type };
256 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
258 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
270 bool lineFeed(std::istream& file)
294 void skipComments(std::istream& file)
302 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
308 bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
317 std::cout<<buffer<<std::endl;
319 if(buffer!=
"%%MatrixMarket"){
321 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
332 if(buffer !=
"matrix")
335 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
349 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
356 if(buffer !=
"array")
358 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
361 mmHeader.type=array_type;
365 if(buffer !=
"coordinate")
367 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
370 mmHeader.type=coordinate_type;
373 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
387 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
393 if(buffer !=
"integer")
395 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
398 mmHeader.ctype=integer_type;
404 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
407 mmHeader.ctype=double_type;
411 if(buffer !=
"complex")
413 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
416 mmHeader.ctype=complex_type;
420 if(buffer !=
"pattern")
422 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
425 mmHeader.ctype=pattern;
428 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
437 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
443 if(buffer !=
"general")
445 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
448 mmHeader.structure=general;
452 if(buffer !=
"hermitian")
454 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
457 mmHeader.structure=hermitian;
460 if(buffer.size()==1){
461 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
469 if(buffer !=
"symmetric")
471 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
474 mmHeader.structure=symmetric;
478 if(buffer !=
"skew-symmetric")
480 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
483 mmHeader.structure=skew_symmetric;
486 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
490 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
493 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
499 void readNextLine(std::istream& file, std::ostringstream& line, LineType&
type)
505 while(index==0&&!file.eof())
508 while(!file.eof() && (c=file.get())==
' ');
511 while(!file.eof() && (c=file.get())==
'\n'){
516 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
524 template<std::
size_t brows, std::
size_t bcols>
525 Dune::tuple<std::size_t, std::size_t, std::size_t>
526 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries,
const MMHeader& header)
528 std::size_t blockrows=rows/brows;
529 std::size_t blockcols=cols/bcols;
530 std::size_t blocksize=brows*bcols;
531 std::size_t blockentries=0;
533 switch(header.structure)
536 blockentries = entries/blocksize;
break;
538 blockentries = 2*entries/blocksize;
break;
540 blockentries = (2*entries-rows)/blocksize;
break;
542 blockentries = (2*entries-rows)/blocksize;
break;
544 throw Dune::NotImplemented();
546 return Dune::make_tuple(blockrows, blockcols, blockentries);
556 struct IndexData :
public T
573 struct NumericWrapper
589 struct NumericWrapper<PatternDummy>
593 std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
595 return is>>num.number;
598 std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
609 bool operator<(const IndexData<T>& i1,
const IndexData<T>& i2)
611 return i1.index<i2.index;
620 std::istream& operator>>(std::istream& is, IndexData<T>& data)
625 return is>>data.number;
634 template<
typename D,
int brows,
int bcols>
635 struct MatrixValuesSetter
643 void operator()(
const std::vector<std::set<IndexData<D> > >& rows,
646 for(
typename M::RowIterator iter=matrix.begin();
647 iter!= matrix.end(); ++iter)
649 for(
typename M::size_type brow=iter.index()*brows,
650 browend=iter.index()*brows+brows;
651 brow<browend; ++brow)
653 typedef typename std::set<IndexData<D> >::const_iterator Siter;
654 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
655 siter != send; ++siter)
656 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
662 template<
int brows,
int bcols>
663 struct MatrixValuesSetter<PatternDummy,brows,bcols>
666 void operator()(
const std::vector<std::set<IndexData<PatternDummy> > >& rows,
671 template<
typename T,
typename A,
int brows,
int bcols,
typename D>
672 void readSparseEntries(
Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,
A>& matrix,
673 std::istream& file, std::size_t entries,
674 const MMHeader& mmHeader,
const D&)
680 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
682 for(; entries>0;--entries){
688 assert(row/bcols<matrix.N());
690 assert(data.index/bcols<matrix.M());
691 rows[
row].insert(data);
695 if(mmHeader.structure!= general)
696 DUNE_THROW(Dune::NotImplemented,
"Only general is supported right now!");
700 for(
typename Matrix::CreateIterator iter=matrix.createbegin();
701 iter!= matrix.createend(); ++iter)
703 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
704 brow<browend; ++brow)
706 typedef typename std::set<IndexData<D> >::const_iterator Siter;
707 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
708 siter != send; ++siter, ++nnz)
709 iter.insert(siter->index/bcols);
716 MatrixValuesSetter<D,brows,bcols> Setter;
718 Setter(rows, matrix);
726 void mm_read_header(std::size_t& rows, std::size_t& cols, MMHeader& header, std::istream& istr,
729 if(!readMatrixMarketBanner(istr, header)){
730 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
731 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
734 istr.seekg(0, std::ios::beg);
736 header.type=array_type;
751 template<
typename T,
typename A,
int entries>
756 for(
int i=0;size>0; ++i, --size){
759 vector[i/entries][i%entries]=val;
770 template<
typename T,
typename A,
int entries>
775 std::size_t rows, cols;
780 if(header.type!=array_type)
783 std::size_t size=rows/entries;
784 if(size*entries!=rows)
789 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
801 template<
typename T,
typename A,
int brows,
int bcols>
809 if(!readMatrixMarketBanner(istr, header)){
810 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
811 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
814 istr.seekg(0, std::ios::beg);
818 std::size_t rows, cols, entries;
834 std::size_t nnz, blockrows, blockcols;
836 Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
838 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
841 matrix.setSize(blockrows, blockcols);
842 matrix.setBuildMode(
Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,
A>::row_wise);
844 if(header.type==array_type)
845 DUNE_THROW(Dune::NotImplemented,
"Array format currently not supported for matrices!");
847 readSparseEntries(matrix, istr, entries, header, NumericWrapper<double>());
854 template<
typename B,
int i,
int j,
typename A>
863 template<
typename B,
int i,
int j>
865 typename FieldMatrix<B,i,j>::size_type rowidx,
866 typename FieldMatrix<B,i,j>::size_type colidx,
869 typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
870 typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
871 for(riterator row=entry.begin(); row != entry.end(); ++
row, ++rowidx){
873 for(citerator
col = row->begin();
col != row->end(); ++
col, ++coli)
874 ostr<< rowidx<<
" "<<coli<<
" "<<*
col<<std::endl;
881 const integral_constant<int,1>&)
883 ostr<<entry<<std::endl;
889 const integral_constant<int,0>&)
892 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
893 typedef typename V::const_iterator VIter;
895 for(VIter i=vector.begin(); i != vector.end(); ++i)
898 integral_constant<int,isnumeric>());
901 template<
typename T,
typename A,
int i>
904 return vector.size()*i;
910 const integral_constant<int,0>&)
913 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
921 const integral_constant<int,1>&)
927 typedef typename M::const_iterator riterator;
928 typedef typename M::ConstColIterator citerator;
929 for(riterator row=matrix.begin(); row != matrix.end(); ++
row)
930 for(citerator
col = row->begin();
col != row->end(); ++
col)
945 mm_header_printer<M>::print(ostr);
946 mm_block_structure_header<M>::print(ostr,matrix);
964 std::string filename)
966 std::ofstream file(filename.c_str());
967 file.setf(std::ios::scientific,std::ios::floatfield);
986 template<
typename M,
typename G,
typename L>
988 std::string filename,
990 bool storeIndices=
true)
995 std::ostringstream rfilename;
996 rfilename<<filename <<
"_"<<rank<<
".mm";
997 std::cout<< rfilename.str()<<std::endl;
998 std::ofstream file(rfilename.str().c_str());
999 file.setf(std::ios::scientific,std::ios::floatfield);
1008 rfilename<<filename<<
"_"<<rank<<
".idx";
1009 file.open(rfilename.str().c_str());
1010 file.setf(std::ios::scientific,std::ios::floatfield);
1012 typedef typename IndexSet::const_iterator Iterator;
1013 for(Iterator iter = comm.
indexSet().begin();
1014 iter != comm.
indexSet().end(); ++iter){
1015 file << iter->global()<<
" "<<(std::size_t)iter->local()<<
" "
1016 <<(int)iter->local().attribute()<<
" "<<(int)iter->local().isPublic()<<std::endl;
1019 file<<
"neighbours:";
1020 const std::set<int>& neighbours=comm.
remoteIndices().getNeighbours();
1021 typedef std::set<int>::const_iterator SIter;
1022 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour){
1023 file<<
" "<< *neighbour;
1042 template<
typename M,
typename G,
typename L>
1044 const std::string& filename,
1046 bool readIndices=
true)
1049 typedef typename LocalIndex::Attribute Attribute;
1053 std::ostringstream rfilename;
1054 rfilename<<filename <<
"_"<<rank<<
".mm";
1056 file.open(rfilename.str().c_str(), std::ios::in);
1058 DUNE_THROW(IOError,
"Could not open file" << rfilename.str().c_str());
1069 IndexSet& pis=comm.pis;
1071 rfilename<<filename<<
"_"<<rank<<
".idx";
1072 file.open(rfilename.str().c_str());
1074 DUNE_THROW(InvalidIndexSetState,
"Index set is not empty!");
1077 while(!file.eof() && file.peek()!=
'n'){
1086 pis.add(g,LocalIndex(l,Attribute(c),b));
1094 if(s!=
"neighbours:")
1103 comm.ri.setNeighbours(nb);
1105 comm.ri.template rebuild<false>();
1120 template<
typename M>
1122 const std::string& filename)
1125 file.open(filename.c_str(), std::ios::in);
1127 DUNE_THROW(IOError,
"Could not open file" << filename);