Dune Core Modules (2.8.0)

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 <algorithm>
7#include <complex>
8#include <cstddef>
9#include <fstream>
10#include <ios>
11#include <iostream>
12#include <istream>
13#include <limits>
14#include <ostream>
15#include <set>
16#include <sstream>
17#include <string>
18#include <tuple>
19#include <type_traits>
20#include <vector>
21
25#include <dune/common/hybridutilities.hh>
28
30#include <dune/istl/bvector.hh>
31#include <dune/istl/matrixutils.hh> // countNonZeros()
33
34namespace Dune
35{
36
62 namespace MatrixMarketImpl
63 {
73 template<class T>
75 enum {
79 is_numeric=false
80 };
81 };
82
83 template<>
84 struct mm_numeric_type<int>
85 {
86 enum {
90 is_numeric=true
91 };
92
93 static std::string str()
94 {
95 return "integer";
96 }
97 };
98
99 template<>
100 struct mm_numeric_type<double>
101 {
102 enum {
106 is_numeric=true
107 };
108
109 static std::string str()
110 {
111 return "real";
112 }
113 };
114
115 template<>
116 struct mm_numeric_type<float>
117 {
118 enum {
122 is_numeric=true
123 };
124
125 static std::string str()
126 {
127 return "real";
128 }
129 };
130
131 template<>
132 struct mm_numeric_type<std::complex<double> >
133 {
134 enum {
138 is_numeric=true
139 };
140
141 static std::string str()
142 {
143 return "complex";
144 }
145 };
146
147 template<>
148 struct mm_numeric_type<std::complex<float> >
149 {
150 enum {
154 is_numeric=true
155 };
156
157 static std::string str()
158 {
159 return "complex";
160 }
161 };
162
171 template<class M>
173
174 template<typename T, typename A>
175 struct mm_header_printer<BCRSMatrix<T,A> >
176 {
177 static void print(std::ostream& os)
178 {
179 os<<"%%MatrixMarket matrix coordinate ";
180 os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<T>::field_type>>::str()<<" general"<<std::endl;
181 }
182 };
183
184 template<typename B, typename A>
185 struct mm_header_printer<BlockVector<B,A> >
186 {
187 static void print(std::ostream& os)
188 {
189 os<<"%%MatrixMarket matrix array ";
190 os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<B>::field_type>>::str()<<" general"<<std::endl;
191 }
192 };
193
194 template<typename T, int j>
195 struct mm_header_printer<FieldVector<T,j> >
196 {
197 static void print(std::ostream& os)
198 {
199 os<<"%%MatrixMarket matrix array ";
200 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
201 }
202 };
203
204 template<typename T, int i, int j>
205 struct mm_header_printer<FieldMatrix<T,i,j> >
206 {
207 static void print(std::ostream& os)
208 {
209 os<<"%%MatrixMarket matrix array ";
210 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
211 }
212 };
213
222 template<class M>
224
225 template<typename T, typename A>
227 {
228 typedef BlockVector<T,A> M;
229 static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
230
231 static void print(std::ostream& os, const M&)
232 {
233 os<<"% ISTL_STRUCT blocked ";
234 os<<"1 1"<<std::endl;
235 }
236 };
237
238 template<typename T, typename A, int i>
239 struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
240 {
241 typedef BlockVector<FieldVector<T,i>,A> M;
242
243 static void print(std::ostream& os, const M&)
244 {
245 os<<"% ISTL_STRUCT blocked ";
246 os<<i<<" "<<1<<std::endl;
247 }
248 };
249
250 template<typename T, typename A>
251 struct mm_block_structure_header<BCRSMatrix<T,A> >
252 {
253 typedef BCRSMatrix<T,A> M;
254 static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
255
256 static void print(std::ostream& os, const M&)
257 {
258 os<<"% ISTL_STRUCT blocked ";
259 os<<"1 1"<<std::endl;
260 }
261 };
262
263 template<typename T, typename A, int i, int j>
264 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
265 {
266 typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M;
267
268 static void print(std::ostream& os, const M&)
269 {
270 os<<"% ISTL_STRUCT blocked ";
271 os<<i<<" "<<j<<std::endl;
272 }
273 };
274
275
276 template<typename T, int i, int j>
277 struct mm_block_structure_header<FieldMatrix<T,i,j> >
278 {
279 typedef FieldMatrix<T,i,j> M;
280
281 static void print(std::ostream& os, const M& m)
282 {}
283 };
284
285 template<typename T, int i>
286 struct mm_block_structure_header<FieldVector<T,i> >
287 {
288 typedef FieldVector<T,i> M;
289
290 static void print(std::ostream& os, const M& m)
291 {}
292 };
293
294 enum LineType { MM_HEADER, MM_ISTLSTRUCT, DATA };
295 enum { MM_MAX_LINE_LENGTH=1025 };
296
297 enum MM_TYPE { coordinate_type, array_type, unknown_type };
298
299 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
300
301 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
302
303 struct MMHeader
304 {
305 MMHeader()
306 : type(coordinate_type), ctype(double_type), structure(general)
307 {}
308 MM_TYPE type;
309 MM_CTYPE ctype;
310 MM_STRUCTURE structure;
311 };
312
313 inline bool lineFeed(std::istream& file)
314 {
315 char c;
316 if(!file.eof())
317 c=file.peek();
318 else
319 return false;
320 // ignore whitespace
321 while(c==' ')
322 {
323 file.get();
324 if(file.eof())
325 return false;
326 c=file.peek();
327 }
328
329 if(c=='\n') {
330 /* eat the line feed */
331 file.get();
332 return true;
333 }
334 return false;
335 }
336
337 inline void skipComments(std::istream& file)
338 {
339 lineFeed(file);
340 char c=file.peek();
341 // ignore comment lines
342 while(c=='%')
343 {
344 /* discard the rest of the line */
346 c=file.peek();
347 }
348 }
349
350
351 inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
352 {
353 std::string buffer;
354 char c;
355 file >> buffer;
356 c=buffer[0];
357 mmHeader=MMHeader();
358 if(c!='%')
359 return false;
360 dverb<<buffer<<std::endl;
361 /* read the banner */
362 if(buffer!="%%MatrixMarket") {
363 /* discard the rest of the line */
365 return false;
366 }
367
368 if(lineFeed(file))
369 /* premature end of line */
370 return false;
371
372 /* read the matrix_type */
373 file >> buffer;
374
375 if(buffer != "matrix")
376 {
377 /* discard the rest of the line */
379 return false;
380 }
381
382 if(lineFeed(file))
383 /* premature end of line */
384 return false;
385
386 /* The type of the matrix */
387 file >> buffer;
388
389 if(buffer.empty())
390 return false;
391
392 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
393 ::tolower);
394
395 switch(buffer[0])
396 {
397 case 'a' :
398 /* sanity check */
399 if(buffer != "array")
400 {
402 return false;
403 }
404 mmHeader.type=array_type;
405 break;
406 case 'c' :
407 /* sanity check */
408 if(buffer != "coordinate")
409 {
411 return false;
412 }
413 mmHeader.type=coordinate_type;
414 break;
415 default :
417 return false;
418 }
419
420 if(lineFeed(file))
421 /* premature end of line */
422 return false;
423
424 /* The numeric type used. */
425 file >> buffer;
426
427 if(buffer.empty())
428 return false;
429
430 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
431 ::tolower);
432 switch(buffer[0])
433 {
434 case 'i' :
435 /* sanity check */
436 if(buffer != "integer")
437 {
439 return false;
440 }
441 mmHeader.ctype=integer_type;
442 break;
443 case 'r' :
444 /* sanity check */
445 if(buffer != "real")
446 {
448 return false;
449 }
450 mmHeader.ctype=double_type;
451 break;
452 case 'c' :
453 /* sanity check */
454 if(buffer != "complex")
455 {
457 return false;
458 }
459 mmHeader.ctype=complex_type;
460 break;
461 case 'p' :
462 /* sanity check */
463 if(buffer != "pattern")
464 {
466 return false;
467 }
468 mmHeader.ctype=pattern;
469 break;
470 default :
472 return false;
473 }
474
475 if(lineFeed(file))
476 return false;
477
478 file >> buffer;
479
480 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
481 ::tolower);
482 switch(buffer[0])
483 {
484 case 'g' :
485 /* sanity check */
486 if(buffer != "general")
487 {
489 return false;
490 }
491 mmHeader.structure=general;
492 break;
493 case 'h' :
494 /* sanity check */
495 if(buffer != "hermitian")
496 {
498 return false;
499 }
500 mmHeader.structure=hermitian;
501 break;
502 case 's' :
503 if(buffer.size()==1) {
505 return false;
506 }
507
508 switch(buffer[1])
509 {
510 case 'y' :
511 /* sanity check */
512 if(buffer != "symmetric")
513 {
515 return false;
516 }
517 mmHeader.structure=symmetric;
518 break;
519 case 'k' :
520 /* sanity check */
521 if(buffer != "skew-symmetric")
522 {
524 return false;
525 }
526 mmHeader.structure=skew_symmetric;
527 break;
528 default :
530 return false;
531 }
532 break;
533 default :
535 return false;
536 }
538 c=file.peek();
539 return true;
540
541 }
542
543 template<std::size_t brows, std::size_t bcols>
544 std::tuple<std::size_t, std::size_t, std::size_t>
545 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
546 {
547 std::size_t blockrows=rows/brows;
548 std::size_t blockcols=cols/bcols;
549 std::size_t blocksize=brows*bcols;
550 std::size_t blockentries=0;
551
552 switch(header.structure)
553 {
554 case general :
555 blockentries = entries/blocksize; break;
556 case skew_symmetric :
557 blockentries = 2*entries/blocksize; break;
558 case symmetric :
559 blockentries = (2*entries-rows)/blocksize; break;
560 case hermitian :
561 blockentries = (2*entries-rows)/blocksize; break;
562 default :
563 throw Dune::NotImplemented();
564 }
565 return std::make_tuple(blockrows, blockcols, blockentries);
566 }
567
568 /*
569 * @brief Storage class for the column index and the numeric value.
570 *
571 * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
572 * for MatrixMarket pattern case.
573 */
574 template<typename T>
575 struct IndexData : public T
576 {
577 std::size_t index;
578 };
579
580
591 template<typename T>
593 {
594 T number;
595 operator T&()
596 {
597 return number;
598 }
599 };
600
605 {};
606
607 template<>
609 {};
610
611 template<typename T>
612 std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
613 {
614 return is>>num.number;
615 }
616
617 inline std::istream& operator>>(std::istream& is, [[maybe_unused]] NumericWrapper<PatternDummy>& num)
618 {
619 return is;
620 }
621
627 template<typename T>
628 bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
629 {
630 return i1.index<i2.index;
631 }
632
638 template<typename T>
639 std::istream& operator>>(std::istream& is, IndexData<T>& data)
640 {
641 is>>data.index;
642 /* MatrixMarket indices are one based. Decrement for C++ */
643 --data.index;
644 return is>>data.number;
645 }
646
652 template<typename T>
653 std::istream& operator>>(std::istream& is, IndexData<NumericWrapper<std::complex<T>>>& data)
654 {
655 is>>data.index;
656 /* MatrixMarket indices are one based. Decrement for C++ */
657 --data.index;
658 // real and imaginary part needs to be read separately as
659 // complex numbers are not provided in pair form. (x,y)
660 NumericWrapper<T> real, imag;
661 is>>real;
662 is>>imag;
663 data.number = {real.number, imag.number};
664 return is;
665 }
666
673 template<typename D, int brows, int bcols>
675 {
681 template<typename T>
682 void operator()(const std::vector<std::set<IndexData<D> > >& rows,
683 BCRSMatrix<T>& matrix)
684 {
685 static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
686 for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
687 {
688 auto brow=iter.index();
689 for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
690 (*iter)[siter->index] = siter->number;
691 }
692 }
693
699 template<typename T>
700 void operator()(const std::vector<std::set<IndexData<D> > >& rows,
702 {
703 for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
704 {
705 for (auto brow=iter.index()*brows,
706 browend=iter.index()*brows+brows;
707 brow<browend; ++brow)
708 {
709 for (auto siter=rows[brow].begin(), send=rows[brow].end();
710 siter != send; ++siter)
711 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
712 }
713 }
714 }
715 };
716
717 template<int brows, int bcols>
718 struct MatrixValuesSetter<PatternDummy,brows,bcols>
719 {
720 template<typename M>
721 void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
722 M& matrix)
723 {}
724 };
725
726 template<class T> struct is_complex : std::false_type {};
727 template<class T> struct is_complex<std::complex<T>> : std::true_type {};
728
729 // wrapper for std::conj. Returns T if T is not complex.
730 template<class T>
731 std::enable_if_t<!is_complex<T>::value, T> conj(const T& r){
732 return r;
733 }
734
735 template<class T>
736 std::enable_if_t<is_complex<T>::value, T> conj(const T& r){
737 return std::conj(r);
738 }
739
740 template<typename M>
741 struct mm_multipliers
742 {};
743
744 template<typename B, typename A>
745 struct mm_multipliers<BCRSMatrix<B,A> >
746 {
747 enum {
748 rows = 1,
749 cols = 1
750 };
751 };
752
753 template<typename B, int i, int j, typename A>
754 struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
755 {
756 enum {
757 rows = i,
758 cols = j
759 };
760 };
761
762 template<typename T, typename A, typename D>
763 void readSparseEntries(Dune::BCRSMatrix<T,A>& matrix,
764 std::istream& file, std::size_t entries,
765 const MMHeader& mmHeader, const D&)
766 {
767 typedef Dune::BCRSMatrix<T,A> Matrix;
768
769 // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
770 constexpr int brows = mm_multipliers<Matrix>::rows;
771 constexpr int bcols = mm_multipliers<Matrix>::cols;
772
773 // First path
774 // store entries together with column index in a separate
775 // data structure
776 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
777
778 auto readloop = [&] (auto symmetryFixup) {
779 for(std::size_t i = 0; i < entries; ++i) {
780 std::size_t row;
781 IndexData<D> data;
782 skipComments(file);
783 file>>row;
784 --row; // Index was 1 based.
785 assert(row/bcols<matrix.N());
786 file>>data;
787 assert(data.index/bcols<matrix.M());
788 rows[row].insert(data);
789 if(row!=data.index)
790 symmetryFixup(row, data);
791 }
792 };
793
794 switch(mmHeader.structure)
795 {
796 case general:
797 readloop([](auto...){});
798 break;
799 case symmetric :
800 readloop([&](auto row, auto data) {
801 IndexData<D> data_sym(data);
802 data_sym.index = row;
803 rows[data.index].insert(data_sym);
804 });
805 break;
806 case skew_symmetric :
807 readloop([&](auto row, auto data) {
808 IndexData<D> data_sym;
809 data_sym.number = -data.number;
810 data_sym.index = row;
811 rows[data.index].insert(data_sym);
812 });
813 break;
814 case hermitian :
815 readloop([&](auto row, auto data) {
816 IndexData<D> data_sym;
817 data_sym.number = conj(data.number);
818 data_sym.index = row;
819 rows[data.index].insert(data_sym);
820 });
821 break;
822 default:
824 "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
825 }
826
827 // Setup the matrix sparsity pattern
828 int nnz=0;
829 for(typename Matrix::CreateIterator iter=matrix.createbegin();
830 iter!= matrix.createend(); ++iter)
831 {
832 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
833 brow<browend; ++brow)
834 {
835 typedef typename std::set<IndexData<D> >::const_iterator Siter;
836 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
837 siter != send; ++siter, ++nnz)
838 iter.insert(siter->index/bcols);
839 }
840 }
841
842 //Set the matrix values
843 matrix=0;
844
845 MatrixValuesSetter<D,brows,bcols> Setter;
846
847 Setter(rows, matrix);
848 }
849 } // end namespace MatrixMarketImpl
850
851 class MatrixMarketFormatError : public Dune::Exception
852 {};
853
854
855 inline void mm_read_header(std::size_t& rows, std::size_t& cols,
856 MatrixMarketImpl::MMHeader& header, std::istream& istr,
857 bool isVector)
858 {
859 using namespace MatrixMarketImpl;
860
861 if(!readMatrixMarketBanner(istr, header)) {
862 std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
863 << "%%MatrixMarket matrix coordinate real general"<<std::endl;
864 // Go to the beginning of the file
865 istr.clear() ;
866 istr.seekg(0, std::ios::beg);
867 if(isVector)
868 header.type=array_type;
869 }
870
871 skipComments(istr);
872
873 if(lineFeed(istr))
874 throw MatrixMarketFormatError();
875
876 istr >> rows;
877
878 if(lineFeed(istr))
879 throw MatrixMarketFormatError();
880 istr >> cols;
881 }
882
883 template<typename T, typename A>
884 void mm_read_vector_entries(Dune::BlockVector<T,A>& vector,
885 std::size_t size,
886 std::istream& istr,
887 size_t lane)
888 {
889 for (int i=0; size>0; ++i, --size)
890 istr>>Simd::lane(lane,vector[i]);
891 }
892
893 template<typename T, typename A, int entries>
894 void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
895 std::size_t size,
896 std::istream& istr,
897 size_t lane)
898 {
899 for(int i=0; size>0; ++i, --size) {
900 Simd::Scalar<T> val;
901 istr>>val;
902 Simd::lane(lane, vector[i/entries][i%entries])=val;
903 }
904 }
905
906
913 template<typename T, typename A>
915 std::istream& istr)
916 {
917 typedef typename Dune::BlockVector<T,A>::field_type field_type;
918 using namespace MatrixMarketImpl;
919
920 MMHeader header;
921 std::size_t rows, cols;
922 mm_read_header(rows,cols,header,istr, true);
923 if(cols!=Simd::lanes<field_type>()) {
924 if(Simd::lanes<field_type>() == 1)
925 DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
926 else
927 DUNE_THROW(MatrixMarketFormatError, "cols does not match the number of lanes in the field_type!");
928 }
929
930 if(header.type!=array_type)
931 DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
932
933
934 if constexpr (Dune::IsNumber<T>())
935 vector.resize(rows);
936 else
937 {
938 T dummy;
939 auto blocksize = dummy.size();
940 std::size_t size=rows/blocksize;
941 if(size*blocksize!=rows)
942 DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
943
944 vector.resize(size);
945 }
946
948 for(size_t l=0;l<Simd::lanes<field_type>();++l){
949 mm_read_vector_entries(vector, rows, istr, l);
950 }
951 }
952
959 template<typename T, typename A>
961 std::istream& istr)
962 {
963 using namespace MatrixMarketImpl;
965
966 MMHeader header;
967 if(!readMatrixMarketBanner(istr, header)) {
968 std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
969 << "%%MatrixMarket matrix coordinate real general"<<std::endl;
970 // Go to the beginning of the file
971 istr.clear() ;
972 istr.seekg(0, std::ios::beg);
973 }
974 skipComments(istr);
975
976 std::size_t rows, cols, entries;
977
978 if(lineFeed(istr))
979 throw MatrixMarketFormatError();
980
981 istr >> rows;
982
983 if(lineFeed(istr))
984 throw MatrixMarketFormatError();
985 istr >> cols;
986
987 if(lineFeed(istr))
988 throw MatrixMarketFormatError();
989
990 istr >>entries;
991
992 std::size_t nnz, blockrows, blockcols;
993
994 // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
995 constexpr int brows = mm_multipliers<Matrix>::rows;
996 constexpr int bcols = mm_multipliers<Matrix>::cols;
997
998 std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
999
1001
1002
1003 matrix.setSize(blockrows, blockcols);
1005
1006 if(header.type==array_type)
1007 DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
1008
1009 readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1010 }
1011
1012 // Print a scalar entry
1013 template<typename B>
1014 void mm_print_entry(const B& entry,
1015 std::size_t rowidx,
1016 std::size_t colidx,
1017 std::ostream& ostr)
1018 {
1019 if constexpr (IsNumber<B>())
1020 ostr << rowidx << " " << colidx << " " << entry << std::endl;
1021 else
1022 {
1023 for (auto row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
1024 int coli=colidx;
1025 for (auto col = row->begin(); col != row->end(); ++col, ++coli)
1026 ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
1027 }
1028 }
1029 }
1030
1031 // Write a vector entry
1032 template<typename V>
1033 void mm_print_vector_entry(const V& entry, std::ostream& ostr,
1034 const std::integral_constant<int,1>&,
1035 size_t lane)
1036 {
1037 ostr<<Simd::lane(lane,entry)<<std::endl;
1038 }
1039
1040 // Write a vector
1041 template<typename V>
1042 void mm_print_vector_entry(const V& vector, std::ostream& ostr,
1043 const std::integral_constant<int,0>&,
1044 size_t lane)
1045 {
1046 using namespace MatrixMarketImpl;
1047
1048 // Is the entry a supported numeric type?
1049 const int isnumeric = mm_numeric_type<Simd::Scalar<typename V::block_type>>::is_numeric;
1050 typedef typename V::const_iterator VIter;
1051
1052 for(VIter i=vector.begin(); i != vector.end(); ++i)
1053
1054 mm_print_vector_entry(*i, ostr,
1055 std::integral_constant<int,isnumeric>(),
1056 lane);
1057 }
1058
1059 template<typename T, typename A>
1060 std::size_t countEntries(const BlockVector<T,A>& vector)
1061 {
1062 return vector.size();
1063 }
1064
1065 template<typename T, typename A, int i>
1066 std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
1067 {
1068 return vector.size()*i;
1069 }
1070
1071 // Version for writing vectors.
1072 template<typename V>
1073 void writeMatrixMarket(const V& vector, std::ostream& ostr,
1074 const std::integral_constant<int,0>&)
1075 {
1076 using namespace MatrixMarketImpl;
1077 typedef typename V::field_type field_type;
1078
1079 ostr<<countEntries(vector)<<" "<<Simd::lanes<field_type>()<<std::endl;
1080 const int isnumeric = mm_numeric_type<Simd::Scalar<V>>::is_numeric;
1081 for(size_t l=0;l<Simd::lanes<field_type>(); ++l){
1082 mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>(), l);
1083 }
1084 }
1085
1086 // Versions for writing matrices
1087 template<typename M>
1088 void writeMatrixMarket(const M& matrix,
1089 std::ostream& ostr,
1090 const std::integral_constant<int,1>&)
1091 {
1092 ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
1093 <<matrix.M()*MatrixMarketImpl::mm_multipliers<M>::cols<<" "
1094 <<countNonZeros(matrix)<<std::endl;
1095
1096 typedef typename M::const_iterator riterator;
1097 typedef typename M::ConstColIterator citerator;
1098 for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1099 for(citerator col = row->begin(); col != row->end(); ++col)
1100 // Matrix Market indexing start with 1!
1101 mm_print_entry(*col, row.index()*MatrixMarketImpl::mm_multipliers<M>::rows+1,
1102 col.index()*MatrixMarketImpl::mm_multipliers<M>::cols+1, ostr);
1103 }
1104
1105
1109 template<typename M>
1110 void writeMatrixMarket(const M& matrix,
1111 std::ostream& ostr)
1112 {
1113 using namespace MatrixMarketImpl;
1114
1115 // Write header information
1116 mm_header_printer<M>::print(ostr);
1117 mm_block_structure_header<M>::print(ostr,matrix);
1118 // Choose the correct function for matrix and vector
1119 writeMatrixMarket(matrix,ostr,std::integral_constant<int,IsMatrix<M>::value>());
1120 }
1121
1122 static const int default_precision = -1;
1134 template<typename M>
1135 void storeMatrixMarket(const M& matrix,
1136 std::string filename,
1137 int prec=default_precision)
1138 {
1139 std::ofstream file(filename.c_str());
1140 file.setf(std::ios::scientific,std::ios::floatfield);
1141 if(prec>0)
1142 file.precision(prec);
1143 writeMatrixMarket(matrix, file);
1144 file.close();
1145 }
1146
1147#if HAVE_MPI
1162 template<typename M, typename G, typename L>
1163 void storeMatrixMarket(const M& matrix,
1164 std::string filename,
1166 bool storeIndices=true,
1167 int prec=default_precision)
1168 {
1169 // Get our rank
1170 int rank = comm.communicator().rank();
1171 // Write the local matrix
1172 std::ostringstream rfilename;
1173 rfilename<<filename <<"_"<<rank<<".mm";
1174 dverb<< rfilename.str()<<std::endl;
1175 std::ofstream file(rfilename.str().c_str());
1176 file.setf(std::ios::scientific,std::ios::floatfield);
1177 if(prec>0)
1178 file.precision(prec);
1179 writeMatrixMarket(matrix, file);
1180 file.close();
1181
1182 if(!storeIndices)
1183 return;
1184
1185 // Write the global to local index mapping
1186 rfilename.str("");
1187 rfilename<<filename<<"_"<<rank<<".idx";
1188 file.open(rfilename.str().c_str());
1189 file.setf(std::ios::scientific,std::ios::floatfield);
1191 typedef typename IndexSet::const_iterator Iterator;
1192 for(Iterator iter = comm.indexSet().begin();
1193 iter != comm.indexSet().end(); ++iter) {
1194 file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1195 <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1196 }
1197 // Store neighbour information for efficient remote indices setup.
1198 file<<"neighbours:";
1199 const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1200 typedef std::set<int>::const_iterator SIter;
1201 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1202 file<<" "<< *neighbour;
1203 }
1204 file.close();
1205 }
1206
1221 template<typename M, typename G, typename L>
1222 void loadMatrixMarket(M& matrix,
1223 const std::string& filename,
1225 bool readIndices=true)
1226 {
1227 using namespace MatrixMarketImpl;
1228
1230 typedef typename LocalIndexT::Attribute Attribute;
1231 // Get our rank
1232 int rank = comm.communicator().rank();
1233 // load local matrix
1234 std::ostringstream rfilename;
1235 rfilename<<filename <<"_"<<rank<<".mm";
1236 std::ifstream file;
1237 file.open(rfilename.str().c_str(), std::ios::in);
1238 if(!file)
1239 DUNE_THROW(IOError, "Could not open file: " << rfilename.str().c_str());
1240 //if(!file.is_open())
1241 //
1242 readMatrixMarket(matrix,file);
1243 file.close();
1244
1245 if(!readIndices)
1246 return;
1247
1248 // read indices
1250 IndexSet& pis=comm.pis;
1251 rfilename.str("");
1252 rfilename<<filename<<"_"<<rank<<".idx";
1253 file.open(rfilename.str().c_str());
1254 if(pis.size()!=0)
1255 DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1256
1257 pis.beginResize();
1258 while(!file.eof() && file.peek()!='n') {
1259 G g;
1260 file >>g;
1261 std::size_t l;
1262 file >>l;
1263 int c;
1264 file >>c;
1265 bool b;
1266 file >> b;
1267 pis.add(g,LocalIndexT(l,Attribute(c),b));
1268 lineFeed(file);
1269 }
1270 pis.endResize();
1271 if(!file.eof()) {
1272 // read neighbours
1273 std::string s;
1274 file>>s;
1275 if(s!="neighbours:")
1276 DUNE_THROW(MatrixMarketFormatError, "was expecting the string: \"neighbours:\"");
1277 std::set<int> nb;
1278 while(!file.eof()) {
1279 int i;
1280 file >> i;
1281 nb.insert(i);
1282 }
1283 file.close();
1284 comm.ri.setNeighbours(nb);
1285 }
1286 comm.ri.template rebuild<false>();
1287 }
1288
1289 #endif
1290
1301 template<typename M>
1302 void loadMatrixMarket(M& matrix,
1303 const std::string& filename)
1304 {
1305 std::ifstream file;
1306 file.open(filename.c_str(), std::ios::in);
1307 if(!file)
1308 DUNE_THROW(IOError, "Could not open file: " << filename);
1309 readMatrixMarket(matrix,file);
1310 file.close();
1311 }
1312
1314}
1315#endif
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:673
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:679
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1101
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1976
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1095
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:831
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:859
A vector of blocks with memory management.
Definition: bvector.hh:393
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:501
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:399
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Default exception class for I/O errors.
Definition: exceptions.hh:229
Index Set Interface base class.
Definition: indexidset.hh:76
auto size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:221
Exception indicating that the index set is not in the expected state.
Definition: indexset.hh:203
A generic dynamic dense matrix.
Definition: matrix.hh:559
Default exception for dummy implementations.
Definition: exceptions.hh:261
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:460
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:469
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
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:237
Stream & operator>>(Stream &stream, std::tuple< Ts... > &t)
Read a std::tuple.
Definition: streamoperators.hh:41
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:914
void storeMatrixMarket(const M &matrix, std::string filename, int prec=default_precision)
Stores a parallel matrix/vector in matrix market format in a file.
Definition: matrixmarket.hh:1135
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:1222
auto countNonZeros(const M &, typename std::enable_if_t< Dune::IsNumber< M >::value > *sfinae=nullptr)
Get the number of nonzero fields in the matrix.
Definition: matrixutils.hh:117
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:635
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition: interface.hh:322
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:114
Some handy generic functions for ISTL matrices.
Dune namespace.
Definition: alignedallocator.hh:11
T lane(std::size_t l, const T &v)
access a lane of a simd vector (scalar version)
Definition: simd.hh:364
STL namespace.
Classes providing communication interfaces for overlapping Schwarz methods.
Include file for users of the SIMD abstraction layer.
Standard Dune debug streams.
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Functor to the data values of the matrix.
Definition: matrixmarket.hh:675
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:682
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:700
a wrapper class of numeric values.
Definition: matrixmarket.hh:593
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:605
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:223
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:172
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:74
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:79
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)