Dune Core Modules (2.4.1)

matrixmarket.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_MATRIXMARKET_HH
4#define DUNE_ISTL_MATRIXMARKET_HH
5
6#include <ostream>
7#include <istream>
8#include <fstream>
9#include <sstream>
10#include <limits>
11#include <ios>
12#include "matrixutils.hh"
13#include "bcrsmatrix.hh"
14#include "owneroverlapcopy.hh"
16#include <dune/common/tuples.hh>
17#include <dune/common/unused.hh>
18
19namespace Dune
20{
21
47 namespace MatrixMarketImpl
48 {
58 template<class T>
60 enum {
64 is_numeric=false
65 };
66 };
67
68 template<>
69 struct mm_numeric_type<int>
70 {
71 enum {
75 is_numeric=true
76 };
77
78 static std::string str()
79 {
80 return "integer";
81 }
82 };
83
84 template<>
85 struct mm_numeric_type<double>
86 {
87 enum {
91 is_numeric=true
92 };
93
94 static std::string str()
95 {
96 return "real";
97 }
98 };
99
100 template<>
101 struct mm_numeric_type<float>
102 {
103 enum {
107 is_numeric=true
108 };
109
110 static std::string str()
111 {
112 return "real";
113 }
114 };
115
116 template<>
117 struct mm_numeric_type<std::complex<double> >
118 {
119 enum {
123 is_numeric=true
124 };
125
126 static std::string str()
127 {
128 return "complex";
129 }
130 };
131
132 template<>
133 struct mm_numeric_type<std::complex<float> >
134 {
135 enum {
139 is_numeric=true
140 };
141
142 static std::string str()
143 {
144 return "complex";
145 }
146 };
147
156 template<class M>
158
159 template<typename T, typename A, int i, int j>
160 struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> >
161 {
162 static void print(std::ostream& os)
163 {
164 os<<"%%MatrixMarket matrix coordinate ";
165 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
166 }
167 };
168
169 template<typename B, typename A>
170 struct mm_header_printer<BlockVector<B,A> >
171 {
172 static void print(std::ostream& os)
173 {
174 os<<"%%MatrixMarket matrix array ";
175 os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl;
176 }
177 };
178
179 template<typename T, int j>
180 struct mm_header_printer<FieldVector<T,j> >
181 {
182 static void print(std::ostream& os)
183 {
184 os<<"%%MatrixMarket matrix array ";
185 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
186 }
187 };
188
189 template<typename T, int i, int j>
190 struct mm_header_printer<FieldMatrix<T,i,j> >
191 {
192 static void print(std::ostream& os)
193 {
194 os<<"%%MatrixMarket matrix array ";
195 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
196 }
197 };
198
207 template<class M>
209
210 template<typename T, typename A, int i>
212 {
213 typedef BlockVector<FieldVector<T,i>,A> M;
214
215 static void print(std::ostream& os, const M&)
216 {
217 os<<"% ISTL_STRUCT blocked ";
218 os<<i<<" "<<1<<std::endl;
219 }
220 };
221
222 template<typename T, typename A, int i, int j>
223 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
224 {
226
227 static void print(std::ostream& os, const M&)
228 {
229 os<<"% ISTL_STRUCT blocked ";
230 os<<i<<" "<<j<<std::endl;
231 }
232 };
233
234
235 template<typename T, int i, int j>
236 struct mm_block_structure_header<FieldMatrix<T,i,j> >
237 {
238 typedef FieldMatrix<T,i,j> M;
239
240 static void print(std::ostream& os, const M& m)
241 {}
242 };
243
244 template<typename T, int i>
245 struct mm_block_structure_header<FieldVector<T,i> >
246 {
247 typedef FieldVector<T,i> M;
248
249 static void print(std::ostream& os, const M& m)
250 {}
251 };
252
253 enum LineType { MM_HEADER, MM_ISTLSTRUCT, DATA };
254 enum { MM_MAX_LINE_LENGTH=1025 };
255
256 enum MM_TYPE { coordinate_type, array_type, unknown_type };
257
258 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
259
260 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
261
262 struct MMHeader
263 {
264 MMHeader()
265 : type(coordinate_type), ctype(double_type), structure(general)
266 {}
267 MM_TYPE type;
268 MM_CTYPE ctype;
269 MM_STRUCTURE structure;
270 };
271
272 inline bool lineFeed(std::istream& file)
273 {
274 char c;
275 if(!file.eof())
276 c=file.peek();
277 else
278 return false;
279 // ignore whitespace
280 while(c==' ')
281 {
282 file.get();
283 if(file.eof())
284 return false;
285 c=file.peek();
286 }
287
288 if(c=='\n') {
289 /* eat the line feed */
290 file.get();
291 return true;
292 }
293 return false;
294 }
295
296 inline void skipComments(std::istream& file)
297 {
298 lineFeed(file);
299 char c=file.peek();
300 // ignore comment lines
301 while(c=='%')
302 {
303 /* disgard the rest of the line */
304 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
305 c=file.peek();
306 }
307 }
308
309
310 inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
311 {
312 std::string buffer;
313 char c;
314 file >> buffer;
315 c=buffer[0];
316 mmHeader=MMHeader();
317 if(c!='%')
318 return false;
319 std::cout<<buffer<<std::endl;
320 /* read the banner */
321 if(buffer!="%%MatrixMarket") {
322 /* disgard the rest of the line */
323 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
324 return false;
325 }
326
327 if(lineFeed(file))
328 /* premature end of line */
329 return false;
330
331 /* read the matrix_type */
332 file >> buffer;
333
334 if(buffer != "matrix")
335 {
336 /* disgard the rest of the line */
337 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
338 return false;
339 }
340
341 if(lineFeed(file))
342 /* premature end of line */
343 return false;
344
345 /* The type of the matrix */
346 file >> buffer;
347
348 if(buffer.empty())
349 return false;
350
351 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
352 tolower);
353
354 switch(buffer[0])
355 {
356 case 'a' :
357 /* sanity check */
358 if(buffer != "array")
359 {
360 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
361 return false;
362 }
363 mmHeader.type=array_type;
364 break;
365 case 'c' :
366 /* sanity check */
367 if(buffer != "coordinate")
368 {
369 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
370 return false;
371 }
372 mmHeader.type=coordinate_type;
373 break;
374 default :
375 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
376 return false;
377 }
378
379 if(lineFeed(file))
380 /* premature end of line */
381 return false;
382
383 /* The numeric type used. */
384 file >> buffer;
385
386 if(buffer.empty())
387 return false;
388
389 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
390 tolower);
391 switch(buffer[0])
392 {
393 case 'i' :
394 /* sanity check */
395 if(buffer != "integer")
396 {
397 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
398 return false;
399 }
400 mmHeader.ctype=integer_type;
401 break;
402 case 'r' :
403 /* sanity check */
404 if(buffer != "real")
405 {
406 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
407 return false;
408 }
409 mmHeader.ctype=double_type;
410 break;
411 case 'c' :
412 /* sanity check */
413 if(buffer != "complex")
414 {
415 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
416 return false;
417 }
418 mmHeader.ctype=complex_type;
419 break;
420 case 'p' :
421 /* sanity check */
422 if(buffer != "pattern")
423 {
424 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
425 return false;
426 }
427 mmHeader.ctype=pattern;
428 break;
429 default :
430 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
431 return false;
432 }
433
434 if(lineFeed(file))
435 return false;
436
437 file >> buffer;
438
439 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
440 tolower);
441 switch(buffer[0])
442 {
443 case 'g' :
444 /* sanity check */
445 if(buffer != "general")
446 {
447 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
448 return false;
449 }
450 mmHeader.structure=general;
451 break;
452 case 'h' :
453 /* sanity check */
454 if(buffer != "hermitian")
455 {
456 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
457 return false;
458 }
459 mmHeader.structure=hermitian;
460 break;
461 case 's' :
462 if(buffer.size()==1) {
463 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
464 return false;
465 }
466
467 switch(buffer[1])
468 {
469 case 'y' :
470 /* sanity check */
471 if(buffer != "symmetric")
472 {
473 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
474 return false;
475 }
476 mmHeader.structure=symmetric;
477 break;
478 case 'k' :
479 /* sanity check */
480 if(buffer != "skew-symmetric")
481 {
482 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
483 return false;
484 }
485 mmHeader.structure=skew_symmetric;
486 break;
487 default :
488 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
489 return false;
490 }
491 default :
492 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
493 return false;
494 }
495 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
496 c=file.peek();
497 return true;
498
499 }
500
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)
504 {
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;
509
510 switch(header.structure)
511 {
512 case general :
513 blockentries = entries/blocksize; break;
514 case skew_symmetric :
515 blockentries = 2*entries/blocksize; break;
516 case symmetric :
517 blockentries = (2*entries-rows)/blocksize; break;
518 case hermitian :
519 blockentries = (2*entries-rows)/blocksize; break;
520 default :
521 throw Dune::NotImplemented();
522 }
523 return Dune::make_tuple(blockrows, blockcols, blockentries);
524 }
525
526 /*
527 * @brief Storage class for the column index and the numeric value.
528 *
529 * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
530 * for MatrixMarket pattern case.
531 */
532 template<typename T>
533 struct IndexData : public T
534 {
535 std::size_t index;
536 };
537
538
549 template<typename T>
551 {
552 T number;
553 operator T&()
554 {
555 return number;
556 }
557 };
558
563 {};
564
565 template<>
567 {};
568
569 template<typename T>
570 std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
571 {
572 return is>>num.number;
573 }
574
575 inline std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
576 {
578 return is;
579 }
580
586 template<typename T>
587 bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
588 {
589 return i1.index<i2.index;
590 }
591
597 template<typename T>
598 std::istream& operator>>(std::istream& is, IndexData<T>& data)
599 {
600 is>>data.index;
601 /* MatrixMarket indices are one based. Decrement for C++ */
602 --data.index;
603 return is>>data.number;
604 }
605
612 template<typename D, int brows, int bcols>
614 {
620 template<typename M>
621 void operator()(const std::vector<std::set<IndexData<D> > >& rows,
622 M& matrix)
623 {
624 for(typename M::RowIterator iter=matrix.begin();
625 iter!= matrix.end(); ++iter)
626 {
627 for(typename M::size_type brow=iter.index()*brows,
628 browend=iter.index()*brows+brows;
629 brow<browend; ++brow)
630 {
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;
635 }
636 }
637 }
638 };
639
640 template<int brows, int bcols>
641 struct MatrixValuesSetter<PatternDummy,brows,bcols>
642 {
643 template<typename M>
644 void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
645 M& matrix)
646 {}
647 };
648
649 template<typename T, typename A, int brows, int bcols, typename D>
650 void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
651 std::istream& file, std::size_t entries,
652 const MMHeader& mmHeader, const D&)
653 {
655 // First path
656 // store entries together with column index in a speparate
657 // data structure
658 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
659
660 for(; entries>0; --entries) {
661 std::size_t row;
662 IndexData<D> data;
663 skipComments(file);
664 file>>row;
665 --row; // Index was 1 based.
666 assert(row/bcols<matrix.N());
667 file>>data;
668 assert(data.index/bcols<matrix.M());
669 rows[row].insert(data);
670 }
671
672 // TODO extend to capture the nongeneral cases.
673 if(mmHeader.structure!= general)
674 DUNE_THROW(Dune::NotImplemented, "Only general is supported right now!");
675
676 // Setup the matrix sparsity pattern
677 int nnz=0;
678 for(typename Matrix::CreateIterator iter=matrix.createbegin();
679 iter!= matrix.createend(); ++iter)
680 {
681 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
682 brow<browend; ++brow)
683 {
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);
688 }
689 }
690
691 //Set the matrix values
692 matrix=0;
693
694 MatrixValuesSetter<D,brows,bcols> Setter;
695
696 Setter(rows, matrix);
697 }
698 } // end namespace MatrixMarketImpl
699
700 class MatrixMarketFormatError : public Dune::Exception
701 {};
702
703
704 inline void mm_read_header(std::size_t& rows, std::size_t& cols,
705 MatrixMarketImpl::MMHeader& header, std::istream& istr,
706 bool isVector)
707 {
708 using namespace MatrixMarketImpl;
709
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;
713 // Go to the beginning of the file
714 istr.clear() ;
715 istr.seekg(0, std::ios::beg);
716 if(isVector)
717 header.type=array_type;
718 }
719
720 skipComments(istr);
721
722 if(lineFeed(istr))
723 throw MatrixMarketFormatError();
724
725 istr >> rows;
726
727 if(lineFeed(istr))
728 throw MatrixMarketFormatError();
729 istr >> cols;
730 }
731
732 template<typename T, typename A, int entries>
733 void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
734 std::size_t size,
735 std::istream& istr)
736 {
737 for(int i=0; size>0; ++i, --size) {
738 T val;
739 istr>>val;
740 vector[i/entries][i%entries]=val;
741 }
742 }
743
744
751 template<typename T, typename A, int entries>
753 std::istream& istr)
754 {
755 using namespace MatrixMarketImpl;
756
757 MMHeader header;
758 std::size_t rows, cols;
759 mm_read_header(rows,cols,header,istr, true);
760 if(cols!=1)
761 DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
762
763 if(header.type!=array_type)
764 DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
765
766 std::size_t size=rows/entries;
767 if(size*entries!=rows)
768 DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct1");
769
770 vector.resize(size);
771
772 istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
773
774 mm_read_vector_entries(vector, rows, istr);
775 }
776
777
784 template<typename T, typename A, int brows, int bcols>
786 std::istream& istr)
787 {
788 using namespace MatrixMarketImpl;
789
790 MMHeader header;
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;
794 // Go to the beginning of the file
795 istr.clear() ;
796 istr.seekg(0, std::ios::beg);
797 }
798 skipComments(istr);
799
800 std::size_t rows, cols, entries;
801
802 if(lineFeed(istr))
803 throw MatrixMarketFormatError();
804
805 istr >> rows;
806
807 if(lineFeed(istr))
808 throw MatrixMarketFormatError();
809 istr >> cols;
810
811 if(lineFeed(istr))
812 throw MatrixMarketFormatError();
813
814 istr >>entries;
815
816 std::size_t nnz, blockrows, blockcols;
817
818 Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
819
820 istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
821
822
823 matrix.setSize(blockrows, blockcols);
824 matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise);
825
826 if(header.type==array_type)
827 DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
828
829 readSparseEntries(matrix, istr, entries, header, NumericWrapper<T>());
830 }
831
832 template<typename M>
833 struct mm_multipliers
834 {};
835
836 template<typename B, int i, int j, typename A>
837 struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
838 {
839 enum {
840 rows = i,
841 cols = j
842 };
843 };
844
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,
849 std::ostream& ostr)
850 {
851 typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
852 typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
853 for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
854 int coli=colidx;
855 for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
856 ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
857 }
858 }
859
860 // Write a vector entry
861 template<typename V>
862 void mm_print_vector_entry(const V& entry, std::ostream& ostr,
863 const integral_constant<int,1>&)
864 {
865 ostr<<entry<<std::endl;
866 }
867
868 // Write a vector
869 template<typename V>
870 void mm_print_vector_entry(const V& vector, std::ostream& ostr,
871 const integral_constant<int,0>&)
872 {
873 using namespace MatrixMarketImpl;
874
875 // Is the entry a supported numeric type?
876 const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
877 typedef typename V::const_iterator VIter;
878
879 for(VIter i=vector.begin(); i != vector.end(); ++i)
880
881 mm_print_vector_entry(*i, ostr,
882 integral_constant<int,isnumeric>());
883 }
884
885 template<typename T, typename A, int i>
886 std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
887 {
888 return vector.size()*i;
889 }
890
891 // Version for writing vectors.
892 template<typename V>
893 void writeMatrixMarket(const V& vector, std::ostream& ostr,
894 const integral_constant<int,0>&)
895 {
896 using namespace MatrixMarketImpl;
897
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>());
901 }
902
903 // Versions for writing matrices
904 template<typename M>
905 void writeMatrixMarket(const M& matrix,
906 std::ostream& ostr,
907 const integral_constant<int,1>&)
908 {
909 ostr<<matrix.N()*mm_multipliers<M>::rows<<" "
910 <<matrix.M()*mm_multipliers<M>::cols<<" "
911 <<countNonZeros(matrix)<<std::endl;
912
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)
917 // Matrix Market indexing start with 1!
918 mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1,
919 col.index()*mm_multipliers<M>::cols+1, ostr);
920 }
921
922
926 template<typename M>
927 void writeMatrixMarket(const M& matrix,
928 std::ostream& ostr)
929 {
930 using namespace MatrixMarketImpl;
931
932 // Write header information
933 mm_header_printer<M>::print(ostr);
934 mm_block_structure_header<M>::print(ostr,matrix);
935 // Choose the correct function for matrix and vector
936 writeMatrixMarket(matrix,ostr,integral_constant<int,IsMatrix<M>::value>());
937 }
938
939
950 template<typename M>
951 void storeMatrixMarket(const M& matrix,
952 std::string filename)
953 {
954 std::ofstream file(filename.c_str());
955 file.setf(std::ios::scientific,std::ios::floatfield);
956 writeMatrixMarket(matrix, file);
957 file.close();
958 }
959
960#if HAVE_MPI
974 template<typename M, typename G, typename L>
975 void storeMatrixMarket(const M& matrix,
976 std::string filename,
978 bool storeIndices=true)
979 {
980 // Get our rank
981 int rank = comm.communicator().rank();
982 // Write the local matrix
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);
989 file.close();
990
991 if(!storeIndices)
992 return;
993
994 // Write the global to local index mapping
995 rfilename.str("");
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;
1001 for(Iterator iter = comm.indexSet().begin();
1002 iter != comm.indexSet().end(); ++iter) {
1003 file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1004 <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1005 }
1006 // Store neighbour information for efficient remote indices setup.
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;
1012 }
1013 file.close();
1014 }
1015
1030 template<typename M, typename G, typename L>
1031 void loadMatrixMarket(M& matrix,
1032 const std::string& filename,
1034 bool readIndices=true)
1035 {
1036 using namespace MatrixMarketImpl;
1037
1039 typedef typename LocalIndex::Attribute Attribute;
1040 // Get our rank
1041 int rank = comm.communicator().rank();
1042 // load local matrix
1043 std::ostringstream rfilename;
1044 rfilename<<filename <<"_"<<rank<<".mm";
1045 std::ifstream file;
1046 file.open(rfilename.str().c_str(), std::ios::in);
1047 if(!file)
1048 DUNE_THROW(IOError, "Could not open file" << rfilename.str().c_str());
1049 //if(!file.is_open())
1050 //
1051 readMatrixMarket(matrix,file);
1052 file.close();
1053
1054 if(!readIndices)
1055 return;
1056
1057 // read indices
1059 IndexSet& pis=comm.pis;
1060 rfilename.str("");
1061 rfilename<<filename<<"_"<<rank<<".idx";
1062 file.open(rfilename.str().c_str());
1063 if(pis.size()!=0)
1064 DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1065
1066 pis.beginResize();
1067 while(!file.eof() && file.peek()!='n') {
1068 G g;
1069 file >>g;
1070 std::size_t l;
1071 file >>l;
1072 int c;
1073 file >>c;
1074 bool b;
1075 file >> b;
1076 pis.add(g,LocalIndex(l,Attribute(c),b));
1077 lineFeed(file);
1078 }
1079 pis.endResize();
1080 if(!file.eof()) {
1081 // read neighbours
1082 std::string s;
1083 file>>s;
1084 if(s!="neighbours:")
1085 DUNE_THROW(MatrixMarketFormatError, "was exspecting the string: \"neighbours:\"");
1086 std::set<int> nb;
1087 while(!file.eof()) {
1088 int i;
1089 file >> i;
1090 nb.insert(i);
1091 }
1092 file.close();
1093 comm.ri.setNeighbours(nb);
1094 }
1095 comm.ri.template rebuild<false>();
1096 }
1097
1098 #endif
1099
1110 template<typename M>
1111 void loadMatrixMarket(M& matrix,
1112 const std::string& filename)
1113 {
1114 std::ifstream file;
1115 file.open(filename.c_str(), std::ios::in);
1116 if(!file)
1117 DUNE_THROW(IOError, "Could not open file" << filename);
1118 readMatrixMarket(matrix,file);
1119 file.close();
1120 }
1121
1123}
1124#endif
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
STL namespace.
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
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:208
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:157
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)