00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__matrix_lib_h
00013 #define __deal2__matrix_lib_h
00014
00015 #include <deal.II/base/subscriptor.h>
00016 #include <deal.II/lac/vector_memory.h>
00017 #include <deal.II/lac/pointer_matrix.h>
00018 #include <deal.II/lac/solver_richardson.h>
00019
00020 DEAL_II_NAMESPACE_OPEN
00021
00022 template<typename number> class Vector;
00023 template<typename number> class BlockVector;
00024 template<typename number> class SparseMatrix;
00025
00043 template<class VECTOR>
00044 class ProductMatrix : public PointerMatrixBase<VECTOR>
00045 {
00046 public:
00053 ProductMatrix();
00054
00060 ProductMatrix(VectorMemory<VECTOR>& mem);
00061
00068 template <class MATRIX1, class MATRIX2>
00069 ProductMatrix(const MATRIX1& m1,
00070 const MATRIX2& m2,
00071 VectorMemory<VECTOR>& mem);
00072
00076 template <class MATRIX1, class MATRIX2>
00077 void reinit(const MATRIX1& m1, const MATRIX2& m2);
00078
00082 template <class MATRIX1, class MATRIX2>
00083 void initialize(const MATRIX1& m1, const MATRIX2& m2,
00084 VectorMemory<VECTOR>& mem);
00085
00089 ~ProductMatrix();
00090
00091
00092 void clear();
00093
00098 virtual void vmult (VECTOR& w,
00099 const VECTOR& v) const;
00100
00106 virtual void Tvmult (VECTOR& w,
00107 const VECTOR& v) const;
00108
00113 virtual void vmult_add (VECTOR& w,
00114 const VECTOR& v) const;
00115
00122 virtual void Tvmult_add (VECTOR& w,
00123 const VECTOR& v) const;
00124
00125 private:
00129 PointerMatrixBase<VECTOR>* m1;
00130
00134 PointerMatrixBase<VECTOR>* m2;
00135
00139 SmartPointer<VectorMemory<VECTOR>,ProductMatrix<VECTOR> > mem;
00144 virtual const void* get() const;
00145 };
00146
00147
00156 template<class VECTOR>
00157 class ScaledMatrix : public Subscriptor
00158 {
00159 public:
00164 ScaledMatrix ();
00168 template <class MATRIX>
00169 ScaledMatrix (const MATRIX& M, const double factor);
00170
00174 ~ScaledMatrix ();
00179 template <class MATRIX>
00180 void initialize (const MATRIX& M, const double factor);
00181
00185 void clear ();
00186
00190 void vmult (VECTOR& w, const VECTOR& v) const;
00191
00196 void Tvmult (VECTOR& w, const VECTOR& v) const;
00197
00198 private:
00202 PointerMatrixBase<VECTOR>* m;
00206 double factor;
00207 };
00208
00209
00210
00223 template<typename number, typename vector_number>
00224 class ProductSparseMatrix : public PointerMatrixBase<Vector<vector_number> >
00225 {
00226 public:
00230 typedef SparseMatrix<number> MatrixType;
00231
00236 typedef Vector<vector_number> VectorType;
00237
00244 ProductSparseMatrix(const MatrixType& m1,
00245 const MatrixType& m2,
00246 VectorMemory<VectorType>& mem);
00247
00255 ProductSparseMatrix();
00256
00257 void initialize(const MatrixType& m1,
00258 const MatrixType& m2,
00259 VectorMemory<VectorType>& mem);
00260
00261
00262 void clear();
00263
00268 virtual void vmult (VectorType& w,
00269 const VectorType& v) const;
00270
00276 virtual void Tvmult (VectorType& w,
00277 const VectorType& v) const;
00278
00283 virtual void vmult_add (VectorType& w,
00284 const VectorType& v) const;
00285
00292 virtual void Tvmult_add (VectorType& w,
00293 const VectorType& v) const;
00294
00295 private:
00299 SmartPointer<const MatrixType,ProductSparseMatrix<number,vector_number> > m1;
00300
00304 SmartPointer<const MatrixType,ProductSparseMatrix<number,vector_number> > m2;
00305
00309 SmartPointer<VectorMemory<VectorType>,ProductSparseMatrix<number,vector_number> > mem;
00314 virtual const void* get() const;
00315 };
00316
00317
00326 class MeanValueFilter : public Subscriptor
00327 {
00328 public:
00333 MeanValueFilter(unsigned int component = numbers::invalid_unsigned_int);
00334
00338 template <typename number>
00339 void filter (Vector<number>& v) const;
00340
00344 template <typename number>
00345 void filter (BlockVector<number>& v) const;
00346
00351 template <typename number>
00352 void vmult (Vector<number>& dst,
00353 const Vector<number>& src) const;
00354
00359 template <typename number>
00360 void vmult_add (Vector<number>& dst,
00361 const Vector<number>& src) const;
00362
00368 template <typename number>
00369 void vmult (BlockVector<number>& dst,
00370 const BlockVector<number>& src) const;
00371
00377 template <typename number>
00378 void vmult_add (BlockVector<number>& dst,
00379 const BlockVector<number>& src) const;
00380
00381
00385 template <typename VECTOR>
00386 void Tvmult(VECTOR&, const VECTOR&) const;
00387
00391 template <typename VECTOR>
00392 void Tvmult_add(VECTOR&, const VECTOR&) const;
00393
00394 private:
00398 unsigned int component;
00399 };
00400
00401
00402
00420 template<class VECTOR>
00421 class InverseMatrixRichardson : public Subscriptor
00422 {
00423 public:
00431 InverseMatrixRichardson (SolverControl& control,
00432 VectorMemory<VECTOR>& mem);
00437 ~InverseMatrixRichardson();
00438
00445 template <class MATRIX, class PRECONDITION>
00446 void initialize (const MATRIX&,
00447 const PRECONDITION&);
00448
00453 SolverControl& control() const;
00457 void vmult (VECTOR&, const VECTOR&) const;
00458
00462 void vmult_add (VECTOR&, const VECTOR&) const;
00463
00467 void Tvmult (VECTOR&, const VECTOR&) const;
00468
00472 void Tvmult_add (VECTOR&, const VECTOR&) const;
00473
00474 private:
00479 VectorMemory<VECTOR> &mem;
00480
00484 mutable SolverRichardson<VECTOR> solver;
00485
00489 PointerMatrixBase<VECTOR>* matrix;
00490
00494 PointerMatrixBase<VECTOR>* precondition;
00495 };
00496
00497
00498
00499
00501
00502
00503
00504 template<class VECTOR>
00505 inline
00506 ScaledMatrix<VECTOR>::ScaledMatrix()
00507 :
00508 m(0)
00509 {}
00510
00511
00512
00513 template<class VECTOR>
00514 template<class MATRIX>
00515 inline
00516 ScaledMatrix<VECTOR>::ScaledMatrix(const MATRIX& mat, const double factor)
00517 :
00518 m(new_pointer_matrix_base(mat, VECTOR())),
00519 factor(factor)
00520 {}
00521
00522
00523
00524 template<class VECTOR>
00525 template<class MATRIX>
00526 inline
00527 void
00528 ScaledMatrix<VECTOR>::initialize(const MATRIX& mat, const double f)
00529 {
00530 if (m) delete m;
00531 m = new_pointer_matrix_base(mat, VECTOR());
00532 factor = f;
00533 }
00534
00535
00536
00537 template<class VECTOR>
00538 inline
00539 void
00540 ScaledMatrix<VECTOR>::clear()
00541 {
00542 if (m) delete m;
00543 m = 0;
00544 }
00545
00546
00547
00548 template<class VECTOR>
00549 inline
00550 ScaledMatrix<VECTOR>::~ScaledMatrix()
00551 {
00552 clear ();
00553 }
00554
00555
00556 template<class VECTOR>
00557 inline
00558 void
00559 ScaledMatrix<VECTOR>::vmult (VECTOR& w, const VECTOR& v) const
00560 {
00561 m->vmult(w, v);
00562 w *= factor;
00563 }
00564
00565
00566 template<class VECTOR>
00567 inline
00568 void
00569 ScaledMatrix<VECTOR>::Tvmult (VECTOR& w, const VECTOR& v) const
00570 {
00571 m->Tvmult(w, v);
00572 w *= factor;
00573 }
00574
00575
00576
00577
00578 template<class VECTOR>
00579 ProductMatrix<VECTOR>::ProductMatrix ()
00580 : m1(0), m2(0), mem(0)
00581 {}
00582
00583
00584 template<class VECTOR>
00585 ProductMatrix<VECTOR>::ProductMatrix (VectorMemory<VECTOR>& m)
00586 : m1(0), m2(0), mem(&m)
00587 {}
00588
00589
00590 template<class VECTOR>
00591 template<class MATRIX1, class MATRIX2>
00592 ProductMatrix<VECTOR>::ProductMatrix (
00593 const MATRIX1& mat1,
00594 const MATRIX2& mat2,
00595 VectorMemory<VECTOR>& m)
00596 : mem(&m)
00597 {
00598 m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1, typeid(*this).name());
00599 m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2, typeid(*this).name());
00600 }
00601
00602
00603 template<class VECTOR>
00604 template<class MATRIX1, class MATRIX2>
00605 void
00606 ProductMatrix<VECTOR>::reinit (
00607 const MATRIX1& mat1,
00608 const MATRIX2& mat2)
00609 {
00610 if (m1) delete m1;
00611 if (m2) delete m2;
00612 m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1, typeid(*this).name());
00613 m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2, typeid(*this).name());
00614 }
00615
00616
00617 template<class VECTOR>
00618 template<class MATRIX1, class MATRIX2>
00619 void
00620 ProductMatrix<VECTOR>::initialize (
00621 const MATRIX1& mat1,
00622 const MATRIX2& mat2,
00623 VectorMemory<VECTOR>& memory)
00624 {
00625 mem = &memory;
00626 if (m1) delete m1;
00627 if (m2) delete m2;
00628 m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1, typeid(*this).name());
00629 m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2, typeid(*this).name());
00630 }
00631
00632
00633 template<class VECTOR>
00634 ProductMatrix<VECTOR>::~ProductMatrix ()
00635 {
00636 if (m1) delete m1;
00637 if (m2) delete m2;
00638 }
00639
00640
00641 template<class VECTOR>
00642 void
00643 ProductMatrix<VECTOR>::clear ()
00644 {
00645 if (m1) delete m1;
00646 m1 = 0;
00647 if (m2) delete m2;
00648 m2 = 0;
00649 }
00650
00651
00652 template<class VECTOR>
00653 void
00654 ProductMatrix<VECTOR>::vmult (VECTOR& dst, const VECTOR& src) const
00655 {
00656 Assert (mem != 0, ExcNotInitialized());
00657 Assert (m1 != 0, ExcNotInitialized());
00658 Assert (m2 != 0, ExcNotInitialized());
00659
00660 VECTOR* v = mem->alloc();
00661 v->reinit(dst);
00662 m2->vmult (*v, src);
00663 m1->vmult (dst, *v);
00664 mem->free(v);
00665 }
00666
00667
00668 template<class VECTOR>
00669 void
00670 ProductMatrix<VECTOR>::vmult_add (VECTOR& dst, const VECTOR& src) const
00671 {
00672 Assert (mem != 0, ExcNotInitialized());
00673 Assert (m1 != 0, ExcNotInitialized());
00674 Assert (m2 != 0, ExcNotInitialized());
00675
00676 VECTOR* v = mem->alloc();
00677 v->reinit(dst);
00678 m2->vmult (*v, src);
00679 m1->vmult_add (dst, *v);
00680 mem->free(v);
00681 }
00682
00683
00684 template<class VECTOR>
00685 void
00686 ProductMatrix<VECTOR>::Tvmult (VECTOR& dst, const VECTOR& src) const
00687 {
00688 Assert (mem != 0, ExcNotInitialized());
00689 Assert (m1 != 0, ExcNotInitialized());
00690 Assert (m2 != 0, ExcNotInitialized());
00691
00692 VECTOR* v = mem->alloc();
00693 v->reinit(dst);
00694 m1->Tvmult (*v, src);
00695 m2->Tvmult (dst, *v);
00696 mem->free(v);
00697 }
00698
00699
00700 template<class VECTOR>
00701 void
00702 ProductMatrix<VECTOR>::Tvmult_add (VECTOR& dst, const VECTOR& src) const
00703 {
00704 Assert (mem != 0, ExcNotInitialized());
00705 Assert (m1 != 0, ExcNotInitialized());
00706 Assert (m2 != 0, ExcNotInitialized());
00707
00708 VECTOR* v = mem->alloc();
00709 v->reinit(dst);
00710 m1->Tvmult (*v, src);
00711 m2->Tvmult_add (dst, *v);
00712 mem->free(v);
00713 }
00714
00715
00716 template<class VECTOR>
00717 const void*
00718 ProductMatrix<VECTOR>::get () const
00719 {
00720 return (void*) m1;
00721 }
00722
00723
00724
00725
00726 template <class VECTOR>
00727 inline void
00728 MeanValueFilter::Tvmult(VECTOR&, const VECTOR&) const
00729 {
00730 Assert(false, ExcNotImplemented());
00731 }
00732
00733
00734 template <class VECTOR>
00735 inline void
00736 MeanValueFilter::Tvmult_add(VECTOR&, const VECTOR&) const
00737 {
00738 Assert(false, ExcNotImplemented());
00739 }
00740
00741
00742
00743 template <class VECTOR>
00744 template <class MATRIX, class PRECONDITION>
00745 inline void
00746 InverseMatrixRichardson<VECTOR>::initialize (const MATRIX& m, const PRECONDITION& p)
00747 {
00748 if (matrix != 0)
00749 delete matrix;
00750 matrix = new PointerMatrix<MATRIX, VECTOR>(&m);
00751 if (precondition != 0)
00752 delete precondition;
00753 precondition = new PointerMatrix<PRECONDITION, VECTOR>(&p);;
00754 }
00755
00756
00757 DEAL_II_NAMESPACE_CLOSE
00758
00759 #endif
00760