dune-istl
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*- 00002 // vi: set et ts=8 sw=4 sts=4: 00003 #ifndef DUNE_MATRIX_HH 00004 #define DUNE_MATRIX_HH 00005 00010 #include <vector> 00011 #include <memory> 00012 00013 #include <dune/istl/vbvector.hh> 00014 00015 namespace Dune { 00016 00022 template<class T, class A=std::allocator<T> > 00023 class Matrix 00024 { 00025 public: 00026 00028 typedef typename T::field_type field_type; 00029 00031 typedef T block_type; 00032 00034 typedef A allocator_type; 00035 00037 typedef typename VariableBlockVector<T,A>::window_type row_type; 00038 00040 typedef typename A::size_type size_type; 00041 00043 typedef typename VariableBlockVector<T,A>::Iterator RowIterator; 00044 00046 typedef typename row_type::iterator ColIterator; 00047 00049 typedef typename VariableBlockVector<T,A>::ConstIterator ConstRowIterator; 00050 00052 typedef typename row_type::const_iterator ConstColIterator; 00053 00054 enum { 00056 blocklevel = T::blocklevel+1 00057 }; 00058 00060 Matrix() : data_(0), cols_(0) 00061 {} 00062 00065 Matrix(size_type rows, size_type cols) : data_(rows,cols), cols_(cols) 00066 {} 00067 00072 void setSize(size_type rows, size_type cols) { 00073 data_.resize(rows,cols); 00074 cols_ = cols; 00075 } 00076 00078 RowIterator begin() 00079 { 00080 return data_.begin(); 00081 } 00082 00084 RowIterator end() 00085 { 00086 return data_.end(); 00087 } 00088 00091 RowIterator beforeEnd () 00092 { 00093 return data_.beforeEnd(); 00094 } 00095 00098 RowIterator beforeBegin () 00099 { 00100 return data_.beforeBegin(); 00101 } 00102 00104 ConstRowIterator begin() const 00105 { 00106 return data_.begin(); 00107 } 00108 00110 ConstRowIterator end() const 00111 { 00112 return data_.end(); 00113 } 00114 00117 ConstRowIterator beforeEnd() const 00118 { 00119 return data_.beforeEnd(); 00120 } 00121 00124 ConstRowIterator beforeBegin () const 00125 { 00126 return data_.beforeBegin(); 00127 } 00128 00130 Matrix& operator= (const field_type& t) 00131 { 00132 data_ = t; 00133 return *this; 00134 } 00135 00137 row_type& operator[](size_type row) { 00138 #ifdef DUNE_ISTL_WITH_CHECKING 00139 if (row<0) 00140 DUNE_THROW(ISTLError, "Can't access negative rows!"); 00141 if (row>=N()) 00142 DUNE_THROW(ISTLError, "Row index out of range!"); 00143 #endif 00144 return data_[row]; 00145 } 00146 00148 const row_type& operator[](size_type row) const { 00149 #ifdef DUNE_ISTL_WITH_CHECKING 00150 if (row<0) 00151 DUNE_THROW(ISTLError, "Can't access negative rows!"); 00152 if (row>=N()) 00153 DUNE_THROW(ISTLError, "Row index out of range!"); 00154 #endif 00155 return data_[row]; 00156 } 00157 00159 size_type N() const { 00160 return data_.N(); 00161 } 00162 00164 size_type M() const { 00165 return cols_; 00166 } 00167 00169 size_type rowdim() const { 00170 #ifdef DUNE_ISTL_WITH_CHECKING 00171 if (M()==0) 00172 DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!"); 00173 #endif 00174 size_type dim = 0; 00175 for (size_type i=0; i<data_.N(); i++) 00176 dim += data_[i][0].rowdim(); 00177 00178 return dim; 00179 } 00180 00182 size_type coldim() const { 00183 #ifdef DUNE_ISTL_WITH_CHECKING 00184 if (N()==0) 00185 DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!"); 00186 #endif 00187 size_type dim = 0; 00188 for (size_type i=0; i<data_[0].size(); i++) 00189 dim += data_[0][i].coldim(); 00190 00191 return dim; 00192 } 00193 00195 size_type rowdim(size_type r) const { 00196 #ifdef DUNE_ISTL_WITH_CHECKING 00197 if (r<0 || r>=N()) 00198 DUNE_THROW(ISTLError, "Rowdim for nonexisting row " << r << " requested!"); 00199 if (M()==0) 00200 DUNE_THROW(ISTLError, "Can't compute rowdim() when there are no columns!"); 00201 #endif 00202 return data_[r][0].rowdim(); 00203 } 00204 00206 size_type coldim(size_type c) const { 00207 #ifdef DUNE_ISTL_WITH_CHECKING 00208 if (c<0 || c>=M()) 00209 DUNE_THROW(ISTLError, "Coldim for nonexisting column " << c << " requested!"); 00210 if (N()==0) 00211 DUNE_THROW(ISTLError, "Can't compute coldim() when there are no rows!"); 00212 #endif 00213 return data_[0][c].coldim(); 00214 } 00215 00217 Matrix<T>& operator*=(const field_type& scalar) { 00218 data_ *= scalar; 00219 return (*this); 00220 } 00221 00223 Matrix<T>& operator/=(const field_type& scalar) { 00224 data_ /= scalar; 00225 return (*this); 00226 } 00227 00233 Matrix& operator+= (const Matrix& b) { 00234 #ifdef DUNE_ISTL_WITH_CHECKING 00235 if(N()!=b.N() || M() != b.M()) 00236 DUNE_THROW(RangeError, "Matrix sizes do not match!"); 00237 #endif 00238 data_ += b.data_; 00239 return (*this); 00240 } 00241 00247 Matrix& operator-= (const Matrix& b) { 00248 #ifdef DUNE_ISTL_WITH_CHECKING 00249 if(N()!=b.N() || M() != b.M()) 00250 DUNE_THROW(RangeError, "Matrix sizes do not match!"); 00251 #endif 00252 data_ -= b.data_; 00253 return (*this); 00254 } 00255 00257 Matrix transpose() const { 00258 Matrix out(N(), M()); 00259 for (size_type i=0; i<M(); i++) 00260 for (size_type j=0; j<N(); j++) 00261 out[j][i] = (*this)[i][j]; 00262 00263 return out; 00264 } 00265 00269 template <class X, class Y> 00270 Y transposedMult(const X& vec) DUNE_DEPRECATED; 00271 00273 friend Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) { 00274 Matrix<T> out(m1.N(), m2.M()); 00275 out = 0; 00276 00277 for (size_type i=0; i<out.N(); i++ ) { 00278 for ( size_type j=0; j<out.M(); j++ ) 00279 for (size_type k=0; k<m1.M(); k++) 00280 out[i][j] += m1[i][k]*m2[k][j]; 00281 } 00282 00283 return out; 00284 } 00285 00287 template <class X, class Y> 00288 friend Y operator*(const Matrix<T>& m, const X& vec) { 00289 #ifdef DUNE_ISTL_WITH_CHECKING 00290 if (m.M()!=vec.size()) 00291 DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix columns!"); 00292 #endif 00293 Y out(m.N()); 00294 out = 0; 00295 00296 for (size_type i=0; i<out.size(); i++ ) { 00297 for ( size_type j=0; j<vec.size(); j++ ) 00298 out[i] += m[i][j]*vec[j]; 00299 } 00300 00301 return out; 00302 } 00303 00305 template <class X, class Y> 00306 void mv(const X& x, Y& y) const 00307 { 00308 #ifdef DUNE_ISTL_WITH_CHECKING 00309 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00310 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00311 #endif 00312 00313 for (size_type i=0; i<data_.N(); i++) { 00314 y[i]=0; 00315 for (size_type j=0; j<cols_; j++) 00316 (*this)[i][j].umv(x[j], y[i]); 00317 00318 } 00319 00320 } 00321 00323 template<class X, class Y> 00324 void mtv (const X& x, Y& y) const 00325 { 00326 #ifdef DUNE_ISTL_WITH_CHECKING 00327 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 00328 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 00329 #endif 00330 00331 for(size_type i=0; i<y.N(); ++i) 00332 y[i]=0; 00333 umtv(x,y); 00334 } 00335 00337 template <class X, class Y> 00338 void umv(const X& x, Y& y) const 00339 { 00340 #ifdef DUNE_ISTL_WITH_CHECKING 00341 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00342 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00343 #endif 00344 00345 for (size_type i=0; i<data_.N(); i++) { 00346 00347 for (size_type j=0; j<cols_; j++) 00348 (*this)[i][j].umv(x[j], y[i]); 00349 00350 } 00351 00352 } 00353 00355 template<class X, class Y> 00356 void mmv (const X& x, Y& y) const 00357 { 00358 #ifdef DUNE_ISTL_WITH_CHECKING 00359 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00360 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00361 #endif 00362 ConstRowIterator endi=end(); 00363 for (ConstRowIterator i=begin(); i!=endi; ++i) 00364 { 00365 ConstColIterator endj = (*i).end(); 00366 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00367 (*j).mmv(x[j.index()],y[i.index()]); 00368 } 00369 } 00370 00372 template <class X, class Y> 00373 void usmv(const field_type& alpha, const X& x, Y& y) const 00374 { 00375 #ifdef DUNE_ISTL_WITH_CHECKING 00376 if (x.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00377 if (y.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00378 #endif 00379 00380 for (size_type i=0; i<data_.N(); i++) { 00381 00382 for (size_type j=0; j<cols_; j++) 00383 (*this)[i][j].usmv(alpha, x[j], y[i]); 00384 00385 } 00386 00387 } 00388 00390 template<class X, class Y> 00391 void umtv (const X& x, Y& y) const 00392 { 00393 #ifdef DUNE_ISTL_WITH_CHECKING 00394 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00395 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00396 #endif 00397 ConstRowIterator endi=end(); 00398 for (ConstRowIterator i=begin(); i!=endi; ++i) 00399 { 00400 ConstColIterator endj = (*i).end(); 00401 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00402 (*j).umtv(x[i.index()],y[j.index()]); 00403 } 00404 } 00405 00407 template<class X, class Y> 00408 void mmtv (const X& x, Y& y) const 00409 { 00410 #ifdef DUNE_ISTL_WITH_CHECKING 00411 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00412 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00413 #endif 00414 ConstRowIterator endi=end(); 00415 for (ConstRowIterator i=begin(); i!=endi; ++i) 00416 { 00417 ConstColIterator endj = (*i).end(); 00418 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00419 (*j).mmtv(x[i.index()],y[j.index()]); 00420 } 00421 } 00422 00424 template<class X, class Y> 00425 void usmtv (const field_type& alpha, const X& x, Y& y) const 00426 { 00427 #ifdef DUNE_ISTL_WITH_CHECKING 00428 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00429 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00430 #endif 00431 ConstRowIterator endi=end(); 00432 for (ConstRowIterator i=begin(); i!=endi; ++i) 00433 { 00434 ConstColIterator endj = (*i).end(); 00435 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00436 (*j).usmtv(alpha,x[i.index()],y[j.index()]); 00437 } 00438 } 00439 00441 template<class X, class Y> 00442 void umhv (const X& x, Y& y) const 00443 { 00444 #ifdef DUNE_ISTL_WITH_CHECKING 00445 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00446 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00447 #endif 00448 ConstRowIterator endi=end(); 00449 for (ConstRowIterator i=begin(); i!=endi; ++i) 00450 { 00451 ConstColIterator endj = (*i).end(); 00452 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00453 (*j).umhv(x[i.index()],y[j.index()]); 00454 } 00455 } 00456 00458 template<class X, class Y> 00459 void mmhv (const X& x, Y& y) const 00460 { 00461 #ifdef DUNE_ISTL_WITH_CHECKING 00462 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00463 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00464 #endif 00465 ConstRowIterator endi=end(); 00466 for (ConstRowIterator i=begin(); i!=endi; ++i) 00467 { 00468 ConstColIterator endj = (*i).end(); 00469 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00470 (*j).mmhv(x[i.index()],y[j.index()]); 00471 } 00472 } 00473 00475 template<class X, class Y> 00476 void usmhv (const field_type& alpha, const X& x, Y& y) const 00477 { 00478 #ifdef DUNE_ISTL_WITH_CHECKING 00479 if (x.N()!=N()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00480 if (y.N()!=M()) DUNE_THROW(ISTLError,"vector/matrix size mismatch!"); 00481 #endif 00482 ConstRowIterator endi=end(); 00483 for (ConstRowIterator i=begin(); i!=endi; ++i) 00484 { 00485 ConstColIterator endj = (*i).end(); 00486 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 00487 (*j).usmhv(alpha,x[i.index()],y[j.index()]); 00488 } 00489 } 00490 00491 //===== norms 00492 00494 double frobenius_norm () const 00495 { 00496 return std::sqrt(frobenius_norm2()); 00497 } 00498 00500 double frobenius_norm2 () const 00501 { 00502 double sum=0; 00503 for (size_type i=0; i<N(); ++i) 00504 for (size_type j=0; j<M(); ++j) 00505 sum += data_[i][j].frobenius_norm2(); 00506 return sum; 00507 } 00508 00510 double infinity_norm () const 00511 { 00512 double max=0; 00513 for (size_type i=0; i<N(); ++i) { 00514 double sum=0; 00515 for (size_type j=0; j<M(); j++) 00516 sum += data_[i][j].infinity_norm(); 00517 max = std::max(max,sum); 00518 } 00519 return max; 00520 } 00521 00523 double infinity_norm_real () const 00524 { 00525 double max=0; 00526 for (size_type i=0; i<N(); ++i) { 00527 double sum=0; 00528 for (size_type j=0; j<M(); j++) 00529 sum += data_[i][j].infinity_norm_real(); 00530 max = std::max(max,sum); 00531 } 00532 return max; 00533 } 00534 00535 //===== query 00536 00538 bool exists (size_type i, size_type j) const 00539 { 00540 #ifdef DUNE_ISTL_WITH_CHECKING 00541 if (i<0 || i>=N()) DUNE_THROW(ISTLError,"row index out of range"); 00542 if (j<0 || i>=M()) DUNE_THROW(ISTLError,"column index out of range"); 00543 #endif 00544 return true; 00545 } 00546 00547 protected: 00548 00554 VariableBlockVector<T,A> data_; 00555 00561 size_type cols_; 00562 }; 00563 00564 template<class T, class A> 00565 template<class X, class Y> 00566 inline Y Matrix<T,A>::transposedMult(const X& vec) 00567 { 00568 #ifdef DUNE_ISTL_WITH_CHECKING 00569 if (N()!=vec.size()) 00570 DUNE_THROW(ISTLError, "Vector size doesn't match the number of matrix rows!"); 00571 #endif 00572 Y out(M()); 00573 out = 0; 00574 00575 for (size_type i=0; i<out.size(); i++ ) { 00576 for ( size_type j=0; j<vec.size(); j++ ) 00577 out[i] += (*this)[j][i]*vec[j]; 00578 } 00579 00580 return out; 00581 } 00582 00584 } // end namespace Dune 00585 00586 #endif