DUNE-FEM (unstable)

matrixmarket.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_MATRIXMARKET_HH
6#define DUNE_ISTL_MATRIXMARKET_HH
7
8#include <algorithm>
9#include <complex>
10#include <cstddef>
11#include <fstream>
12#include <ios>
13#include <iostream>
14#include <istream>
15#include <limits>
16#include <ostream>
17#include <set>
18#include <sstream>
19#include <string>
20#include <tuple>
21#include <type_traits>
22#include <vector>
23
27#include <dune/common/hybridutilities.hh>
30
32#include <dune/istl/bvector.hh>
33#include <dune/istl/matrixutils.hh> // countNonZeros()
35
36namespace Dune
37{
38
64 namespace MatrixMarketImpl
65 {
75 template<class T>
77 enum {
81 is_numeric=false
82 };
83 };
84
85 template<>
86 struct mm_numeric_type<int>
87 {
88 enum {
92 is_numeric=true
93 };
94
95 static std::string str()
96 {
97 return "integer";
98 }
99 };
100
101 template<>
102 struct mm_numeric_type<double>
103 {
104 enum {
108 is_numeric=true
109 };
110
111 static std::string str()
112 {
113 return "real";
114 }
115 };
116
117 template<>
118 struct mm_numeric_type<float>
119 {
120 enum {
124 is_numeric=true
125 };
126
127 static std::string str()
128 {
129 return "real";
130 }
131 };
132
133 template<>
134 struct mm_numeric_type<std::complex<double> >
135 {
136 enum {
140 is_numeric=true
141 };
142
143 static std::string str()
144 {
145 return "complex";
146 }
147 };
148
149 template<>
150 struct mm_numeric_type<std::complex<float> >
151 {
152 enum {
156 is_numeric=true
157 };
158
159 static std::string str()
160 {
161 return "complex";
162 }
163 };
164
173 template<class M>
175
176 template<typename T, typename A>
177 struct mm_header_printer<BCRSMatrix<T,A> >
178 {
179 static void print(std::ostream& os)
180 {
181 os<<"%%MatrixMarket matrix coordinate ";
182 os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<T>::field_type>>::str()<<" general"<<std::endl;
183 }
184 };
185
186 template<typename B, typename A>
187 struct mm_header_printer<BlockVector<B,A> >
188 {
189 static void print(std::ostream& os)
190 {
191 os<<"%%MatrixMarket matrix array ";
192 os<<mm_numeric_type<Simd::Scalar<typename Imp::BlockTraits<B>::field_type>>::str()<<" general"<<std::endl;
193 }
194 };
195
196 template<typename T, int j>
197 struct mm_header_printer<FieldVector<T,j> >
198 {
199 static void print(std::ostream& os)
200 {
201 os<<"%%MatrixMarket matrix array ";
202 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
203 }
204 };
205
206 template<typename T, int i, int j>
207 struct mm_header_printer<FieldMatrix<T,i,j> >
208 {
209 static void print(std::ostream& os)
210 {
211 os<<"%%MatrixMarket matrix array ";
212 os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
213 }
214 };
215
224 template<class M>
226
227 template<typename T, typename A>
229 {
230 typedef BlockVector<T,A> M;
231 static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
232
233 static void print(std::ostream& os, const M&)
234 {
235 os<<"% ISTL_STRUCT blocked ";
236 os<<"1 1"<<std::endl;
237 }
238 };
239
240 template<typename T, typename A, int i>
241 struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
242 {
243 typedef BlockVector<FieldVector<T,i>,A> M;
244
245 static void print(std::ostream& os, const M&)
246 {
247 os<<"% ISTL_STRUCT blocked ";
248 os<<i<<" "<<1<<std::endl;
249 }
250 };
251
252 template<typename T, typename A>
253 struct mm_block_structure_header<BCRSMatrix<T,A> >
254 {
255 typedef BCRSMatrix<T,A> M;
256 static_assert(IsNumber<T>::value, "Only scalar entries are expected here!");
257
258 static void print(std::ostream& os, const M&)
259 {
260 os<<"% ISTL_STRUCT blocked ";
261 os<<"1 1"<<std::endl;
262 }
263 };
264
265 template<typename T, typename A, int i, int j>
266 struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
267 {
268 typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M;
269
270 static void print(std::ostream& os, const M&)
271 {
272 os<<"% ISTL_STRUCT blocked ";
273 os<<i<<" "<<j<<std::endl;
274 }
275 };
276
277
278 template<typename T, int i, int j>
279 struct mm_block_structure_header<FieldMatrix<T,i,j> >
280 {
281 typedef FieldMatrix<T,i,j> M;
282
283 static void print(std::ostream& os, const M& m)
284 {}
285 };
286
287 template<typename T, int i>
288 struct mm_block_structure_header<FieldVector<T,i> >
289 {
290 typedef FieldVector<T,i> M;
291
292 static void print(std::ostream& os, const M& m)
293 {}
294 };
295
296 enum LineType { MM_HEADER, MM_ISTLSTRUCT, DATA };
297 enum { MM_MAX_LINE_LENGTH=1025 };
298
299 enum MM_TYPE { coordinate_type, array_type, unknown_type };
300
301 enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
302
303 enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure };
304
305 struct MMHeader
306 {
307 MMHeader()
308 : type(coordinate_type), ctype(double_type), structure(general)
309 {}
310 MM_TYPE type;
311 MM_CTYPE ctype;
312 MM_STRUCTURE structure;
313 };
314
315 inline bool lineFeed(std::istream& file)
316 {
317 char c;
318 if(!file.eof())
319 c=file.peek();
320 else
321 return false;
322 // ignore whitespace
323 while(c==' ')
324 {
325 file.get();
326 if(file.eof())
327 return false;
328 c=file.peek();
329 }
330
331 if(c=='\n') {
332 /* eat the line feed */
333 file.get();
334 return true;
335 }
336 return false;
337 }
338
339 inline void skipComments(std::istream& file)
340 {
341 lineFeed(file);
342 char c=file.peek();
343 // ignore comment lines
344 while(c=='%')
345 {
346 /* discard the rest of the line */
348 c=file.peek();
349 }
350 }
351
352
353 inline bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
354 {
355 std::string buffer;
356 char c;
357 file >> buffer;
358 c=buffer[0];
359 mmHeader=MMHeader();
360 if(c!='%')
361 return false;
362 dverb<<buffer<<std::endl;
363 /* read the banner */
364 if(buffer!="%%MatrixMarket") {
365 /* discard the rest of the line */
367 return false;
368 }
369
370 if(lineFeed(file))
371 /* premature end of line */
372 return false;
373
374 /* read the matrix_type */
375 file >> buffer;
376
377 if(buffer != "matrix")
378 {
379 /* discard the rest of the line */
381 return false;
382 }
383
384 if(lineFeed(file))
385 /* premature end of line */
386 return false;
387
388 /* The type of the matrix */
389 file >> buffer;
390
391 if(buffer.empty())
392 return false;
393
394 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
395 ::tolower);
396
397 switch(buffer[0])
398 {
399 case 'a' :
400 /* sanity check */
401 if(buffer != "array")
402 {
404 return false;
405 }
406 mmHeader.type=array_type;
407 break;
408 case 'c' :
409 /* sanity check */
410 if(buffer != "coordinate")
411 {
413 return false;
414 }
415 mmHeader.type=coordinate_type;
416 break;
417 default :
419 return false;
420 }
421
422 if(lineFeed(file))
423 /* premature end of line */
424 return false;
425
426 /* The numeric type used. */
427 file >> buffer;
428
429 if(buffer.empty())
430 return false;
431
432 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
433 ::tolower);
434 switch(buffer[0])
435 {
436 case 'i' :
437 /* sanity check */
438 if(buffer != "integer")
439 {
441 return false;
442 }
443 mmHeader.ctype=integer_type;
444 break;
445 case 'r' :
446 /* sanity check */
447 if(buffer != "real")
448 {
450 return false;
451 }
452 mmHeader.ctype=double_type;
453 break;
454 case 'c' :
455 /* sanity check */
456 if(buffer != "complex")
457 {
459 return false;
460 }
461 mmHeader.ctype=complex_type;
462 break;
463 case 'p' :
464 /* sanity check */
465 if(buffer != "pattern")
466 {
468 return false;
469 }
470 mmHeader.ctype=pattern;
471 break;
472 default :
474 return false;
475 }
476
477 if(lineFeed(file))
478 return false;
479
480 file >> buffer;
481
482 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
483 ::tolower);
484 switch(buffer[0])
485 {
486 case 'g' :
487 /* sanity check */
488 if(buffer != "general")
489 {
491 return false;
492 }
493 mmHeader.structure=general;
494 break;
495 case 'h' :
496 /* sanity check */
497 if(buffer != "hermitian")
498 {
500 return false;
501 }
502 mmHeader.structure=hermitian;
503 break;
504 case 's' :
505 if(buffer.size()==1) {
507 return false;
508 }
509
510 switch(buffer[1])
511 {
512 case 'y' :
513 /* sanity check */
514 if(buffer != "symmetric")
515 {
517 return false;
518 }
519 mmHeader.structure=symmetric;
520 break;
521 case 'k' :
522 /* sanity check */
523 if(buffer != "skew-symmetric")
524 {
526 return false;
527 }
528 mmHeader.structure=skew_symmetric;
529 break;
530 default :
532 return false;
533 }
534 break;
535 default :
537 return false;
538 }
540 c=file.peek();
541 return true;
542
543 }
544
545 template<std::size_t brows, std::size_t bcols>
546 std::tuple<std::size_t, std::size_t, std::size_t>
547 calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
548 {
549 std::size_t blockrows=rows/brows;
550 std::size_t blockcols=cols/bcols;
551 std::size_t blocksize=brows*bcols;
552 std::size_t blockentries=0;
553
554 switch(header.structure)
555 {
556 case general :
557 blockentries = entries/blocksize; break;
558 case skew_symmetric :
559 blockentries = 2*entries/blocksize; break;
560 case symmetric :
561 blockentries = (2*entries-rows)/blocksize; break;
562 case hermitian :
563 blockentries = (2*entries-rows)/blocksize; break;
564 default :
565 throw Dune::NotImplemented();
566 }
567 return std::make_tuple(blockrows, blockcols, blockentries);
568 }
569
570 /*
571 * @brief Storage class for the column index and the numeric value.
572 *
573 * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
574 * for MatrixMarket pattern case.
575 */
576 template<typename T>
577 struct IndexData : public T
578 {
579 std::size_t index = {};
580 };
581
582
593 template<typename T>
595 {
596 T number = {};
597 operator T&()
598 {
599 return number;
600 }
601 };
602
607 {};
608
609 template<>
611 {};
612
613 template<typename T>
614 std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
615 {
616 return is>>num.number;
617 }
618
619 inline std::istream& operator>>(std::istream& is, [[maybe_unused]] NumericWrapper<PatternDummy>& num)
620 {
621 return is;
622 }
623
629 template<typename T>
630 bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
631 {
632 return i1.index<i2.index;
633 }
634
640 template<typename T>
641 std::istream& operator>>(std::istream& is, IndexData<T>& data)
642 {
643 is>>data.index;
644 /* MatrixMarket indices are one based. Decrement for C++ */
645 --data.index;
646 return is>>data.number;
647 }
648
654 template<typename T>
655 std::istream& operator>>(std::istream& is, IndexData<NumericWrapper<std::complex<T>>>& data)
656 {
657 is>>data.index;
658 /* MatrixMarket indices are one based. Decrement for C++ */
659 --data.index;
660 // real and imaginary part needs to be read separately as
661 // complex numbers are not provided in pair form. (x,y)
662 NumericWrapper<T> real, imag;
663 is>>real;
664 is>>imag;
665 data.number = {real.number, imag.number};
666 return is;
667 }
668
675 template<typename D, int brows, int bcols>
677 {
683 template<typename T>
684 void operator()(const std::vector<std::set<IndexData<D> > >& rows,
685 BCRSMatrix<T>& matrix)
686 {
687 static_assert(IsNumber<T>::value && brows==1 && bcols==1, "Only scalar entries are expected here!");
688 for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
689 {
690 auto brow=iter.index();
691 for (auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
692 (*iter)[siter->index] = siter->number;
693 }
694 }
695
701 template<typename T>
702 void operator()(const std::vector<std::set<IndexData<D> > >& rows,
704 {
705 for (auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
706 {
707 for (auto brow=iter.index()*brows,
708 browend=iter.index()*brows+brows;
709 brow<browend; ++brow)
710 {
711 for (auto siter=rows[brow].begin(), send=rows[brow].end();
712 siter != send; ++siter)
713 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
714 }
715 }
716 }
717 };
718
719 template<int brows, int bcols>
720 struct MatrixValuesSetter<PatternDummy,brows,bcols>
721 {
722 template<typename M>
723 void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
724 M& matrix)
725 {}
726 };
727
728 template<class T> struct is_complex : std::false_type {};
729 template<class T> struct is_complex<std::complex<T>> : std::true_type {};
730
731 // wrapper for std::conj. Returns T if T is not complex.
732 template<class T>
733 std::enable_if_t<!is_complex<T>::value, T> conj(const T& r){
734 return r;
735 }
736
737 template<class T>
738 std::enable_if_t<is_complex<T>::value, T> conj(const T& r){
739 return std::conj(r);
740 }
741
742 template<typename M>
743 struct mm_multipliers
744 {};
745
746 template<typename B, typename A>
747 struct mm_multipliers<BCRSMatrix<B,A> >
748 {
749 enum {
750 rows = 1,
751 cols = 1
752 };
753 };
754
755 template<typename B, int i, int j, typename A>
756 struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
757 {
758 enum {
759 rows = i,
760 cols = j
761 };
762 };
763
764 template<typename T, typename A, typename D>
765 void readSparseEntries(Dune::BCRSMatrix<T,A>& matrix,
766 std::istream& file, std::size_t entries,
767 const MMHeader& mmHeader, const D&)
768 {
769 typedef Dune::BCRSMatrix<T,A> Matrix;
770
771 // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
772 constexpr int brows = mm_multipliers<Matrix>::rows;
773 constexpr int bcols = mm_multipliers<Matrix>::cols;
774
775 // First path
776 // store entries together with column index in a separate
777 // data structure
778 std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
779
780 auto readloop = [&] (auto symmetryFixup) {
781 for(std::size_t i = 0; i < entries; ++i) {
782 std::size_t row;
783 IndexData<D> data;
784 skipComments(file);
785 file>>row;
786 --row; // Index was 1 based.
787 assert(row/bcols<matrix.N());
788 file>>data;
789 assert(data.index/bcols<matrix.M());
790 rows[row].insert(data);
791 if(row!=data.index)
792 symmetryFixup(row, data);
793 }
794 };
795
796 switch(mmHeader.structure)
797 {
798 case general:
799 readloop([](auto...){});
800 break;
801 case symmetric :
802 readloop([&](auto row, auto data) {
803 IndexData<D> data_sym(data);
804 data_sym.index = row;
805 rows[data.index].insert(data_sym);
806 });
807 break;
808 case skew_symmetric :
809 readloop([&](auto row, auto data) {
810 IndexData<D> data_sym;
811 data_sym.number = -data.number;
812 data_sym.index = row;
813 rows[data.index].insert(data_sym);
814 });
815 break;
816 case hermitian :
817 readloop([&](auto row, auto data) {
818 IndexData<D> data_sym;
819 data_sym.number = conj(data.number);
820 data_sym.index = row;
821 rows[data.index].insert(data_sym);
822 });
823 break;
824 default:
826 "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
827 }
828
829 // Setup the matrix sparsity pattern
830 int nnz=0;
831 for(typename Matrix::CreateIterator iter=matrix.createbegin();
832 iter!= matrix.createend(); ++iter)
833 {
834 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
835 brow<browend; ++brow)
836 {
837 typedef typename std::set<IndexData<D> >::const_iterator Siter;
838 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
839 siter != send; ++siter, ++nnz)
840 iter.insert(siter->index/bcols);
841 }
842 }
843
844 //Set the matrix values
845 matrix=0;
846
847 MatrixValuesSetter<D,brows,bcols> Setter;
848
849 Setter(rows, matrix);
850 }
851
852 inline std::tuple<std::string, std::string> splitFilename(const std::string& filename) {
853 std::size_t lastdot = filename.find_last_of(".");
854 if(lastdot == std::string::npos)
855 return std::make_tuple(filename, "");
856 else {
857 std::string potentialFileExtension = filename.substr(lastdot);
858 if (potentialFileExtension == ".mm" || potentialFileExtension == ".mtx")
859 return std::make_tuple(filename.substr(0, lastdot), potentialFileExtension);
860 else
861 return std::make_tuple(filename, "");
862 }
863 }
864
865 } // end namespace MatrixMarketImpl
866
867 class MatrixMarketFormatError : public Dune::Exception
868 {};
869
870
871 inline void mm_read_header(std::size_t& rows, std::size_t& cols,
872 MatrixMarketImpl::MMHeader& header, std::istream& istr,
873 bool isVector)
874 {
875 using namespace MatrixMarketImpl;
876
877 if(!readMatrixMarketBanner(istr, header)) {
878 std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
879 << "%%MatrixMarket matrix coordinate real general"<<std::endl;
880 // Go to the beginning of the file
881 istr.clear() ;
882 istr.seekg(0, std::ios::beg);
883 if(isVector)
884 header.type=array_type;
885 }
886
887 skipComments(istr);
888
889 if(lineFeed(istr))
890 throw MatrixMarketFormatError();
891
892 istr >> rows;
893
894 if(lineFeed(istr))
895 throw MatrixMarketFormatError();
896 istr >> cols;
897 }
898
899 template<typename T, typename A>
900 void mm_read_vector_entries(Dune::BlockVector<T,A>& vector,
901 std::size_t size,
902 std::istream& istr,
903 size_t lane)
904 {
905 for (int i=0; size>0; ++i, --size)
906 istr>>Simd::lane(lane,vector[i]);
907 }
908
909 template<typename T, typename A, int entries>
910 void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
911 std::size_t size,
912 std::istream& istr,
913 size_t lane)
914 {
915 for(int i=0; size>0; ++i, --size) {
916 Simd::Scalar<T> val;
917 istr>>val;
918 Simd::lane(lane, vector[i/entries][i%entries])=val;
919 }
920 }
921
922
929 template<typename T, typename A>
931 std::istream& istr)
932 {
933 typedef typename Dune::BlockVector<T,A>::field_type field_type;
934 using namespace MatrixMarketImpl;
935
936 MMHeader header;
937 std::size_t rows, cols;
938 mm_read_header(rows,cols,header,istr, true);
939 if(cols!=Simd::lanes<field_type>()) {
940 if(Simd::lanes<field_type>() == 1)
941 DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
942 else
943 DUNE_THROW(MatrixMarketFormatError, "cols does not match the number of lanes in the field_type!");
944 }
945
946 if(header.type!=array_type)
947 DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
948
949
950 if constexpr (Dune::IsNumber<T>())
951 vector.resize(rows);
952 else
953 {
954 T dummy;
955 auto blocksize = dummy.size();
956 std::size_t size=rows/blocksize;
957 if(size*blocksize!=rows)
958 DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct!");
959
960 vector.resize(size);
961 }
962
964 for(size_t l=0;l<Simd::lanes<field_type>();++l){
965 mm_read_vector_entries(vector, rows, istr, l);
966 }
967 }
968
975 template<typename T, typename A>
977 std::istream& istr)
978 {
979 using namespace MatrixMarketImpl;
981
982 MMHeader header;
983 if(!readMatrixMarketBanner(istr, header)) {
984 std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
985 << "%%MatrixMarket matrix coordinate real general"<<std::endl;
986 // Go to the beginning of the file
987 istr.clear() ;
988 istr.seekg(0, std::ios::beg);
989 }
990 skipComments(istr);
991
992 std::size_t rows, cols, entries;
993
994 if(lineFeed(istr))
995 throw MatrixMarketFormatError();
996
997 istr >> rows;
998
999 if(lineFeed(istr))
1000 throw MatrixMarketFormatError();
1001 istr >> cols;
1002
1003 if(lineFeed(istr))
1004 throw MatrixMarketFormatError();
1005
1006 istr >>entries;
1007
1008 std::size_t nnz, blockrows, blockcols;
1009
1010 // Number of rows and columns of T, if it is a matrix (1x1 otherwise)
1011 constexpr int brows = mm_multipliers<Matrix>::rows;
1012 constexpr int bcols = mm_multipliers<Matrix>::cols;
1013
1014 std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
1015
1017
1018
1019 matrix.setSize(blockrows, blockcols, nnz);
1021
1022 if(header.type==array_type)
1023 DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
1024
1025 readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1026 }
1027
1028 // Print a scalar entry
1029 template<typename B>
1030 void mm_print_entry(const B& entry,
1031 std::size_t rowidx,
1032 std::size_t colidx,
1033 std::ostream& ostr)
1034 {
1035 if constexpr (IsNumber<B>())
1036 ostr << rowidx << " " << colidx << " " << entry << std::endl;
1037 else
1038 {
1039 for (auto row=entry.begin(); row != entry.end(); ++row, ++rowidx) {
1040 int coli=colidx;
1041 for (auto col = row->begin(); col != row->end(); ++col, ++coli)
1042 ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
1043 }
1044 }
1045 }
1046
1047 // Write a vector entry
1048 template<typename V>
1049 void mm_print_vector_entry(const V& entry, std::ostream& ostr,
1050 const std::integral_constant<int,1>&,
1051 size_t lane)
1052 {
1053 ostr<<Simd::lane(lane,entry)<<std::endl;
1054 }
1055
1056 // Write a vector
1057 template<typename V>
1058 void mm_print_vector_entry(const V& vector, std::ostream& ostr,
1059 const std::integral_constant<int,0>&,
1060 size_t lane)
1061 {
1062 using namespace MatrixMarketImpl;
1063
1064 // Is the entry a supported numeric type?
1065 const int isnumeric = mm_numeric_type<Simd::Scalar<typename V::block_type>>::is_numeric;
1066 typedef typename V::const_iterator VIter;
1067
1068 for(VIter i=vector.begin(); i != vector.end(); ++i)
1069
1070 mm_print_vector_entry(*i, ostr,
1071 std::integral_constant<int,isnumeric>(),
1072 lane);
1073 }
1074
1075 template<typename T, typename A>
1076 std::size_t countEntries(const BlockVector<T,A>& vector)
1077 {
1078 return vector.size();
1079 }
1080
1081 template<typename T, typename A, int i>
1082 std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
1083 {
1084 return vector.size()*i;
1085 }
1086
1087 // Version for writing vectors.
1088 template<typename V>
1089 void writeMatrixMarket(const V& vector, std::ostream& ostr,
1090 const std::integral_constant<int,0>&)
1091 {
1092 using namespace MatrixMarketImpl;
1093 typedef typename V::field_type field_type;
1094
1095 ostr<<countEntries(vector)<<" "<<Simd::lanes<field_type>()<<std::endl;
1096 const int isnumeric = mm_numeric_type<Simd::Scalar<V>>::is_numeric;
1097 for(size_t l=0;l<Simd::lanes<field_type>(); ++l){
1098 mm_print_vector_entry(vector,ostr, std::integral_constant<int,isnumeric>(), l);
1099 }
1100 }
1101
1102 // Versions for writing matrices
1103 template<typename M>
1104 void writeMatrixMarket(const M& matrix,
1105 std::ostream& ostr,
1106 const std::integral_constant<int,1>&)
1107 {
1108 ostr<<matrix.N()*MatrixMarketImpl::mm_multipliers<M>::rows<<" "
1109 <<matrix.M()*MatrixMarketImpl::mm_multipliers<M>::cols<<" "
1110 <<countNonZeros(matrix)<<std::endl;
1111
1112 typedef typename M::const_iterator riterator;
1113 typedef typename M::ConstColIterator citerator;
1114 for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1115 for(citerator col = row->begin(); col != row->end(); ++col)
1116 // Matrix Market indexing start with 1!
1117 mm_print_entry(*col, row.index()*MatrixMarketImpl::mm_multipliers<M>::rows+1,
1118 col.index()*MatrixMarketImpl::mm_multipliers<M>::cols+1, ostr);
1119 }
1120
1121
1125 template<typename M>
1126 void writeMatrixMarket(const M& matrix,
1127 std::ostream& ostr)
1128 {
1129 using namespace MatrixMarketImpl;
1130
1131 // Write header information
1132 mm_header_printer<M>::print(ostr);
1133 mm_block_structure_header<M>::print(ostr,matrix);
1134 // Choose the correct function for matrix and vector
1135 writeMatrixMarket(matrix,ostr,std::integral_constant<int,IsMatrix<M>::value>());
1136 }
1137
1138 static const int default_precision = -1;
1150 template<typename M>
1151 void storeMatrixMarket(const M& matrix,
1152 std::string filename,
1153 int prec=default_precision)
1154 {
1155 auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1156 std::string rfilename;
1157 std::ofstream file;
1158 if (extension != "") {
1159 rfilename = pureFilename + extension;
1160 file.open(rfilename.c_str());
1161 if(!file)
1162 DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1163 }
1164 else {
1165 // only try .mm so we do not ignore potential errors
1166 rfilename = pureFilename + ".mm";
1167 file.open(rfilename.c_str());
1168 if(!file)
1169 DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1170 }
1171
1172 file.setf(std::ios::scientific,std::ios::floatfield);
1173 if(prec>0)
1174 file.precision(prec);
1175 writeMatrixMarket(matrix, file);
1176 file.close();
1177 }
1178
1179#if HAVE_MPI
1194 template<typename M, typename G, typename L>
1195 void storeMatrixMarket(const M& matrix,
1196 std::string filename,
1198 bool storeIndices=true,
1199 int prec=default_precision)
1200 {
1201 // Get our rank
1202 int rank = comm.communicator().rank();
1203 // Write the local matrix
1204 auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1205 std::string rfilename;
1206 std::ofstream file;
1207 if (extension != "") {
1208 rfilename = pureFilename + "_" + std::to_string(rank) + extension;
1209 file.open(rfilename.c_str());
1210 dverb<< rfilename <<std::endl;
1211 if(!file)
1212 DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1213 }
1214 else {
1215 // only try .mm so we do not ignore potential errors
1216 rfilename = pureFilename + "_" + std::to_string(rank) + ".mm";
1217 file.open(rfilename.c_str());
1218 dverb<< rfilename <<std::endl;
1219 if(!file)
1220 DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1221 }
1222 file.setf(std::ios::scientific,std::ios::floatfield);
1223 if(prec>0)
1224 file.precision(prec);
1225 writeMatrixMarket(matrix, file);
1226 file.close();
1227
1228 if(!storeIndices)
1229 return;
1230
1231 // Write the global to local index mapping
1232 rfilename = pureFilename + "_" + std::to_string(rank) + ".idx";
1233 file.open(rfilename.c_str());
1234 if(!file)
1235 DUNE_THROW(IOError, "Could not open file for storage: " << rfilename.c_str());
1236 file.setf(std::ios::scientific,std::ios::floatfield);
1238 typedef typename IndexSet::const_iterator Iterator;
1239 for(Iterator iter = comm.indexSet().begin();
1240 iter != comm.indexSet().end(); ++iter) {
1241 file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
1242 <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
1243 }
1244 // Store neighbour information for efficient remote indices setup.
1245 file<<"neighbours:";
1246 const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1247 typedef std::set<int>::const_iterator SIter;
1248 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1249 file<<" "<< *neighbour;
1250 }
1251 file.close();
1252 }
1253
1268 template<typename M, typename G, typename L>
1269 void loadMatrixMarket(M& matrix,
1270 const std::string& filename,
1272 bool readIndices=true)
1273 {
1274 using namespace MatrixMarketImpl;
1275
1277 typedef typename LocalIndexT::Attribute Attribute;
1278 // Get our rank
1279 int rank = comm.communicator().rank();
1280 // load local matrix
1281 auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1282 std::string rfilename;
1283 std::ifstream file;
1284 if (extension != "") {
1285 rfilename = pureFilename + "_" + std::to_string(rank) + extension;
1286 file.open(rfilename.c_str(), std::ios::in);
1287 dverb<< rfilename <<std::endl;
1288 if(!file)
1289 DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1290 }
1291 else {
1292 // try both .mm and .mtx
1293 rfilename = pureFilename + "_" + std::to_string(rank) + ".mm";
1294 file.open(rfilename.c_str(), std::ios::in);
1295 if(!file) {
1296 rfilename = pureFilename + "_" + std::to_string(rank) + ".mtx";
1297 file.open(rfilename.c_str(), std::ios::in);
1298 dverb<< rfilename <<std::endl;
1299 if(!file)
1300 DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1301 }
1302 }
1303 readMatrixMarket(matrix,file);
1304 file.close();
1305
1306 if(!readIndices)
1307 return;
1308
1309 // read indices
1311 IndexSet& pis=comm.pis;
1312 rfilename = pureFilename + "_" + std::to_string(rank) + ".idx";
1313 file.open(rfilename.c_str());
1314 if(!file)
1315 DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1316 if(pis.size()!=0)
1317 DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
1318
1319 pis.beginResize();
1320 while(!file.eof() && file.peek()!='n') {
1321 G g;
1322 file >>g;
1323 std::size_t l;
1324 file >>l;
1325 int c;
1326 file >>c;
1327 bool b;
1328 file >> b;
1329 pis.add(g,LocalIndexT(l,Attribute(c),b));
1330 lineFeed(file);
1331 }
1332 pis.endResize();
1333 if(!file.eof()) {
1334 // read neighbours
1335 std::string s;
1336 file>>s;
1337 if(s!="neighbours:")
1338 DUNE_THROW(MatrixMarketFormatError, "was expecting the string: \"neighbours:\"");
1339 std::set<int> nb;
1340 while(!file.eof()) {
1341 int i;
1342 file >> i;
1343 nb.insert(i);
1344 }
1345 file.close();
1346 comm.ri.setNeighbours(nb);
1347 }
1348 comm.ri.template rebuild<false>();
1349 }
1350
1351 #endif
1352
1363 template<typename M>
1364 void loadMatrixMarket(M& matrix,
1365 const std::string& filename)
1366 {
1367 auto [pureFilename, extension] = MatrixMarketImpl::splitFilename(filename);
1368 std::string rfilename;
1369 std::ifstream file;
1370 if (extension != "") {
1371 rfilename = pureFilename + extension;
1372 file.open(rfilename.c_str());
1373 if(!file)
1374 DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1375 }
1376 else {
1377 // try both .mm and .mtx
1378 rfilename = pureFilename + ".mm";
1379 file.open(rfilename.c_str(), std::ios::in);
1380 if(!file) {
1381 rfilename = pureFilename + ".mtx";
1382 file.open(rfilename.c_str(), std::ios::in);
1383 if(!file)
1384 DUNE_THROW(IOError, "Could not open file: " << rfilename.c_str());
1385 }
1386 }
1387 readMatrixMarket(matrix,file);
1388 file.close();
1389 }
1390
1392}
1393#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:466
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:671
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:677
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1100
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:2007
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1094
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:2001
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:830
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:858
A vector of blocks with memory management.
Definition: bvector.hh:392
void resize(size_type size)
Resize the vector.
Definition: bvector.hh:496
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bvector.hh:398
int rank() const
Return rank, is between 0 and size()-1.
Definition: mpicommunication.hh:133
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception class for I/O errors.
Definition: exceptions.hh:231
Index Set Interface base class.
Definition: indexidset.hh:78
auto size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:223
Exception indicating that the index set is not in the expected state.
Definition: indexset.hh:205
A generic dynamic dense matrix.
Definition: matrix.hh:561
Default exception for dummy implementations.
Definition: exceptions.hh:263
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:174
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:462
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:471
A few common exception classes.
Classes providing communication interfaces for overlapping Schwarz methods.
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:239
Stream & operator>>(Stream &stream, std::tuple< Ts... > &t)
Read a std::tuple.
Definition: streamoperators.hh:43
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
void readMatrixMarket(Dune::BlockVector< T, A > &vector, std::istream &istr)
Reads a BlockVector from a matrix market file.
Definition: matrixmarket.hh:930
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:1151
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:1269
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:119
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:638
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition: interface.hh:324
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:117
Some handy generic functions for ISTL matrices.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
T lane(std::size_t l, const T &v)
access a lane of a simd vector (scalar version)
Definition: simd.hh:366
STL namespace.
Include file for users of the SIMD abstraction layer.
Standard Dune debug streams.
Test whether a type is an ISTL Matrix.
Definition: matrixutils.hh:504
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:677
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< T > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:684
void operator()(const std::vector< std::set< IndexData< D > > > &rows, BCRSMatrix< FieldMatrix< T, brows, bcols > > &matrix)
Sets the matrix values.
Definition: matrixmarket.hh:702
a wrapper class of numeric values.
Definition: matrixmarket.hh:595
Utility class for marking the pattern type of the MatrixMarket matrices.
Definition: matrixmarket.hh:607
Metaprogram for writing the ISTL block structure header.
Definition: matrixmarket.hh:225
Meta program to write the correct Matrix Market header.
Definition: matrixmarket.hh:174
Helper metaprogram to get the matrix market string representation of the numeric type.
Definition: matrixmarket.hh:76
@ is_numeric
Whether T is a supported numeric type.
Definition: matrixmarket.hh:81
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)