dune-istl 2.8.0
Loading...
Searching...
No Matches
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
22#include <dune/common/exceptions.hh>
23#include <dune/common/fmatrix.hh>
24#include <dune/common/fvector.hh>
25#include <dune/common/hybridutilities.hh>
26#include <dune/common/stdstreams.hh>
27#include <dune/common/simd/simd.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>
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>
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>
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 {
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 {
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>
252 {
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>
265 {
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>
278 {
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
295 enum { MM_MAX_LINE_LENGTH=1025 };
296
298
300
302
303 struct MMHeader
304 {
307 {}
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 */
345 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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 */
364 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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 */
378 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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 {
401 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
402 return false;
403 }
404 mmHeader.type=array_type;
405 break;
406 case 'c' :
407 /* sanity check */
408 if(buffer != "coordinate")
409 {
410 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
411 return false;
412 }
413 mmHeader.type=coordinate_type;
414 break;
415 default :
416 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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 {
438 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
439 return false;
440 }
441 mmHeader.ctype=integer_type;
442 break;
443 case 'r' :
444 /* sanity check */
445 if(buffer != "real")
446 {
447 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
448 return false;
449 }
450 mmHeader.ctype=double_type;
451 break;
452 case 'c' :
453 /* sanity check */
454 if(buffer != "complex")
455 {
456 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
457 return false;
458 }
459 mmHeader.ctype=complex_type;
460 break;
461 case 'p' :
462 /* sanity check */
463 if(buffer != "pattern")
464 {
465 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
466 return false;
467 }
468 mmHeader.ctype=pattern;
469 break;
470 default :
471 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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 {
488 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
489 return false;
490 }
491 mmHeader.structure=general;
492 break;
493 case 'h' :
494 /* sanity check */
495 if(buffer != "hermitian")
496 {
497 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
498 return false;
499 }
500 mmHeader.structure=hermitian;
501 break;
502 case 's' :
503 if(buffer.size()==1) {
504 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
505 return false;
506 }
507
508 switch(buffer[1])
509 {
510 case 'y' :
511 /* sanity check */
512 if(buffer != "symmetric")
513 {
514 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
515 return false;
516 }
517 mmHeader.structure=symmetric;
518 break;
519 case 'k' :
520 /* sanity check */
521 if(buffer != "skew-symmetric")
522 {
523 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
524 return false;
525 }
526 mmHeader.structure=skew_symmetric;
527 break;
528 default :
529 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
530 return false;
531 }
532 break;
533 default :
534 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
535 return false;
536 }
537 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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 {
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>
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>
742 {};
743
744 template<typename B, typename A>
746 {
747 enum {
748 rows = 1,
749 cols = 1
750 };
751 };
752
753 template<typename B, int i, int j, typename A>
755 {
756 enum {
757 rows = i,
758 cols = j
759 };
760 };
761
762 template<typename T, typename A, typename D>
764 std::istream& file, std::size_t entries,
765 const MMHeader& mmHeader, const D&)
766 {
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:
823 DUNE_THROW(Dune::NotImplemented,
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
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))
875
876 istr >> rows;
877
878 if(lineFeed(istr))
880 istr >> cols;
881 }
882
883 template<typename T, typename A>
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
947 istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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))
980
981 istr >> rows;
982
983 if(lineFeed(istr))
985 istr >> cols;
986
987 if(lineFeed(istr))
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
1000 istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
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<<" "
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!
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
Some handy generic functions for ISTL matrices.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Classes providing communication interfaces for overlapping Schwarz methods.
Implementation of the BCRSMatrix class.
Col col
Definition: matrixmatrix.hh:349
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
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
std::size_t countEntries(const BlockVector< T, A > &vector)
Definition: matrixmarket.hh:1060
void writeMatrixMarket(const V &vector, std::ostream &ostr, const std::integral_constant< int, 0 > &)
Definition: matrixmarket.hh:1073
void mm_print_vector_entry(const V &entry, std::ostream &ostr, const std::integral_constant< int, 1 > &, size_t lane)
Definition: matrixmarket.hh:1033
static const int default_precision
Definition: matrixmarket.hh:1122
void mm_read_vector_entries(Dune::BlockVector< T, A > &vector, std::size_t size, std::istream &istr, size_t lane)
Definition: matrixmarket.hh:884
void mm_read_header(std::size_t &rows, std::size_t &cols, MatrixMarketImpl::MMHeader &header, std::istream &istr, bool isVector)
Definition: matrixmarket.hh:855
void mm_print_entry(const B &entry, std::size_t rowidx, std::size_t colidx, std::ostream &ostr)
Definition: matrixmarket.hh:1014
Definition: allocator.hh:9
std::tuple< std::size_t, std::size_t, std::size_t > calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader &header)
Definition: matrixmarket.hh:545
bool operator<(const IndexData< T > &i1, const IndexData< T > &i2)
LessThan operator.
Definition: matrixmarket.hh:628
LineType
Definition: matrixmarket.hh:294
@ DATA
Definition: matrixmarket.hh:294
@ MM_HEADER
Definition: matrixmarket.hh:294
@ MM_ISTLSTRUCT
Definition: matrixmarket.hh:294
bool readMatrixMarketBanner(std::istream &file, MMHeader &mmHeader)
Definition: matrixmarket.hh:351
void readSparseEntries(Dune::BCRSMatrix< T, A > &matrix, std::istream &file, std::size_t entries, const MMHeader &mmHeader, const D &)
Definition: matrixmarket.hh:763
MM_TYPE
Definition: matrixmarket.hh:297
@ array_type
Definition: matrixmarket.hh:297
@ coordinate_type
Definition: matrixmarket.hh:297
@ unknown_type
Definition: matrixmarket.hh:297
std::istream & operator>>(std::istream &is, NumericWrapper< T > &num)
Definition: matrixmarket.hh:612
void skipComments(std::istream &file)
Definition: matrixmarket.hh:337
bool lineFeed(std::istream &file)
Definition: matrixmarket.hh:313
@ MM_MAX_LINE_LENGTH
Definition: matrixmarket.hh:295
MM_STRUCTURE
Definition: matrixmarket.hh:301
@ skew_symmetric
Definition: matrixmarket.hh:301
@ general
Definition: matrixmarket.hh:301
@ hermitian
Definition: matrixmarket.hh:301
@ unknown_structure
Definition: matrixmarket.hh:301
@ symmetric
Definition: matrixmarket.hh:301
MM_CTYPE
Definition: matrixmarket.hh:299
@ unknown_ctype
Definition: matrixmarket.hh:299
@ pattern
Definition: matrixmarket.hh:299
@ complex_type
Definition: matrixmarket.hh:299
@ double_type
Definition: matrixmarket.hh:299
@ integer_type
Definition: matrixmarket.hh:299
std::enable_if_t<!is_complex< T >::value, T > conj(const T &r)
Definition: matrixmarket.hh:731
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
A generic dynamic dense matrix.
Definition: matrix.hh:559
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
static std::string str()
Definition: matrixmarket.hh:93
static std::string str()
Definition: matrixmarket.hh:109
static std::string str()
Definition: matrixmarket.hh:125
static std::string str()
Definition: matrixmarket.hh:141
static std::string str()
Definition: matrixmarket.hh:157
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:172
static void print(std::ostream &os)
Definition: matrixmarket.hh:177
static void print(std::ostream &os)
Definition: matrixmarket.hh:187
static void print(std::ostream &os)
Definition: matrixmarket.hh:197
static void print(std::ostream &os)
Definition: matrixmarket.hh:207
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:223
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:231
BlockVector< T, A > M
Definition: matrixmarket.hh:228
BlockVector< FieldVector< T, i >, A > M
Definition: matrixmarket.hh:241
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:243
BCRSMatrix< T, A > M
Definition: matrixmarket.hh:253
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:256
BCRSMatrix< FieldMatrix< T, i, j >, A > M
Definition: matrixmarket.hh:266
static void print(std::ostream &os, const M &)
Definition: matrixmarket.hh:268
static void print(std::ostream &os, const M &m)
Definition: matrixmarket.hh:281
FieldMatrix< T, i, j > M
Definition: matrixmarket.hh:279
static void print(std::ostream &os, const M &m)
Definition: matrixmarket.hh:290
FieldVector< T, i > M
Definition: matrixmarket.hh:288
Definition: matrixmarket.hh:304
MM_STRUCTURE structure
Definition: matrixmarket.hh:310
MM_TYPE type
Definition: matrixmarket.hh:308
MMHeader()
Definition: matrixmarket.hh:305
MM_CTYPE ctype
Definition: matrixmarket.hh:309
Definition: matrixmarket.hh:576
std::size_t index
Definition: matrixmarket.hh:577
a wrapper class of numeric values.
Definition: matrixmarket.hh:593
T number
Definition: matrixmarket.hh:594
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:605
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
void operator()(const std::vector< std::set< IndexData< PatternDummy > > > &rows, M &matrix)
Definition: matrixmarket.hh:721
Definition: matrixmarket.hh:726
Definition: matrixmarket.hh:742
Definition: matrixmarket.hh:852
Definition: matrixutils.hh:25
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:502
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 CollectiveCommunication< MPI_Comm > & communicator() const
Definition: owneroverlapcopy.hh:297
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:469
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:447