00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__mg_matrix_h
00013 #define __deal2__mg_matrix_h
00014
00015 #include <deal.II/lac/vector.h>
00016 #include <deal.II/lac/pointer_matrix.h>
00017 #include <deal.II/lac/sparse_matrix.h>
00018 #include <deal.II/multigrid/mg_base.h>
00019 #include <deal.II/base/mg_level_object.h>
00020 #include <deal.II/base/std_cxx1x/shared_ptr.h>
00021
00022 DEAL_II_NAMESPACE_OPEN
00023
00026
00027 namespace mg
00028 {
00037 template <class VECTOR = Vector<double> >
00038 class Matrix
00039 : public MGMatrixBase<VECTOR>
00040 {
00041 public:
00046 Matrix();
00047
00053 template <class MATRIX>
00054 Matrix(const MGLevelObject<MATRIX>& M);
00055
00062 template <class MATRIX>
00063 void
00064 initialize(const MGLevelObject<MATRIX>& M);
00065
00069 const PointerMatrixBase<VECTOR>& operator[] (unsigned int level) const;
00070
00071 virtual void vmult (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
00072 virtual void vmult_add (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
00073 virtual void Tvmult (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
00074 virtual void Tvmult_add (const unsigned int level, VECTOR& dst, const VECTOR& src) const;
00075
00079 std::size_t memory_consumption () const;
00080 private:
00081 MGLevelObject<std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > > matrices;
00082 };
00083
00084 }
00085
00095 template <class MATRIX = SparseMatrix<double>, class VECTOR = Vector<double> >
00096 class MGMatrix : public MGMatrixBase<VECTOR>
00097 {
00098 public:
00108 MGMatrix (MGLevelObject<MATRIX>* = 0);
00109
00117 void set_matrix (MGLevelObject<MATRIX>* M);
00118
00123 virtual void vmult (const unsigned int level,
00124 VECTOR& dst,
00125 const VECTOR& src) const;
00126
00131 virtual void vmult_add (const unsigned int level,
00132 VECTOR& dst,
00133 const VECTOR& src) const;
00134
00140 virtual void Tvmult (const unsigned int level,
00141 VECTOR& dst,
00142 const VECTOR& src) const;
00143
00149 virtual void Tvmult_add (const unsigned int level,
00150 VECTOR& dst,
00151 const VECTOR& src) const;
00152
00156 std::size_t memory_consumption () const;
00157
00158
00159 private:
00163 SmartPointer<MGLevelObject<MATRIX>,MGMatrix<MATRIX,VECTOR> > matrix;
00164 };
00165
00166
00178 template <class MATRIX, typename number>
00179 class MGMatrixSelect : public MGMatrixBase<Vector<number> >
00180 {
00181 public:
00189 MGMatrixSelect (const unsigned int row = 0,
00190 const unsigned int col = 0,
00191 MGLevelObject<MATRIX>* = 0);
00192
00200 void set_matrix (MGLevelObject<MATRIX>* M);
00201
00206 void select_block (const unsigned int row,
00207 const unsigned int col);
00208
00213 virtual void vmult (const unsigned int level,
00214 Vector<number>& dst,
00215 const Vector<number>& src) const;
00216
00221 virtual void vmult_add (const unsigned int level,
00222 Vector<number>& dst,
00223 const Vector<number>& src) const;
00224
00230 virtual void Tvmult (const unsigned int level,
00231 Vector<number>& dst,
00232 const Vector<number>& src) const;
00233
00239 virtual void Tvmult_add (const unsigned int level,
00240 Vector<number>& dst,
00241 const Vector<number>& src) const;
00242
00243 private:
00247 SmartPointer<MGLevelObject<MATRIX>,MGMatrixSelect<MATRIX,number> > matrix;
00251 unsigned int row;
00255 unsigned int col;
00256
00257 };
00258
00261
00262
00263 namespace mg
00264 {
00265 template <class VECTOR>
00266 template <class MATRIX>
00267 inline
00268 void
00269 Matrix<VECTOR>::initialize (const MGLevelObject<MATRIX>& p)
00270 {
00271 matrices.resize(p.min_level(), p.max_level());
00272 for (unsigned int level=p.min_level();level <= p.max_level();++level)
00273 matrices[level] = std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > (new_pointer_matrix_base(p[level], VECTOR()));
00274 }
00275
00276
00277 template <class VECTOR>
00278 template <class MATRIX>
00279 inline
00280 Matrix<VECTOR>::Matrix (const MGLevelObject<MATRIX>& p)
00281 {
00282 initialize(p);
00283 }
00284
00285
00286 template <class VECTOR>
00287 inline
00288 Matrix<VECTOR>::Matrix ()
00289 {}
00290
00291
00292 template <class VECTOR>
00293 inline
00294 const PointerMatrixBase<VECTOR>&
00295 Matrix<VECTOR>::operator[] (unsigned int level) const
00296 {
00297 return *matrices[level];
00298 }
00299
00300
00301 template <class VECTOR>
00302 void
00303 Matrix<VECTOR>::vmult (
00304 const unsigned int level,
00305 VECTOR& dst,
00306 const VECTOR& src) const
00307 {
00308 matrices[level]->vmult(dst, src);
00309 }
00310
00311
00312 template <class VECTOR>
00313 void
00314 Matrix<VECTOR>::vmult_add (
00315 const unsigned int level,
00316 VECTOR& dst,
00317 const VECTOR& src) const
00318 {
00319 matrices[level]->vmult_add(dst, src);
00320 }
00321
00322
00323 template <class VECTOR>
00324 void
00325 Matrix<VECTOR>::Tvmult (const unsigned int level,
00326 VECTOR& dst,
00327 const VECTOR& src) const
00328 {
00329 matrices[level]->Tvmult(dst, src);
00330 }
00331
00332
00333 template <class VECTOR>
00334 void
00335 Matrix<VECTOR>::Tvmult_add (const unsigned int level,
00336 VECTOR& dst,
00337 const VECTOR& src) const
00338 {
00339 matrices[level]->Tvmult_add(dst, src);
00340 }
00341
00342
00343 template <class VECTOR>
00344 inline
00345 std::size_t
00346 Matrix<VECTOR>::memory_consumption () const
00347 {
00348 return sizeof(*this) + matrices->memory_consumption();
00349 }
00350 }
00351
00352
00353
00354 template <class MATRIX, class VECTOR>
00355 MGMatrix<MATRIX, VECTOR>::MGMatrix (MGLevelObject<MATRIX>* p)
00356 :
00357 matrix (p, typeid(*this).name())
00358 {}
00359
00360
00361 template <class MATRIX, class VECTOR>
00362 void
00363 MGMatrix<MATRIX, VECTOR>::set_matrix (MGLevelObject<MATRIX>* p)
00364 {
00365 matrix = p;
00366 }
00367
00368
00369 template <class MATRIX, class VECTOR>
00370 void
00371 MGMatrix<MATRIX, VECTOR>::vmult (const unsigned int level,
00372 VECTOR& dst,
00373 const VECTOR& src) const
00374 {
00375 Assert(matrix != 0, ExcNotInitialized());
00376
00377 const MGLevelObject<MATRIX>& m = *matrix;
00378 m[level].vmult(dst, src);
00379 }
00380
00381
00382 template <class MATRIX, class VECTOR>
00383 void
00384 MGMatrix<MATRIX, VECTOR>::vmult_add (const unsigned int level,
00385 VECTOR& dst,
00386 const VECTOR& src) const
00387 {
00388 Assert(matrix != 0, ExcNotInitialized());
00389
00390 const MGLevelObject<MATRIX>& m = *matrix;
00391 m[level].vmult_add(dst, src);
00392 }
00393
00394
00395 template <class MATRIX, class VECTOR>
00396 void
00397 MGMatrix<MATRIX, VECTOR>::Tvmult (const unsigned int level,
00398 VECTOR& dst,
00399 const VECTOR& src) const
00400 {
00401 Assert(matrix != 0, ExcNotInitialized());
00402
00403 const MGLevelObject<MATRIX>& m = *matrix;
00404 m[level].Tvmult(dst, src);
00405 }
00406
00407
00408 template <class MATRIX, class VECTOR>
00409 void
00410 MGMatrix<MATRIX, VECTOR>::Tvmult_add (const unsigned int level,
00411 VECTOR& dst,
00412 const VECTOR& src) const
00413 {
00414 Assert(matrix != 0, ExcNotInitialized());
00415
00416 const MGLevelObject<MATRIX>& m = *matrix;
00417 m[level].Tvmult_add(dst, src);
00418 }
00419
00420
00421 template <class MATRIX, class VECTOR>
00422 std::size_t
00423 MGMatrix<MATRIX, VECTOR>::memory_consumption () const
00424 {
00425 return sizeof(*this) + matrix->memory_consumption();
00426 }
00427
00428
00429
00430 template <class MATRIX, typename number>
00431 MGMatrixSelect<MATRIX, number>::
00432 MGMatrixSelect (const unsigned int row,
00433 const unsigned int col,
00434 MGLevelObject<MATRIX>* p)
00435 :
00436 matrix (p, typeid(*this).name()),
00437 row(row),
00438 col(col)
00439 {}
00440
00441
00442
00443 template <class MATRIX, typename number>
00444 void
00445 MGMatrixSelect<MATRIX, number>::set_matrix (MGLevelObject<MATRIX>* p)
00446 {
00447 matrix = p;
00448 }
00449
00450
00451 template <class MATRIX, typename number>
00452 void
00453 MGMatrixSelect<MATRIX, number>::
00454 select_block (const unsigned int brow,
00455 const unsigned int bcol)
00456 {
00457 row = brow;
00458 col = bcol;
00459 }
00460
00461
00462 template <class MATRIX, typename number>
00463 void
00464 MGMatrixSelect<MATRIX, number>::
00465 vmult (const unsigned int level,
00466 Vector<number>& dst,
00467 const Vector<number>& src) const
00468 {
00469 Assert(matrix != 0, ExcNotInitialized());
00470
00471 const MGLevelObject<MATRIX>& m = *matrix;
00472 m[level].block(row, col).vmult(dst, src);
00473 }
00474
00475
00476 template <class MATRIX, typename number>
00477 void
00478 MGMatrixSelect<MATRIX, number>::
00479 vmult_add (const unsigned int level,
00480 Vector<number>& dst,
00481 const Vector<number>& src) const
00482 {
00483 Assert(matrix != 0, ExcNotInitialized());
00484
00485 const MGLevelObject<MATRIX>& m = *matrix;
00486 m[level].block(row, col).vmult_add(dst, src);
00487 }
00488
00489
00490 template <class MATRIX, typename number>
00491 void
00492 MGMatrixSelect<MATRIX, number>::
00493 Tvmult (const unsigned int level,
00494 Vector<number>& dst,
00495 const Vector<number>& src) const
00496 {
00497 Assert(matrix != 0, ExcNotInitialized());
00498
00499 const MGLevelObject<MATRIX>& m = *matrix;
00500 m[level].block(row, col).Tvmult(dst, src);
00501 }
00502
00503
00504 template <class MATRIX, typename number>
00505 void
00506 MGMatrixSelect<MATRIX, number>::
00507 Tvmult_add (const unsigned int level,
00508 Vector<number>& dst,
00509 const Vector<number>& src) const
00510 {
00511 Assert(matrix != 0, ExcNotInitialized());
00512
00513 const MGLevelObject<MATRIX>& m = *matrix;
00514 m[level].block(row, col).Tvmult_add(dst, src);
00515 }
00516
00517 DEAL_II_NAMESPACE_CLOSE
00518
00519 #endif
00520