dune-istl  2.8.0
io.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_IO_HH
4 #define DUNE_ISTL_IO_HH
5 
6 #include <cmath>
7 #include <complex>
8 #include <limits>
9 #include <ios>
10 #include <iomanip>
11 #include <fstream>
12 #include <string>
13 
14 #include "matrixutils.hh"
15 #include "istlexception.hh"
16 #include <dune/common/fvector.hh>
17 #include <dune/common/fmatrix.hh>
18 #include <dune/common/hybridutilities.hh>
19 
20 #include <dune/istl/bcrsmatrix.hh>
21 
22 
23 namespace Dune {
24 
37  //
38  // pretty printing of vectors
39  //
40 
48  template<class V>
49  void recursive_printvector (std::ostream& s, const V& v, std::string rowtext,
50  int& counter, int columns, int width)
51  {
52  if constexpr (IsNumber<V>())
53  {
54  // Print one number
55  if (counter%columns==0)
56  {
57  s << rowtext; // start a new row
58  s << " "; // space in front of each entry
59  s.width(4); // set width for counter
60  s << counter; // number of first entry in a line
61  }
62  s << " "; // space in front of each entry
63  s.width(width); // set width for each entry anew
64  s << v; // yeah, the number !
65  counter++; // increment the counter
66  if (counter%columns==0)
67  s << std::endl; // start a new line
68  }
69  else
70  {
71  // Recursively print a vector
72  for (const auto& entry : v)
73  recursive_printvector(s,entry,rowtext,counter,columns,width);
74  }
75  }
76 
77 
85  template<class V>
86  void printvector (std::ostream& s, const V& v, std::string title,
87  std::string rowtext, int columns=1, int width=10,
88  int precision=2)
89  {
90  // count the numbers printed to make columns
91  int counter=0;
92 
93  // remember old flags
94  std::ios_base::fmtflags oldflags = s.flags();
95 
96  // set the output format
97  s.setf(std::ios_base::scientific, std::ios_base::floatfield);
98  int oldprec = s.precision();
99  s.precision(precision);
100 
101  // print title
102  s << title << " [blocks=" << v.N() << ",dimension=" << v.dim() << "]"
103  << std::endl;
104 
105  // print data from all blocks
106  recursive_printvector(s,v,rowtext,counter,columns,width);
107 
108  // check if new line is required
109  if (counter%columns!=0)
110  s << std::endl;
111 
112  // reset the output format
113  s.flags(oldflags);
114  s.precision(oldprec);
115  }
116 
117 
119  //
120  // pretty printing of matrices
121  //
122 
130  inline void fill_row (std::ostream& s, int m, int width, [[maybe_unused]] int precision)
131  {
132  for (int j=0; j<m; j++)
133  {
134  s << " "; // space in front of each entry
135  s.width(width); // set width for each entry anew
136  s << "."; // yeah, the number !
137  }
138  }
139 
147  template<class K>
148  void print_row (std::ostream& s, const K& value,
149  [[maybe_unused]] typename FieldMatrix<K,1,1>::size_type I,
150  [[maybe_unused]] typename FieldMatrix<K,1,1>::size_type J,
151  [[maybe_unused]] typename FieldMatrix<K,1,1>::size_type therow,
152  int width,
153  [[maybe_unused]] int precision,
154  typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
155  {
156  s << " "; // space in front of each entry
157  s.width(width); // set width for each entry anew
158  s << value;
159  }
160 
168  template<class M>
169  void print_row (std::ostream& s, const M& A, typename M::size_type I,
170  typename M::size_type J, typename M::size_type therow,
171  int width, int precision,
172  typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
173  {
174  typename M::size_type i0=I;
175  for (typename M::size_type i=0; i<A.N(); i++)
176  {
177  if (therow>=i0 && therow<i0+MatrixDimension<M>::rowdim(A,i))
178  {
179  // the row is in this block row !
180  typename M::size_type j0=J;
181  for (typename M::size_type j=0; j<A.M(); j++)
182  {
183  // find this block
184  typename M::ConstColIterator it = A[i].find(j);
185 
186  // print row or filler
187  if (it!=A[i].end())
188  print_row(s,*it,i0,j0,therow,width,precision);
189  else
190  fill_row(s,MatrixDimension<M>::coldim(A,j),width,precision);
191 
192  // advance columns
193  j0 += MatrixDimension<M>::coldim(A,j);
194  }
195  }
196  // advance rows
197  i0 += MatrixDimension<M>::rowdim(A,i);
198  }
199  }
200 
209  template<class M>
210  void printmatrix (std::ostream& s, const M& A, std::string title,
211  std::string rowtext, int width=10, int precision=2)
212  {
213 
214  // remember old flags
215  std::ios_base::fmtflags oldflags = s.flags();
216 
217  // set the output format
218  s.setf(std::ios_base::scientific, std::ios_base::floatfield);
219  int oldprec = s.precision();
220  s.precision(precision);
221 
222  // print title
223  s << title
224  << " [n=" << A.N()
225  << ",m=" << A.M()
226  << ",rowdim=" << MatrixDimension<M>::rowdim(A)
227  << ",coldim=" << MatrixDimension<M>::coldim(A)
228  << "]" << std::endl;
229 
230  // print all rows
231  for (typename M::size_type i=0; i<MatrixDimension<M>::rowdim(A); i++)
232  {
233  s << rowtext; // start a new row
234  s << " "; // space in front of each entry
235  s.width(4); // set width for counter
236  s << i; // number of first entry in a line
237  print_row(s,A,0,0,i,width,precision); // generic print
238  s << std::endl; // start a new line
239  }
240 
241  // reset the output format
242  s.flags(oldflags);
243  s.precision(oldprec);
244  }
245 
267  template<class B, int n, int m, class A>
268  void printSparseMatrix(std::ostream& s,
270  std::string title, std::string rowtext,
271  int width=3, int precision=2)
272  {
274  // remember old flags
275  std::ios_base::fmtflags oldflags = s.flags();
276  // set the output format
277  s.setf(std::ios_base::scientific, std::ios_base::floatfield);
278  int oldprec = s.precision();
279  s.precision(precision);
280  // print title
281  s << title
282  << " [n=" << mat.N()
283  << ",m=" << mat.M()
284  << ",rowdim=" << MatrixDimension<Matrix>::rowdim(mat)
285  << ",coldim=" << MatrixDimension<Matrix>::coldim(mat)
286  << "]" << std::endl;
287 
288  typedef typename Matrix::ConstRowIterator Row;
289 
290  for(Row row=mat.begin(); row != mat.end(); ++row) {
291  int skipcols=0;
292  bool reachedEnd=false;
293 
294  while(!reachedEnd) {
295  for(int innerrow=0; innerrow<n; ++innerrow) {
296  int count=0;
297  typedef typename Matrix::ConstColIterator Col;
298  Col col=row->begin();
299  for(; col != row->end(); ++col,++count) {
300  if(count<skipcols)
301  continue;
302  if(count>=skipcols+width)
303  break;
304  if(innerrow==0) {
305  if(count==skipcols) {
306  s << rowtext; // start a new row
307  s << " "; // space in front of each entry
308  s.width(4); // set width for counter
309  s << row.index()<<": "; // number of first entry in a line
310  }
311  s.width(4);
312  s<<col.index()<<": |";
313  } else {
314  if(count==skipcols) {
315  for(typename std::string::size_type i=0; i < rowtext.length(); i++)
316  s<<" ";
317  s<<" ";
318  }
319  s<<" |";
320  }
321  for(int innercol=0; innercol < m; ++innercol) {
322  s.width(9);
323  s<<(*col)[innerrow][innercol]<<" ";
324  }
325 
326  s<<"|";
327  }
328  if(innerrow==n-1 && col==row->end())
329  reachedEnd = true;
330  else
331  s << std::endl;
332  }
333  skipcols += width;
334  s << std::endl;
335  }
336  s << std::endl;
337  }
338 
339  // reset the output format
340  s.flags(oldflags);
341  s.precision(oldprec);
342  }
343 
344  namespace
345  {
346  template<typename T>
347  struct MatlabPODWriter
348  {
349  static std::ostream& write(const T& t, std::ostream& s)
350  {
351  s << t;
352  return s;
353  }
354  };
355  template<typename T>
356  struct MatlabPODWriter<std::complex<T> >
357  {
358  static std::ostream& write(const std::complex<T>& t, std::ostream& s)
359  {
360  s << t.real() << " " << t.imag();
361  return s;
362  }
363  };
364  } // anonymous namespace
365 
375  template <class FieldType>
376  void writeMatrixToMatlabHelper(const FieldType& value,
377  int rowOffset, int colOffset,
378  std::ostream& s,
379  typename std::enable_if_t<Dune::IsNumber<FieldType>::value>* sfinae = nullptr)
380  {
381  //+1 for Matlab numbering
382  s << rowOffset + 1 << " " << colOffset + 1 << " ";
383  MatlabPODWriter<FieldType>::write(value, s)<< std::endl;
384  }
385 
393  template <class MatrixType>
394  void writeMatrixToMatlabHelper(const MatrixType& matrix,
395  int externalRowOffset, int externalColOffset,
396  std::ostream& s,
397  typename std::enable_if_t<!Dune::IsNumber<MatrixType>::value>* sfinae = nullptr)
398  {
399  // Precompute the accumulated sizes of the columns
400  std::vector<typename MatrixType::size_type> colOffset(matrix.M());
401  if (colOffset.size() > 0)
402  colOffset[0] = 0;
403 
404  for (typename MatrixType::size_type i=0; i<matrix.M()-1; i++)
405  colOffset[i+1] = colOffset[i] +
407 
408  typename MatrixType::size_type rowOffset = 0;
409 
410  // Loop over all matrix rows
411  for (typename MatrixType::size_type rowIdx=0; rowIdx<matrix.N(); rowIdx++)
412  {
413  auto cIt = matrix[rowIdx].begin();
414  auto cEndIt = matrix[rowIdx].end();
415 
416  // Loop over all columns in this row
417  for (; cIt!=cEndIt; ++cIt)
419  externalRowOffset+rowOffset,
420  externalColOffset + colOffset[cIt.index()],
421  s);
422 
423  rowOffset += MatrixDimension<MatrixType>::rowdim(matrix, rowIdx);
424  }
425 
426  }
427 
447  template <class MatrixType>
448  void writeMatrixToMatlab(const MatrixType& matrix,
449  const std::string& filename, int outputPrecision = 18)
450  {
451  std::ofstream outStream(filename.c_str());
452  int oldPrecision = outStream.precision();
453  outStream.precision(outputPrecision);
454 
455  writeMatrixToMatlabHelper(matrix, 0, 0, outStream);
456  outStream.precision(oldPrecision);
457  }
458 
459  // Recursively write vector entries to a stream
460  template<class V>
461  void writeVectorToMatlabHelper (const V& v, std::ostream& stream)
462  {
463  if constexpr (IsNumber<V>()) {
464  stream << v << std::endl;
465  } else {
466  for (const auto& entry : v)
467  writeVectorToMatlabHelper(entry, stream);
468  }
469  }
470 
488  template <class VectorType>
489  void writeVectorToMatlab(const VectorType& vector,
490  const std::string& filename, int outputPrecision = 18)
491  {
492  std::ofstream outStream(filename.c_str());
493  int oldPrecision = outStream.precision();
494  outStream.precision(outputPrecision);
495 
496  writeVectorToMatlabHelper(vector, outStream);
497  outStream.precision(oldPrecision);
498  }
499 
502 } // namespace Dune
503 
504 #endif
Implementation of the BCRSMatrix class.
Some handy generic functions for ISTL matrices.
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
void fill_row(std::ostream &s, int m, int width, [[maybe_unused]] int precision)
Print a row of zeros for a non-existing block.
Definition: io.hh:130
void writeMatrixToMatlab(const MatrixType &matrix, const std::string &filename, int outputPrecision=18)
Writes sparse matrix in a Matlab-readable format.
Definition: io.hh:448
void print_row(std::ostream &s, const K &value, [[maybe_unused]] typename FieldMatrix< K, 1, 1 >::size_type I, [[maybe_unused]] typename FieldMatrix< K, 1, 1 >::size_type J, [[maybe_unused]] typename FieldMatrix< K, 1, 1 >::size_type therow, int width, [[maybe_unused]] int precision, typename std::enable_if_t< Dune::IsNumber< K >::value > *sfinae=nullptr)
Print one row of a matrix, specialization for number types.
Definition: io.hh:148
void printmatrix(std::ostream &s, const M &A, std::string title, std::string rowtext, int width=10, int precision=2)
Print a generic block matrix.
Definition: io.hh:210
void printvector(std::ostream &s, const V &v, std::string title, std::string rowtext, int columns=1, int width=10, int precision=2)
Print an ISTL vector.
Definition: io.hh:86
void writeMatrixToMatlabHelper(const FieldType &value, int rowOffset, int colOffset, std::ostream &s, typename std::enable_if_t< Dune::IsNumber< FieldType >::value > *sfinae=nullptr)
Helper method for the writeMatrixToMatlab routine.
Definition: io.hh:376
void writeVectorToMatlabHelper(const V &v, std::ostream &stream)
Definition: io.hh:461
void writeVectorToMatlab(const VectorType &vector, const std::string &filename, int outputPrecision=18)
Writes vectors in a Matlab-readable format.
Definition: io.hh:489
void recursive_printvector(std::ostream &s, const V &v, std::string rowtext, int &counter, int columns, int width)
Recursively print a vector.
Definition: io.hh:49
void printSparseMatrix(std::ostream &s, const BCRSMatrix< FieldMatrix< B, n, m >, A > &mat, std::string title, std::string rowtext, int width=3, int precision=2)
Prints a BCRSMatrix with fixed sized blocks.
Definition: io.hh:268
Definition: allocator.hh:9
Definition: matrixutils.hh:209
static auto coldim(const M &A)
Definition: matrixutils.hh:217
static auto rowdim(const M &A)
Definition: matrixutils.hh:212
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:618
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:612
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
Definition: matrixutils.hh:25