dune-istl
2.2.0
|
00001 #ifndef DUNE_MATRIXMARKET_HH 00002 #define DUNE_MATRIXMARKET_HH 00003 00004 #include<ostream> 00005 #include<istream> 00006 #include<fstream> 00007 #include<sstream> 00008 #include<limits> 00009 #include<ios> 00010 #include"matrixutils.hh" 00011 #include "bcrsmatrix.hh" 00012 #include"owneroverlapcopy.hh" 00013 #include<dune/common/fmatrix.hh> 00014 #include<dune/common/tuples.hh> 00015 #include<dune/common/misc.hh> 00016 00017 namespace Dune 00018 { 00019 00045 namespace 00046 { 00056 template<class T> 00057 struct mm_numeric_type{ 00058 enum{ 00062 is_numeric=false 00063 }; 00064 }; 00065 00066 template<> 00067 struct mm_numeric_type<int> 00068 { 00069 enum{ 00073 is_numeric=true 00074 }; 00075 00076 static std::string str() 00077 { 00078 return "integer"; 00079 } 00080 }; 00081 00082 template<> 00083 struct mm_numeric_type<double> 00084 { 00085 enum{ 00089 is_numeric=true 00090 }; 00091 00092 static std::string str() 00093 { 00094 return "real"; 00095 } 00096 }; 00097 00098 template<> 00099 struct mm_numeric_type<float> 00100 { 00101 enum{ 00105 is_numeric=true 00106 }; 00107 00108 static std::string str() 00109 { 00110 return "real"; 00111 } 00112 }; 00113 00114 template<> 00115 struct mm_numeric_type<std::complex<double> > 00116 { 00117 enum{ 00121 is_numeric=true 00122 }; 00123 00124 static std::string str() 00125 { 00126 return "complex"; 00127 } 00128 }; 00129 00130 template<> 00131 struct mm_numeric_type<std::complex<float> > 00132 { 00133 enum{ 00137 is_numeric=true 00138 }; 00139 00140 static std::string str() 00141 { 00142 return "complex"; 00143 } 00144 }; 00145 00154 template<class M> 00155 struct mm_header_printer; 00156 00157 template<typename T, typename A, int i, int j> 00158 struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> > 00159 { 00160 static void print(std::ostream& os) 00161 { 00162 os<<"%%MatrixMarket matrix coordinate "; 00163 os<<mm_numeric_type<T>::str()<<" general"<<std::endl; 00164 } 00165 }; 00166 00167 template<typename B, typename A> 00168 struct mm_header_printer<BlockVector<B,A> > 00169 { 00170 static void print(std::ostream& os) 00171 { 00172 os<<"%%MatrixMarket matrix array "; 00173 os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl; 00174 } 00175 }; 00176 00177 template<typename T, int j> 00178 struct mm_header_printer<FieldVector<T,j> > 00179 { 00180 static void print(std::ostream& os) 00181 { 00182 os<<"%%MatrixMarket matrix array "; 00183 os<<mm_numeric_type<T>::str()<<" general"<<std::endl; 00184 } 00185 }; 00186 00187 template<typename T, int i, int j> 00188 struct mm_header_printer<FieldMatrix<T,i,j> > 00189 { 00190 static void print(std::ostream& os) 00191 { 00192 os<<"%%MatrixMarket matrix array "; 00193 os<<mm_numeric_type<T>::str()<<" general"<<std::endl; 00194 } 00195 }; 00196 00205 template<class M> 00206 struct mm_block_structure_header; 00207 00208 template<typename T, typename A, int i> 00209 struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> > 00210 { 00211 typedef BlockVector<FieldVector<T,i>,A> M; 00212 00213 static void print(std::ostream& os, const M& m) 00214 { 00215 os<<"% ISTL_STRUCT blocked "; 00216 os<<i<<" "<<1<<std::endl; 00217 } 00218 }; 00219 00220 template<typename T, typename A, int i, int j> 00221 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> > 00222 { 00223 typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M; 00224 00225 static void print(std::ostream& os, const M& m) 00226 { 00227 os<<"% ISTL_STRUCT blocked "; 00228 os<<i<<" "<<j<<std::endl; 00229 } 00230 }; 00231 00232 00233 template<typename T, int i, int j> 00234 struct mm_block_structure_header<FieldMatrix<T,i,j> > 00235 { 00236 typedef FieldMatrix<T,i,j> M; 00237 00238 static void print(std::ostream& os, const M& m) 00239 {} 00240 }; 00241 00242 template<typename T, int i> 00243 struct mm_block_structure_header<FieldVector<T,i> > 00244 { 00245 typedef FieldVector<T,i> M; 00246 00247 static void print(std::ostream& os, const M& m) 00248 {} 00249 }; 00250 00251 enum LineType{ MM_HEADER, MM_ISTLSTRUCT, DATA }; 00252 enum{ MM_MAX_LINE_LENGTH=1025 }; 00253 00254 enum MM_TYPE { coordinate_type, array_type, unknown_type }; 00255 00256 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype }; 00257 00258 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure }; 00259 00260 struct MMHeader 00261 { 00262 MMHeader() 00263 : type(coordinate_type), ctype(double_type), structure(general) 00264 {} 00265 MM_TYPE type; 00266 MM_CTYPE ctype; 00267 MM_STRUCTURE structure; 00268 }; 00269 00270 bool lineFeed(std::istream& file) 00271 { 00272 char c; 00273 if(!file.eof()) 00274 c=file.peek(); 00275 else 00276 return false; 00277 // ignore whitespace 00278 while(c==' ') 00279 { 00280 file.get(); 00281 if(file.eof()) 00282 return false; 00283 c=file.peek(); 00284 } 00285 00286 if(c=='\n'){ 00287 /* eat the line feed */ 00288 file.get(); 00289 return true; 00290 } 00291 return false; 00292 } 00293 00294 void skipComments(std::istream& file) 00295 { 00296 lineFeed(file); 00297 char c=file.peek(); 00298 // ignore comment lines 00299 while(c=='%') 00300 { 00301 /* disgard the rest of the line */ 00302 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00303 c=file.peek(); 00304 } 00305 } 00306 00307 00308 bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader) 00309 { 00310 std::string buffer; 00311 char c; 00312 file >> buffer; 00313 c=buffer[0]; 00314 mmHeader=MMHeader(); 00315 if(c!='%') 00316 return false; 00317 std::cout<<buffer<<std::endl; 00318 /* read the banner */ 00319 if(buffer!="%%MatrixMarket"){ 00320 /* disgard the rest of the line */ 00321 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00322 return false; 00323 } 00324 00325 if(lineFeed(file)) 00326 /* premature end of line */ 00327 return false; 00328 00329 /* read the matrix_type */ 00330 file >> buffer; 00331 00332 if(buffer != "matrix") 00333 { 00334 /* disgard the rest of the line */ 00335 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00336 return false; 00337 } 00338 00339 if(lineFeed(file)) 00340 /* premature end of line */ 00341 return false; 00342 00343 /* The type of the matrix */ 00344 file >> buffer; 00345 00346 if(buffer.empty()) 00347 return false; 00348 00349 std::transform(buffer.begin(), buffer.end(), buffer.begin(), 00350 tolower); 00351 00352 switch(buffer[0]) 00353 { 00354 case 'a': 00355 /* sanity check */ 00356 if(buffer != "array") 00357 { 00358 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00359 return false; 00360 } 00361 mmHeader.type=array_type; 00362 break; 00363 case 'c': 00364 /* sanity check */ 00365 if(buffer != "coordinate") 00366 { 00367 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00368 return false; 00369 } 00370 mmHeader.type=coordinate_type; 00371 break; 00372 default: 00373 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00374 return false; 00375 } 00376 00377 if(lineFeed(file)) 00378 /* premature end of line */ 00379 return false; 00380 00381 /* The numeric type used. */ 00382 file >> buffer; 00383 00384 if(buffer.empty()) 00385 return false; 00386 00387 std::transform(buffer.begin(), buffer.end(), buffer.begin(), 00388 tolower); 00389 switch(buffer[0]) 00390 { 00391 case 'i': 00392 /* sanity check */ 00393 if(buffer != "integer") 00394 { 00395 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00396 return false; 00397 } 00398 mmHeader.ctype=integer_type; 00399 break; 00400 case 'r': 00401 /* sanity check */ 00402 if(buffer != "real") 00403 { 00404 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00405 return false; 00406 } 00407 mmHeader.ctype=double_type; 00408 break; 00409 case 'c': 00410 /* sanity check */ 00411 if(buffer != "complex") 00412 { 00413 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00414 return false; 00415 } 00416 mmHeader.ctype=complex_type; 00417 break; 00418 case 'p': 00419 /* sanity check */ 00420 if(buffer != "pattern") 00421 { 00422 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00423 return false; 00424 } 00425 mmHeader.ctype=pattern; 00426 break; 00427 default: 00428 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00429 return false; 00430 } 00431 00432 if(lineFeed(file)) 00433 return false; 00434 00435 file >> buffer; 00436 00437 std::transform(buffer.begin(), buffer.end(), buffer.begin(), 00438 tolower); 00439 switch(buffer[0]) 00440 { 00441 case 'g': 00442 /* sanity check */ 00443 if(buffer != "general") 00444 { 00445 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00446 return false; 00447 } 00448 mmHeader.structure=general; 00449 break; 00450 case 'h': 00451 /* sanity check */ 00452 if(buffer != "hermitian") 00453 { 00454 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00455 return false; 00456 } 00457 mmHeader.structure=hermitian; 00458 break; 00459 case 's': 00460 if(buffer.size()==1){ 00461 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00462 return false; 00463 } 00464 00465 switch(buffer[1]) 00466 { 00467 case 'y': 00468 /* sanity check */ 00469 if(buffer != "symmetric") 00470 { 00471 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00472 return false; 00473 } 00474 mmHeader.structure=symmetric; 00475 break; 00476 case 'k': 00477 /* sanity check */ 00478 if(buffer != "skew-symmetric") 00479 { 00480 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00481 return false; 00482 } 00483 mmHeader.structure=skew_symmetric; 00484 break; 00485 default: 00486 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00487 return false; 00488 } 00489 default: 00490 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00491 return false; 00492 } 00493 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00494 c=file.peek(); 00495 return true; 00496 00497 } 00498 00499 void readNextLine(std::istream& file, std::ostringstream& line, LineType& type) 00500 { 00501 char c; 00502 std::size_t index=0; 00503 00504 //empty lines will be disgarded and we will simply read the next line 00505 while(index==0&&!file.eof()) 00506 { 00507 // strip spaces 00508 while(!file.eof() && (c=file.get())==' '); 00509 00510 //read the rest of the line until comment 00511 while(!file.eof() && (c=file.get())=='\n'){ 00512 switch(c) 00513 { 00514 case '%': 00515 /* disgard the rest of the line */ 00516 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00517 } 00518 } 00519 } 00520 00521 // buffer[index]='\0'; 00522 } 00523 00524 template<std::size_t brows, std::size_t bcols> 00525 Dune::tuple<std::size_t, std::size_t, std::size_t> 00526 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header) 00527 { 00528 std::size_t blockrows=rows/brows; 00529 std::size_t blockcols=cols/bcols; 00530 std::size_t blocksize=brows*bcols; 00531 std::size_t blockentries=0; 00532 00533 switch(header.structure) 00534 { 00535 case general: 00536 blockentries = entries/blocksize; break; 00537 case skew_symmetric: 00538 blockentries = 2*entries/blocksize; break; 00539 case symmetric: 00540 blockentries = (2*entries-rows)/blocksize; break; 00541 case hermitian: 00542 blockentries = (2*entries-rows)/blocksize; break; 00543 default: 00544 throw Dune::NotImplemented(); 00545 } 00546 return Dune::make_tuple(blockrows, blockcols, blockentries); 00547 } 00548 00549 /* 00550 * @brief Storage class for the column index and the numeric value. 00551 * 00552 * \tparam T Either a NumericWrapper of the numeric type or PatternDummy 00553 * for MatrixMarket pattern case. 00554 */ 00555 template<typename T> 00556 struct IndexData : public T 00557 { 00558 std::size_t index; 00559 }; 00560 00561 00572 template<typename T> 00573 struct NumericWrapper 00574 { 00575 T number; 00576 operator T&() 00577 { 00578 return number; 00579 } 00580 }; 00581 00585 struct PatternDummy 00586 {}; 00587 00588 template<> 00589 struct NumericWrapper<PatternDummy> 00590 {}; 00591 00592 template<typename T> 00593 std::istream& operator>>(std::istream& is, NumericWrapper<T>& num) 00594 { 00595 return is>>num.number; 00596 } 00597 00598 std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num) 00599 { 00600 return is; 00601 } 00602 00608 template<typename T> 00609 bool operator<(const IndexData<T>& i1, const IndexData<T>& i2) 00610 { 00611 return i1.index<i2.index; 00612 }; 00613 00619 template<typename T> 00620 std::istream& operator>>(std::istream& is, IndexData<T>& data) 00621 { 00622 is>>data.index; 00623 /* MatrixMarket indices are one based. Decrement for C++ */ 00624 --data.index; 00625 return is>>data.number; 00626 } 00627 00634 template<typename D, int brows, int bcols> 00635 struct MatrixValuesSetter 00636 { 00642 template<typename M> 00643 void operator()(const std::vector<std::set<IndexData<D> > >& rows, 00644 M& matrix) 00645 { 00646 for(typename M::RowIterator iter=matrix.begin(); 00647 iter!= matrix.end(); ++iter) 00648 { 00649 for(typename M::size_type brow=iter.index()*brows, 00650 browend=iter.index()*brows+brows; 00651 brow<browend; ++brow) 00652 { 00653 typedef typename std::set<IndexData<D> >::const_iterator Siter; 00654 for(Siter siter=rows[brow].begin(), send=rows[brow].end(); 00655 siter != send; ++siter) 00656 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number; 00657 } 00658 } 00659 } 00660 }; 00661 00662 template<int brows, int bcols> 00663 struct MatrixValuesSetter<PatternDummy,brows,bcols> 00664 { 00665 template<typename M> 00666 void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows, 00667 M& matrix) 00668 {} 00669 }; 00670 00671 template<typename T, typename A, int brows, int bcols, typename D> 00672 void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix, 00673 std::istream& file, std::size_t entries, 00674 const MMHeader& mmHeader, const D&) 00675 { 00676 typedef Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A> Matrix; 00677 // First path 00678 // store entries together with column index in a speparate 00679 // data structure 00680 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows); 00681 00682 for(; entries>0;--entries){ 00683 std::size_t row; 00684 IndexData<D> data; 00685 skipComments(file); 00686 file>>row; 00687 --row; // Index was 1 based. 00688 assert(row/bcols<matrix.N()); 00689 file>>data; 00690 assert(data.index/bcols<matrix.M()); 00691 rows[row].insert(data); 00692 } 00693 00694 // TODO extend to capture the nongeneral cases. 00695 if(mmHeader.structure!= general) 00696 DUNE_THROW(Dune::NotImplemented, "Only general is supported right now!"); 00697 00698 // Setup the matrix sparsity pattern 00699 int nnz=0; 00700 for(typename Matrix::CreateIterator iter=matrix.createbegin(); 00701 iter!= matrix.createend(); ++iter) 00702 { 00703 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows; 00704 brow<browend; ++brow) 00705 { 00706 typedef typename std::set<IndexData<D> >::const_iterator Siter; 00707 for(Siter siter=rows[brow].begin(), send=rows[brow].end(); 00708 siter != send; ++siter, ++nnz) 00709 iter.insert(siter->index/bcols); 00710 } 00711 } 00712 00713 //Set the matrix values 00714 matrix=0; 00715 00716 MatrixValuesSetter<D,brows,bcols> Setter; 00717 00718 Setter(rows, matrix); 00719 } 00720 } // end anonymous namespace 00721 00722 class MatrixMarketFormatError : public Dune::Exception 00723 {}; 00724 00725 00726 void mm_read_header(std::size_t& rows, std::size_t& cols, MMHeader& header, std::istream& istr, 00727 bool isVector) 00728 { 00729 if(!readMatrixMarketBanner(istr, header)){ 00730 std::cerr << "First line was not a correct Matrix Market banner. Using default:\n" 00731 << "%%MatrixMarket matrix coordinate real general"<<std::endl; 00732 // Go to the beginning of the file 00733 istr.clear() ; 00734 istr.seekg(0, std::ios::beg); 00735 if(isVector) 00736 header.type=array_type; 00737 } 00738 00739 skipComments(istr); 00740 00741 if(lineFeed(istr)) 00742 throw MatrixMarketFormatError(); 00743 00744 istr >> rows; 00745 00746 if(lineFeed(istr)) 00747 throw MatrixMarketFormatError(); 00748 istr >> cols; 00749 } 00750 00751 template<typename T, typename A, int entries> 00752 void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector, 00753 std::size_t size, 00754 std::istream& istr) 00755 { 00756 for(int i=0;size>0; ++i, --size){ 00757 T val; 00758 istr>>val; 00759 vector[i/entries][i%entries]=val; 00760 } 00761 } 00762 00763 00770 template<typename T, typename A, int entries> 00771 void readMatrixMarket(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector, 00772 std::istream& istr) 00773 { 00774 MMHeader header; 00775 std::size_t rows, cols; 00776 mm_read_header(rows,cols,header,istr, true); 00777 if(cols!=1) 00778 DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!"); 00779 00780 if(header.type!=array_type) 00781 DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!"); 00782 00783 std::size_t size=rows/entries; 00784 if(size*entries!=rows) 00785 DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct1"); 00786 00787 vector.resize(size); 00788 00789 istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00790 00791 mm_read_vector_entries(vector, rows, istr); 00792 } 00793 00794 00801 template<typename T, typename A, int brows, int bcols> 00802 void readMatrixMarket(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix, 00803 std::istream& istr) 00804 { 00805 00806 typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,brows,bcols> > Matrix; 00807 00808 MMHeader header; 00809 if(!readMatrixMarketBanner(istr, header)){ 00810 std::cerr << "First line was not a correct Matrix Market banner. Using default:\n" 00811 << "%%MatrixMarket matrix coordinate real general"<<std::endl; 00812 // Go to the beginning of the file 00813 istr.clear() ; 00814 istr.seekg(0, std::ios::beg); 00815 } 00816 skipComments(istr); 00817 00818 std::size_t rows, cols, entries; 00819 00820 if(lineFeed(istr)) 00821 throw MatrixMarketFormatError(); 00822 00823 istr >> rows; 00824 00825 if(lineFeed(istr)) 00826 throw MatrixMarketFormatError(); 00827 istr >> cols; 00828 00829 if(lineFeed(istr)) 00830 throw MatrixMarketFormatError(); 00831 00832 istr >>entries; 00833 00834 std::size_t nnz, blockrows, blockcols; 00835 00836 Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header); 00837 00838 istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n'); 00839 00840 00841 matrix.setSize(blockrows, blockcols); 00842 matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise); 00843 00844 if(header.type==array_type) 00845 DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!"); 00846 00847 readSparseEntries(matrix, istr, entries, header, NumericWrapper<double>()); 00848 } 00849 00850 template<typename M> 00851 struct mm_multipliers 00852 {}; 00853 00854 template<typename B, int i, int j, typename A> 00855 struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> > 00856 { 00857 enum{ 00858 rows = i, 00859 cols = j 00860 }; 00861 }; 00862 00863 template<typename B, int i, int j> 00864 void mm_print_entry(const FieldMatrix<B,i,j>& entry, 00865 typename FieldMatrix<B,i,j>::size_type rowidx, 00866 typename FieldMatrix<B,i,j>::size_type colidx, 00867 std::ostream& ostr) 00868 { 00869 typedef typename FieldMatrix<B,i,j>::const_iterator riterator; 00870 typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator; 00871 for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx){ 00872 int coli=colidx; 00873 for(citerator col = row->begin(); col != row->end(); ++col, ++coli) 00874 ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl; 00875 } 00876 } 00877 00878 // Write a vector entry 00879 template<typename V> 00880 void mm_print_vector_entry(const V& entry, std::ostream& ostr, 00881 const integral_constant<int,1>&) 00882 { 00883 ostr<<entry<<std::endl; 00884 } 00885 00886 // Write a vector 00887 template<typename V> 00888 void mm_print_vector_entry(const V& vector, std::ostream& ostr, 00889 const integral_constant<int,0>&) 00890 { 00891 // Is the entry a supported numeric type? 00892 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric; 00893 typedef typename V::const_iterator VIter; 00894 00895 for(VIter i=vector.begin(); i != vector.end(); ++i) 00896 00897 mm_print_vector_entry(*i, ostr, 00898 integral_constant<int,isnumeric>()); 00899 } 00900 00901 template<typename T, typename A, int i> 00902 std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector) 00903 { 00904 return vector.size()*i; 00905 } 00906 00907 // Version for writing vectors. 00908 template<typename V> 00909 void writeMatrixMarket(const V& vector, std::ostream& ostr, 00910 const integral_constant<int,0>&) 00911 { 00912 ostr<<countEntries(vector)<<" "<<1<<std::endl; 00913 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric; 00914 mm_print_vector_entry(vector,ostr, integral_constant<int,isnumeric>()); 00915 } 00916 00917 // Versions for writing matrices 00918 template<typename M> 00919 void writeMatrixMarket(const M& matrix, 00920 std::ostream& ostr, 00921 const integral_constant<int,1>&) 00922 { 00923 ostr<<matrix.M()*mm_multipliers<M>::rows<<" " 00924 <<matrix.N()*mm_multipliers<M>::cols<<" " 00925 <<countNonZeros(matrix)<<std::endl; 00926 00927 typedef typename M::const_iterator riterator; 00928 typedef typename M::ConstColIterator citerator; 00929 for(riterator row=matrix.begin(); row != matrix.end(); ++row) 00930 for(citerator col = row->begin(); col != row->end(); ++col) 00931 // Matrix Market indexing start with 1! 00932 mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1, 00933 col.index()*mm_multipliers<M>::cols+1, ostr); 00934 } 00935 00936 00940 template<typename M> 00941 void writeMatrixMarket(const M& matrix, 00942 std::ostream& ostr) 00943 { 00944 // Write header information 00945 mm_header_printer<M>::print(ostr); 00946 mm_block_structure_header<M>::print(ostr,matrix); 00947 // Choose the correct function for matrix and vector 00948 writeMatrixMarket(matrix,ostr,integral_constant<int,IsMatrix<M>::value>()); 00949 } 00950 00951 00962 template<typename M> 00963 void storeMatrixMarket(const M& matrix, 00964 std::string filename) 00965 { 00966 std::ofstream file(filename.c_str()); 00967 file.setf(std::ios::scientific,std::ios::floatfield); 00968 writeMatrixMarket(matrix, file); 00969 file.close(); 00970 } 00971 00972 #if HAVE_MPI 00973 00986 template<typename M, typename G, typename L> 00987 void storeMatrixMarket(const M& matrix, 00988 std::string filename, 00989 const OwnerOverlapCopyCommunication<G,L>& comm, 00990 bool storeIndices=true) 00991 { 00992 // Get our rank 00993 int rank = comm.communicator().rank(); 00994 // Write the local matrix 00995 std::ostringstream rfilename; 00996 rfilename<<filename <<"_"<<rank<<".mm"; 00997 std::cout<< rfilename.str()<<std::endl; 00998 std::ofstream file(rfilename.str().c_str()); 00999 file.setf(std::ios::scientific,std::ios::floatfield); 01000 writeMatrixMarket(matrix, file); 01001 file.close(); 01002 01003 if(!storeIndices) 01004 return; 01005 01006 // Write the global to local index mapping 01007 rfilename.str(""); 01008 rfilename<<filename<<"_"<<rank<<".idx"; 01009 file.open(rfilename.str().c_str()); 01010 file.setf(std::ios::scientific,std::ios::floatfield); 01011 typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet; 01012 typedef typename IndexSet::const_iterator Iterator; 01013 for(Iterator iter = comm.indexSet().begin(); 01014 iter != comm.indexSet().end(); ++iter){ 01015 file << iter->global()<<" "<<(std::size_t)iter->local()<<" " 01016 <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl; 01017 } 01018 // Store neighbour information for efficient remote indices setup. 01019 file<<"neighbours:"; 01020 const std::set<int>& neighbours=comm.remoteIndices().getNeighbours(); 01021 typedef std::set<int>::const_iterator SIter; 01022 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour){ 01023 file<<" "<< *neighbour; 01024 } 01025 file.close(); 01026 } 01027 01042 template<typename M, typename G, typename L> 01043 void loadMatrixMarket(M& matrix, 01044 const std::string& filename, 01045 OwnerOverlapCopyCommunication<G,L>& comm, 01046 bool readIndices=true) 01047 { 01048 typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet::LocalIndex LocalIndex; 01049 typedef typename LocalIndex::Attribute Attribute; 01050 // Get our rank 01051 int rank = comm.communicator().rank(); 01052 // load local matrix 01053 std::ostringstream rfilename; 01054 rfilename<<filename <<"_"<<rank<<".mm"; 01055 std::ifstream file; 01056 file.open(rfilename.str().c_str(), std::ios::in); 01057 if(!file) 01058 DUNE_THROW(IOError, "Could not open file" << rfilename.str().c_str()); 01059 //if(!file.is_open()) 01060 // 01061 readMatrixMarket(matrix,file); 01062 file.close(); 01063 01064 if(!readIndices) 01065 return; 01066 01067 // read indices 01068 typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet; 01069 IndexSet& pis=comm.pis; 01070 rfilename.str(""); 01071 rfilename<<filename<<"_"<<rank<<".idx"; 01072 file.open(rfilename.str().c_str()); 01073 if(pis.size()!=0) 01074 DUNE_THROW(InvalidIndexSetState, "Index set is not empty!"); 01075 01076 pis.beginResize(); 01077 while(!file.eof() && file.peek()!='n'){ 01078 G g; 01079 file >>g; 01080 std::size_t l; 01081 file >>l; 01082 int c; 01083 file >>c; 01084 bool b; 01085 file >> b; 01086 pis.add(g,LocalIndex(l,Attribute(c),b)); 01087 lineFeed(file); 01088 } 01089 pis.endResize(); 01090 if(!file.eof()){ 01091 // read neighbours 01092 std::string s; 01093 file>>s; 01094 if(s!="neighbours:") 01095 DUNE_THROW(MatrixMarketFormatError, "was exspecting the string: \"neighbours:\""); 01096 std::set<int> nb; 01097 while(!file.eof()){ 01098 int i; 01099 file >> i; 01100 nb.insert(i); 01101 } 01102 file.close(); 01103 comm.ri.setNeighbours(nb); 01104 } 01105 comm.ri.template rebuild<false>(); 01106 } 01107 01108 #endif 01109 01120 template<typename M> 01121 void loadMatrixMarket(M& matrix, 01122 const std::string& filename) 01123 { 01124 std::ifstream file; 01125 file.open(filename.c_str(), std::ios::in); 01126 if(!file) 01127 DUNE_THROW(IOError, "Could not open file" << filename); 01128 readMatrixMarket(matrix,file); 01129 file.close(); 01130 } 01131 01133 } 01134 #endif