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 25 /** 26 @addtogroup ISTL_IO 27 @{ 28 */ 29 30 31 /** \file 32 33 \brief Some generic functions for pretty printing vectors and matrices 34 */ 35 36 //////////////////////////////////////////////////////////////////////// 37 // 38 // pretty printing of vectors 39 // 40 41 /** 42 * \brief Recursively print a vector 43 * 44 * \code 45 * #include <dune/istl/io.hh> 46 * \endcode 47 */ 48 template<class V> recursive_printvector(std::ostream & s,const V & v,std::string rowtext,int & counter,int columns,int width)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 78 /** 79 * \brief Print an ISTL vector 80 * 81 * \code 82 * #include <dune/istl/io.hh> 83 * \endcode 84 */ 85 template<class V> printvector(std::ostream & s,const V & v,std::string title,std::string rowtext,int columns=1,int width=10,int precision=2)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 118 //////////////////////////////////////////////////////////////////////// 119 // 120 // pretty printing of matrices 121 // 122 123 /** 124 * \brief Print a row of zeros for a non-existing block 125 * 126 * \code 127 * #include <dune/istl/io.hh> 128 * \endcode 129 */ fill_row(std::ostream & s,int m,int width,int precision)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 140 /** 141 * \brief Print one row of a matrix, specialization for number types 142 * 143 * \code 144 * #include <dune/istl/io.hh> 145 * \endcode 146 */ 147 template<class K> print_row(std::ostream & s,const K & value,typename FieldMatrix<K,1,1>::size_type I,typename FieldMatrix<K,1,1>::size_type J,typename FieldMatrix<K,1,1>::size_type therow,int width,int precision,typename std::enable_if_t<Dune::IsNumber<K>::value> * sfinae=nullptr)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 161 /** 162 * \brief Print one row of a matrix 163 * 164 * \code 165 * #include <dune/istl/io.hh> 166 * \endcode 167 */ 168 template<class M> print_row(std::ostream & s,const M & A,typename M::size_type I,typename M::size_type J,typename M::size_type therow,int width,int precision,typename std::enable_if_t<!Dune::IsNumber<M>::value> * sfinae=nullptr)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 201 /** 202 * \brief Print a generic block matrix 203 * 204 * \code 205 * #include <dune/istl/io.hh> 206 * \endcode 207 * \bug Empty rows and columns are omitted by this method. (FlySpray #7) 208 */ 209 template<class M> printmatrix(std::ostream & s,const M & A,std::string title,std::string rowtext,int width=10,int precision=2)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 246 /** 247 * \brief Prints a BCRSMatrix with fixed sized blocks. 248 * 249 * \code 250 * #include <dune/istl/io.hh> 251 * \endcode 252 * 253 * Only the nonzero entries will be printed as matrix blocks 254 * together with their 255 * corresponding column index and all others will be omitted. 256 * 257 * This might be preferable over printmatrix in the case of big 258 * sparse matrices with nonscalar blocks. 259 * 260 * @param s The ostream to print to. 261 * @param mat The matrix to print. 262 * @param title The title for the matrix. 263 * @param rowtext The text to prepend to each print out of a matrix row. 264 * @param width The number of nonzero blocks to print in one line. 265 * @param precision The precision to use when printing the numbers. 266 */ 267 template<class B, int n, int m, class A> printSparseMatrix(std::ostream & s,const BCRSMatrix<FieldMatrix<B,n,m>,A> & mat,std::string title,std::string rowtext,int width=3,int precision=2)268 void printSparseMatrix(std::ostream& s, 269 const BCRSMatrix<FieldMatrix<B,n,m>,A>& mat, 270 std::string title, std::string rowtext, 271 int width=3, int precision=2) 272 { 273 typedef BCRSMatrix<FieldMatrix<B,n,m>,A> Matrix; 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 { writeDune::__anone0b50e8f0111::MatlabPODWriter349 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 { writeDune::__anone0b50e8f0111::MatlabPODWriter358 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 366 /** 367 * \brief Helper method for the writeMatrixToMatlab routine. 368 * 369 * \code 370 * #include <dune/istl/io.hh> 371 * \endcode 372 * 373 * This specialization for numbers ends the recursion 374 */ 375 template <class FieldType> writeMatrixToMatlabHelper(const FieldType & value,int rowOffset,int colOffset,std::ostream & s,typename std::enable_if_t<Dune::IsNumber<FieldType>::value> * sfinae=nullptr)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 386 /** 387 * \brief Helper method for the writeMatrixToMatlab routine. 388 * 389 * \code 390 * #include <dune/istl/io.hh> 391 * \endcode 392 */ 393 template <class MatrixType> writeMatrixToMatlabHelper(const MatrixType & matrix,int externalRowOffset,int externalColOffset,std::ostream & s,typename std::enable_if_t<!Dune::IsNumber<MatrixType>::value> * sfinae=nullptr)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] + 406 MatrixDimension<MatrixType>::coldim(matrix,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) 418 writeMatrixToMatlabHelper(*cIt, 419 externalRowOffset+rowOffset, 420 externalColOffset + colOffset[cIt.index()], 421 s); 422 423 rowOffset += MatrixDimension<MatrixType>::rowdim(matrix, rowIdx); 424 } 425 426 } 427 428 /** 429 * \brief Writes sparse matrix in a Matlab-readable format 430 * 431 * \code 432 * #include <dune/istl/io.hh> 433 * \endcode 434 * This routine writes the argument BCRSMatrix to a file with the name given 435 * by the filename argument. The file format is ASCII, with no header, and 436 * three data columns. Each row describes a scalar matrix entry and 437 * consists of the matrix row and column numbers (both counted starting from 438 * 1), and the matrix entry. Such a file can be read from Matlab using the 439 * command 440 * \code 441 * new_mat = spconvert(load('filename')); 442 * \endcode 443 * @param matrix reference to matrix 444 * @param filename 445 * @param outputPrecision (number of digits) which is used to write the output file 446 */ 447 template <class MatrixType> writeMatrixToMatlab(const MatrixType & matrix,const std::string & filename,int outputPrecision=18)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> writeVectorToMatlabHelper(const V & v,std::ostream & stream)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 471 /** 472 * \brief Writes vectors in a Matlab-readable format 473 * 474 * \code 475 * #include <dune/istl/io.hh> 476 * \endcode 477 * This routine writes the argument block vector to a file with the name given 478 * by the filename argument. The file format is ASCII, with no header, and 479 * a single data column. Such a file can be read from Matlab using the 480 * command 481 * \code 482 * new_vec = load('filename'); 483 * \endcode 484 * \param vector reference to vector to be printed to output file 485 * \param filename filename of output file 486 * \param outputPrecision (number of digits) which is used to write the output file 487 */ 488 template <class VectorType> writeVectorToMatlab(const VectorType & vector,const std::string & filename,int outputPrecision=18)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 500 /** @} end documentation */ 501 502 } // namespace Dune 503 504 #endif 505