dune-istl
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=4 sw=2 sts=2: 00003 #ifndef DUNE_ISTLIO_HH 00004 #define DUNE_ISTLIO_HH 00005 00006 #include <cmath> 00007 #include <complex> 00008 #include <limits> 00009 #include <ios> 00010 #include <iomanip> 00011 #include <fstream> 00012 #include <string> 00013 00014 #include "matrixutils.hh" 00015 #include "istlexception.hh" 00016 #include <dune/common/fvector.hh> 00017 #include <dune/common/fmatrix.hh> 00018 00019 #include <dune/istl/matrix.hh> 00020 #include "bcrsmatrix.hh" 00021 00022 00023 namespace Dune { 00024 00036 00037 // 00038 // pretty printing of vectors 00039 // 00040 00041 // recursively print all the blocks 00047 template<class V> 00048 void recursive_printvector (std::ostream& s, const V& v, std::string rowtext, 00049 int& counter, int columns, int width, 00050 int precision) 00051 { 00052 for (typename V::ConstIterator i=v.begin(); i!=v.end(); ++i) 00053 recursive_printvector(s,*i,rowtext,counter,columns,width,precision); 00054 } 00055 00056 // recursively print all the blocks -- specialization for FieldVector 00062 template<class K, int n> 00063 void recursive_printvector (std::ostream& s, const FieldVector<K,n>& v, 00064 std::string rowtext, int& counter, int columns, 00065 int width, int precision) 00066 { 00067 // we now can print n numbers 00068 for (int i=0; i<n; i++) 00069 { 00070 if (counter%columns==0) 00071 { 00072 s << rowtext; // start a new row 00073 s << " "; // space in front of each entry 00074 s.width(4); // set width for counter 00075 s << counter; // number of first entry in a line 00076 } 00077 s << " "; // space in front of each entry 00078 s.width(width); // set width for each entry anew 00079 s << v[i]; // yeah, the number ! 00080 counter++; // increment the counter 00081 if (counter%columns==0) 00082 s << std::endl; // start a new line 00083 } 00084 } 00085 00086 00088 00093 template<class V> 00094 void printvector (std::ostream& s, const V& v, std::string title, 00095 std::string rowtext, int columns=1, int width=10, 00096 int precision=2) 00097 { 00098 // count the numbers printed to make columns 00099 int counter=0; 00100 00101 // remember old flags 00102 std::ios_base::fmtflags oldflags = s.flags(); 00103 00104 // set the output format 00105 s.setf(std::ios_base::scientific, std::ios_base::floatfield); 00106 int oldprec = s.precision(); 00107 s.precision(precision); 00108 00109 // print title 00110 s << title << " [blocks=" << v.N() << ",dimension=" << v.dim() << "]" 00111 << std::endl; 00112 00113 // print data from all blocks 00114 recursive_printvector(s,v,rowtext,counter,columns,width,precision); 00115 00116 // check if new line is required 00117 if (counter%columns!=0) 00118 s << std::endl; 00119 00120 // reset the output format 00121 s.flags(oldflags); 00122 s.precision(oldprec); 00123 } 00124 00125 00127 // 00128 // pretty printing of matrices 00129 // 00130 00132 00137 inline void fill_row (std::ostream& s, int m, int width, int precision) 00138 { 00139 for (int j=0; j<m; j++) 00140 { 00141 s << " "; // space in front of each entry 00142 s.width(width); // set width for each entry anew 00143 s << "."; // yeah, the number ! 00144 } 00145 } 00146 00148 00153 template<class M> 00154 void print_row (std::ostream& s, const M& A, typename M::size_type I, 00155 typename M::size_type J, typename M::size_type therow, 00156 int width, int precision) 00157 { 00158 typename M::size_type i0=I; 00159 for (typename M::size_type i=0; i<A.N(); i++) 00160 { 00161 if (therow>=i0 && therow<i0+MatrixDimension<M>::rowdim(A,i)) 00162 { 00163 // the row is in this block row ! 00164 typename M::size_type j0=J; 00165 for (typename M::size_type j=0; j<A.M(); j++) 00166 { 00167 // find this block 00168 typename M::ConstColIterator it = A[i].find(j); 00169 00170 // print row or filler 00171 if (it!=A[i].end()) 00172 print_row(s,*it,i0,j0,therow,width,precision); 00173 else 00174 fill_row(s,MatrixDimension<M>::coldim(A,j),width,precision); 00175 00176 // advance columns 00177 j0 += MatrixDimension<M>::coldim(A,j); 00178 } 00179 } 00180 // advance rows 00181 i0 += MatrixDimension<M>::rowdim(A,i); 00182 } 00183 } 00184 00186 00191 template<class K, int n, int m> 00192 void print_row (std::ostream& s, const FieldMatrix<K,n,m>& A, 00193 typename FieldMatrix<K,n,m>::size_type I, 00194 typename FieldMatrix<K,n,m>::size_type J, 00195 typename FieldMatrix<K,n,m>::size_type therow, int width, 00196 int precision) 00197 { 00198 typedef typename FieldMatrix<K,n,m>::size_type size_type; 00199 00200 for (size_type i=0; i<n; i++) 00201 if (I+i==therow) 00202 for (int j=0; j<m; j++) 00203 { 00204 s << " "; // space in front of each entry 00205 s.width(width); // set width for each entry anew 00206 s << A[i][j]; // yeah, the number ! 00207 } 00208 } 00209 00211 00216 template<class K> 00217 void print_row (std::ostream& s, const FieldMatrix<K,1,1>& A, 00218 typename FieldMatrix<K,1,1>::size_type I, 00219 typename FieldMatrix<K,1,1>::size_type J, 00220 typename FieldMatrix<K,1,1>::size_type therow, 00221 int width, int precision) 00222 { 00223 if (I==therow) 00224 { 00225 s << " "; // space in front of each entry 00226 s.width(width); // set width for each entry anew 00227 s << static_cast<K>(A); // yeah, the number ! 00228 } 00229 } 00230 00232 00238 template<class M> 00239 void printmatrix (std::ostream& s, const M& A, std::string title, 00240 std::string rowtext, int width=10, int precision=2) 00241 { 00242 00243 // remember old flags 00244 std::ios_base::fmtflags oldflags = s.flags(); 00245 00246 // set the output format 00247 s.setf(std::ios_base::scientific, std::ios_base::floatfield); 00248 int oldprec = s.precision(); 00249 s.precision(precision); 00250 00251 // print title 00252 s << title 00253 << " [n=" << A.N() 00254 << ",m=" << A.M() 00255 << ",rowdim=" << MatrixDimension<M>::rowdim(A) 00256 << ",coldim=" << MatrixDimension<M>::coldim(A) 00257 << "]" << std::endl; 00258 00259 // print all rows 00260 for (typename M::size_type i=0; i<MatrixDimension<M>::rowdim(A); i++) 00261 { 00262 s << rowtext; // start a new row 00263 s << " "; // space in front of each entry 00264 s.width(4); // set width for counter 00265 s << i; // number of first entry in a line 00266 print_row(s,A,0,0,i,width,precision); // generic print 00267 s << std::endl;// start a new line 00268 } 00269 00270 // reset the output format 00271 s.flags(oldflags); 00272 s.precision(oldprec); 00273 } 00274 00276 00295 template<class B, int n, int m, class A> 00296 void printSparseMatrix(std::ostream& s, 00297 const BCRSMatrix<FieldMatrix<B,n,m>,A>& mat, 00298 std::string title, std::string rowtext, 00299 int width=3, int precision=2) 00300 { 00301 typedef BCRSMatrix<FieldMatrix<B,n,m>,A> Matrix; 00302 // remember old flags 00303 std::ios_base::fmtflags oldflags = s.flags(); 00304 // set the output format 00305 s.setf(std::ios_base::scientific, std::ios_base::floatfield); 00306 int oldprec = s.precision(); 00307 s.precision(precision); 00308 // print title 00309 s << title 00310 << " [n=" << mat.N() 00311 << ",m=" << mat.M() 00312 << ",rowdim=" << MatrixDimension<Matrix>::rowdim(mat) 00313 << ",coldim=" << MatrixDimension<Matrix>::coldim(mat) 00314 << "]" << std::endl; 00315 00316 typedef typename Matrix::ConstRowIterator Row; 00317 00318 for(Row row=mat.begin(); row != mat.end();++row) { 00319 int skipcols=0; 00320 bool reachedEnd=false; 00321 00322 while(!reachedEnd) { 00323 for(int innerrow=0; innerrow<n; ++innerrow) { 00324 int count=0; 00325 typedef typename Matrix::ConstColIterator Col; 00326 Col col=row->begin(); 00327 for(; col != row->end(); ++col,++count) { 00328 if(count<skipcols) 00329 continue; 00330 if(count>=skipcols+width) 00331 break; 00332 if(innerrow==0) { 00333 if(count==skipcols) { 00334 s << rowtext; // start a new row 00335 s << " "; // space in front of each entry 00336 s.width(4); // set width for counter 00337 s << row.index()<<": "; // number of first entry in a line 00338 } 00339 s.width(4); 00340 s<<col.index()<<": |"; 00341 } else { 00342 if(count==skipcols){ 00343 for(typename std::string::size_type i=0; i < rowtext.length(); i++) 00344 s<<" "; 00345 s<<" "; 00346 } 00347 s<<" |"; 00348 } 00349 for(int innercol=0; innercol < m; ++innercol) { 00350 s.width(9); 00351 s<<(*col)[innerrow][innercol]<<" "; 00352 } 00353 00354 s<<"|"; 00355 } 00356 if(innerrow==n-1 && col==row->end()) 00357 reachedEnd = true; 00358 else 00359 s << std::endl; 00360 } 00361 skipcols += width; 00362 s << std::endl; 00363 } 00364 s << std::endl; 00365 } 00366 00367 // reset the output format 00368 s.flags(oldflags); 00369 s.precision(oldprec); 00370 } 00371 00372 namespace 00373 { 00374 template<typename T> 00375 struct MatlabPODWriter 00376 { 00377 static std::ostream& write(const T& t, std::ostream& s) 00378 { 00379 s << t; 00380 return s; 00381 } 00382 }; 00383 template<typename T> 00384 struct MatlabPODWriter<std::complex<T> > 00385 { 00386 static std::ostream& write(const std::complex<T>& t, std::ostream& s) 00387 { 00388 s << t.real() << " " << t.imag(); 00389 return s; 00390 } 00391 }; 00392 } // anonymous namespace 00393 00395 00402 template <class FieldType, int rows, int cols> 00403 void writeMatrixToMatlabHelper 00404 ( const FieldMatrix<FieldType,rows,cols>& matrix, int rowOffset, 00405 int colOffset, std::ostream& s) 00406 { 00407 for (int i=0; i<rows; i++) 00408 for (int j=0; j<cols; j++){ 00409 //+1 for Matlab numbering 00410 s << rowOffset + i + 1 << " " << colOffset + j + 1 << " "; 00411 MatlabPODWriter<FieldType>::write(matrix[i][j], s)<< std::endl; 00412 } 00413 } 00414 00416 00421 template <class MatrixType> 00422 void writeMatrixToMatlabHelper(const MatrixType& matrix, 00423 int externalRowOffset, int externalColOffset, 00424 std::ostream& s) 00425 { 00426 // Precompute the accumulated sizes of the columns 00427 std::vector<typename MatrixType::size_type> colOffset(matrix.M()); 00428 if (colOffset.size() > 0) 00429 colOffset[0] = 0; 00430 00431 for (typename MatrixType::size_type i=0; i<matrix.M()-1; i++) 00432 colOffset[i+1] = colOffset[i] + 00433 MatrixDimension<MatrixType>::coldim(matrix,i); 00434 00435 typename MatrixType::size_type rowOffset = 0; 00436 00437 // Loop over all matrix rows 00438 for (typename MatrixType::size_type rowIdx=0; rowIdx<matrix.N(); rowIdx++) 00439 { 00440 00441 const typename MatrixType::row_type& row = matrix[rowIdx]; 00442 00443 typename MatrixType::row_type::ConstIterator cIt = row.begin(); 00444 typename MatrixType::row_type::ConstIterator cEndIt = row.end(); 00445 00446 // Loop over all columns in this row 00447 for (; cIt!=cEndIt; ++cIt) 00448 writeMatrixToMatlabHelper(*cIt, 00449 externalRowOffset+rowOffset, 00450 externalColOffset + colOffset[cIt.index()], 00451 s); 00452 00453 rowOffset += MatrixDimension<MatrixType>::rowdim(matrix, rowIdx); 00454 } 00455 00456 } 00457 00459 00476 template <class MatrixType> 00477 void writeMatrixToMatlab(const MatrixType& matrix, 00478 const std::string& filename, int outputPrecision = 18) 00479 { 00480 std::ofstream outStream(filename.c_str()); 00481 int oldPrecision = outStream.precision(); 00482 outStream.precision(outputPrecision); 00483 00484 writeMatrixToMatlabHelper(matrix, 0, 0, outStream); 00485 outStream.precision(oldPrecision); 00486 } 00487 00490 } // namespace Dune 00491 00492 #endif