3 #ifndef DUNE_ISTL_MATRIXMARKET_HH
4 #define DUNE_ISTL_MATRIXMARKET_HH
19 #include <type_traits>
22 #include <dune/common/exceptions.hh>
23 #include <dune/common/fmatrix.hh>
24 #include <dune/common/fvector.hh>
25 #include <dune/common/unused.hh>
26 #include <dune/common/hybridutilities.hh>
27 #include <dune/common/stdstreams.hh>
62 namespace MatrixMarketImpl
93 static std::string
str()
174 template<
typename T,
typename A>
179 os<<
"%%MatrixMarket matrix coordinate ";
180 os<<mm_numeric_type<typename Imp::BlockTraits<T>::field_type>::str()<<
" general"<<std::endl;
184 template<
typename B,
typename A>
189 os<<
"%%MatrixMarket matrix array ";
190 os<<mm_numeric_type<typename Imp::BlockTraits<B>::field_type>::str()<<
" general"<<std::endl;
194 template<
typename T,
int j>
199 os<<
"%%MatrixMarket matrix array ";
200 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
204 template<
typename T,
int i,
int j>
209 os<<
"%%MatrixMarket matrix array ";
210 os<<mm_numeric_type<T>::str()<<
" general"<<std::endl;
225 template<
typename T,
typename A>
229 static_assert(IsNumber<T>::value,
"Only scalar entries are expected here!");
231 static void print(std::ostream& os,
const M&)
233 os<<
"% ISTL_STRUCT blocked ";
234 os<<
"1 1"<<std::endl;
238 template<
typename T,
typename A,
int i>
243 static void print(std::ostream& os,
const M&)
245 os<<
"% ISTL_STRUCT blocked ";
246 os<<i<<
" "<<1<<std::endl;
250 template<
typename T,
typename A>
254 static_assert(IsNumber<T>::value,
"Only scalar entries are expected here!");
256 static void print(std::ostream& os,
const M&)
258 os<<
"% ISTL_STRUCT blocked ";
259 os<<
"1 1"<<std::endl;
263 template<
typename T,
typename A,
int i,
int j>
268 static void print(std::ostream& os,
const M&)
270 os<<
"% ISTL_STRUCT blocked ";
271 os<<i<<
" "<<j<<std::endl;
276 template<
typename T,
int i,
int j>
281 static void print(std::ostream& os,
const M& m)
285 template<
typename T,
int i>
288 typedef FieldVector<T,i>
M;
290 static void print(std::ostream& os,
const M& m)
345 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
360 dverb<<buffer<<std::endl;
362 if(buffer!=
"%%MatrixMarket") {
364 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
375 if(buffer !=
"matrix")
378 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
392 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
399 if(buffer !=
"array")
401 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
408 if(buffer !=
"coordinate")
410 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
416 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
430 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
436 if(buffer !=
"integer")
438 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
447 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
454 if(buffer !=
"complex")
456 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
463 if(buffer !=
"pattern")
465 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
471 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
480 std::transform(buffer.begin(), buffer.end(), buffer.begin(),
486 if(buffer !=
"general")
488 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
495 if(buffer !=
"hermitian")
497 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
503 if(buffer.size()==1) {
504 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
512 if(buffer !=
"symmetric")
514 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
521 if(buffer !=
"skew-symmetric")
523 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
529 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
534 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
537 file.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
543 template<std::
size_t brows, std::
size_t bcols>
544 std::tuple<std::size_t, std::size_t, std::size_t>
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;
555 blockentries = entries/blocksize;
break;
557 blockentries = 2*entries/blocksize;
break;
559 blockentries = (2*entries-rows)/blocksize;
break;
561 blockentries = (2*entries-rows)/blocksize;
break;
563 throw Dune::NotImplemented();
565 return std::make_tuple(blockrows, blockcols, blockentries);
619 DUNE_UNUSED_PARAMETER(num);
645 return is>>data.number;
664 data.
number = {real.number, imag.number};
674 template<
typename D,
int brows,
int bcols>
686 static_assert(IsNumber<T>::value && brows==1 && bcols==1,
"Only scalar entries are expected here!");
687 for (
auto iter=matrix.
begin(); iter!= matrix.
end(); ++iter)
689 auto brow=iter.index();
690 for (
auto siter=rows[brow].begin(); siter != rows[brow].end(); ++siter)
691 (*iter)[siter->index] = siter->number;
704 for (
auto iter=matrix.begin(); iter!= matrix.end(); ++iter)
706 for (
auto brow=iter.index()*brows,
707 browend=iter.index()*brows+brows;
708 brow<browend; ++brow)
710 for (
auto siter=rows[brow].begin(), send=rows[brow].end();
711 siter != send; ++siter)
712 (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
718 template<
int brows,
int bcols>
728 template<
class T>
struct is_complex<std::complex<T>> : std::true_type {};
732 std::enable_if_t<!is_complex<T>::value, T>
conj(
const T& r){
737 std::enable_if_t<is_complex<T>::value, T>
conj(
const T& r){
745 template<
typename B,
typename A>
754 template<
typename B,
int i,
int j,
typename A>
763 template<
typename T,
typename A,
typename D>
765 std::istream& file, std::size_t entries,
777 std::vector<std::set<IndexData<D> > > rows(matrix.
N()*brows);
779 auto readloop = [&] (
auto symmetryFixup) {
780 for(std::size_t i = 0; i < entries; ++i) {
786 assert(row/bcols<matrix.
N());
788 assert(data.index/bcols<matrix.
M());
789 rows[row].insert(data);
791 symmetryFixup(row, data);
798 readloop([](
auto...){});
801 readloop([&](
auto row,
auto data) {
803 data_sym.
index = row;
804 rows[data.index].insert(data_sym);
808 readloop([&](
auto row,
auto data) {
810 data_sym.number = -data.number;
811 data_sym.
index = row;
812 rows[data.index].insert(data_sym);
816 readloop([&](
auto row,
auto data) {
818 data_sym.number =
conj(data.number);
819 data_sym.
index = row;
820 rows[data.index].insert(data_sym);
824 DUNE_THROW(Dune::NotImplemented,
825 "Only general, symmetric, skew-symmetric and hermitian is supported right now!");
830 for(
typename Matrix::CreateIterator iter=matrix.
createbegin();
833 for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
834 brow<browend; ++brow)
836 typedef typename std::set<IndexData<D> >::const_iterator Siter;
837 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
838 siter != send; ++siter, ++nnz)
839 iter.insert(siter->index/bcols);
848 Setter(rows, matrix);
860 using namespace MatrixMarketImpl;
863 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
864 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
867 istr.seekg(0, std::ios::beg);
884 template<
typename T,
typename A>
889 for (
int i=0; size>0; ++i, --size)
893 template<
typename T,
typename A,
int entries>
898 for(
int i=0; size>0; ++i, --size) {
901 vector[i/entries][i%entries]=val;
912 template<
typename T,
typename A>
916 using namespace MatrixMarketImpl;
919 std::size_t rows, cols;
928 Dune::Hybrid::ifElse(Dune::IsNumber<T>(),
934 auto blocksize = id(dummy).size();
935 std::size_t size=rows/blocksize;
936 if(size*blocksize!=rows)
942 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
952 template<
typename T,
typename A>
956 using namespace MatrixMarketImpl;
961 std::cerr <<
"First line was not a correct Matrix Market banner. Using default:\n"
962 <<
"%%MatrixMarket matrix coordinate real general"<<std::endl;
965 istr.seekg(0, std::ios::beg);
969 std::size_t rows, cols, entries;
985 std::size_t nnz, blockrows, blockcols;
988 constexpr
int brows = mm_multipliers<Matrix>::rows;
989 constexpr
int bcols = mm_multipliers<Matrix>::cols;
991 std::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
993 istr.ignore(std::numeric_limits<std::streamsize>::max(),
'\n');
996 matrix.
setSize(blockrows, blockcols);
1000 DUNE_THROW(Dune::NotImplemented,
"Array format currently not supported for matrices!");
1002 readSparseEntries(matrix, istr, entries, header, NumericWrapper<typename Matrix::field_type>());
1006 template<
typename B>
1012 Hybrid::ifElse(IsNumber<B>(),
1014 ostr << rowidx <<
" " << colidx <<
" " << entry << std::endl;
1017 for (
auto row=
id(entry).begin(); row != id(entry).end(); ++row, ++rowidx) {
1019 for (
auto col = row->begin();
col != row->end(); ++
col, ++coli)
1020 ostr<< rowidx<<
" "<<coli<<
" "<<*
col<<std::endl;
1026 template<
typename V>
1028 const std::integral_constant<int,1>&)
1030 ostr<<entry<<std::endl;
1034 template<
typename V>
1036 const std::integral_constant<int,0>&)
1038 using namespace MatrixMarketImpl;
1042 typedef typename V::const_iterator VIter;
1044 for(VIter i=vector.begin(); i != vector.end(); ++i)
1047 std::integral_constant<int,isnumeric>());
1050 template<
typename T,
typename A>
1053 return vector.size();
1056 template<
typename T,
typename A,
int i>
1059 return vector.size()*i;
1063 template<
typename V>
1065 const std::integral_constant<int,0>&)
1067 using namespace MatrixMarketImpl;
1075 template<
typename M>
1078 const std::integral_constant<int,1>&)
1084 typedef typename M::const_iterator riterator;
1085 typedef typename M::ConstColIterator citerator;
1086 for(riterator row=matrix.begin(); row != matrix.end(); ++row)
1087 for(citerator
col = row->begin();
col != row->end(); ++
col)
1097 template<
typename M>
1101 using namespace MatrixMarketImpl;
1104 mm_header_printer<M>::print(ostr);
1105 mm_block_structure_header<M>::print(ostr,matrix);
1121 template<
typename M>
1123 std::string filename)
1125 std::ofstream file(filename.c_str());
1126 file.setf(std::ios::scientific,std::ios::floatfield);
1145 template<
typename M,
typename G,
typename L>
1147 std::string filename,
1148 const OwnerOverlapCopyCommunication<G,L>& comm,
1149 bool storeIndices=
true)
1152 int rank = comm.communicator().rank();
1154 std::ostringstream rfilename;
1155 rfilename<<filename <<
"_"<<rank<<
".mm";
1156 dverb<< rfilename.str()<<std::endl;
1157 std::ofstream file(rfilename.str().c_str());
1158 file.setf(std::ios::scientific,std::ios::floatfield);
1167 rfilename<<filename<<
"_"<<rank<<
".idx";
1168 file.open(rfilename.str().c_str());
1169 file.setf(std::ios::scientific,std::ios::floatfield);
1170 typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1171 typedef typename IndexSet::const_iterator Iterator;
1172 for(Iterator iter = comm.indexSet().begin();
1173 iter != comm.indexSet().end(); ++iter) {
1174 file << iter->global()<<
" "<<(std::size_t)iter->local()<<
" "
1175 <<(int)iter->local().attribute()<<
" "<<(int)iter->local().isPublic()<<std::endl;
1178 file<<
"neighbours:";
1179 const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
1180 typedef std::set<int>::const_iterator SIter;
1181 for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour) {
1182 file<<
" "<< *neighbour;
1201 template<
typename M,
typename G,
typename L>
1203 const std::string& filename,
1204 OwnerOverlapCopyCommunication<G,L>& comm,
1205 bool readIndices=
true)
1207 using namespace MatrixMarketImpl;
1209 typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet::LocalIndex LocalIndex;
1210 typedef typename LocalIndex::Attribute Attribute;
1212 int rank = comm.communicator().rank();
1214 std::ostringstream rfilename;
1215 rfilename<<filename <<
"_"<<rank<<
".mm";
1217 file.open(rfilename.str().c_str(), std::ios::in);
1219 DUNE_THROW(IOError,
"Could not open file: " << rfilename.str().c_str());
1229 typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
1230 IndexSet& pis=comm.pis;
1232 rfilename<<filename<<
"_"<<rank<<
".idx";
1233 file.open(rfilename.str().c_str());
1235 DUNE_THROW(InvalidIndexSetState,
"Index set is not empty!");
1238 while(!file.eof() && file.peek()!=
'n') {
1247 pis.add(g,LocalIndex(l,Attribute(c),b));
1255 if(s!=
"neighbours:")
1256 DUNE_THROW(MatrixMarketFormatError,
"was expecting the string: \"neighbours:\"");
1258 while(!file.eof()) {
1264 comm.ri.setNeighbours(nb);
1266 comm.ri.template rebuild<false>();
1281 template<
typename M>
1283 const std::string& filename)
1286 file.open(filename.c_str(), std::ios::in);
1288 DUNE_THROW(IOError,
"Could not open file: " << filename);