00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__pointer_matrix_h
00013 #define __deal2__pointer_matrix_h
00014
00015 #include <deal.II/base/subscriptor.h>
00016 #include <deal.II/base/smartpointer.h>
00017 #include <deal.II/lac/vector.h>
00018 #include <deal.II/lac/vector_memory.h>
00019
00020 DEAL_II_NAMESPACE_OPEN
00021
00022 template<class VECTOR> class VectorMemory;
00023
00024 class IdentityMatrix;
00025 template <typename number> class FullMatrix;
00026 template <typename number> class LAPACKFullMatrix;
00027 template <typename number> class SparseMatrix;
00028 template <typename number> class BlockSparseMatrix;
00029 template <typename number> class SparseMatrixEZ;
00030 template <typename number> class BlockSparseMatrixEZ;
00031 template <typename number> class BlockMatrixArray;
00032 template <typename number> class TridiagonalMatrix;
00033
00034
00047 template<class VECTOR>
00048 class PointerMatrixBase : public Subscriptor
00049 {
00050 public:
00065 typedef typename VECTOR::value_type value_type;
00066
00073 virtual ~PointerMatrixBase ();
00074
00079 virtual void clear () = 0;
00080
00085 bool operator == (const PointerMatrixBase<VECTOR>&) const;
00086
00091 bool operator != (const PointerMatrixBase<VECTOR>&) const;
00092
00097 bool operator < (const PointerMatrixBase<VECTOR>&) const;
00098
00103 bool operator <= (const PointerMatrixBase<VECTOR>&) const;
00104
00109 bool operator > (const PointerMatrixBase<VECTOR>&) const;
00110
00115 bool operator >= (const PointerMatrixBase<VECTOR>&) const;
00116
00117
00121 virtual void vmult (VECTOR& dst,
00122 const VECTOR& src) const = 0;
00123
00127 virtual void Tvmult (VECTOR& dst,
00128 const VECTOR& src) const = 0;
00129
00134 virtual void vmult_add (VECTOR& dst,
00135 const VECTOR& src) const = 0;
00136
00141 virtual void Tvmult_add (VECTOR& dst,
00142 const VECTOR& src) const = 0;
00143
00144 private:
00148 virtual const void* get() const = 0;
00149 };
00150
00161 template<class MATRIX, class VECTOR>
00162 class PointerMatrix : public PointerMatrixBase<VECTOR>
00163 {
00164 public:
00175 PointerMatrix (const MATRIX* M=0);
00176
00182 PointerMatrix(const char* name);
00183
00193 PointerMatrix(const MATRIX* M,
00194 const char* name);
00195
00196
00197 virtual void clear();
00198
00203 bool empty () const;
00204
00211 const PointerMatrix& operator= (const MATRIX* M);
00212
00216 virtual void vmult (VECTOR& dst,
00217 const VECTOR& src) const;
00218
00222 virtual void Tvmult (VECTOR& dst,
00223 const VECTOR& src) const;
00224
00229 virtual void vmult_add (VECTOR& dst,
00230 const VECTOR& src) const;
00231
00236 virtual void Tvmult_add (VECTOR& dst,
00237 const VECTOR& src) const;
00238
00239 private:
00244 virtual const void* get() const;
00245
00249 SmartPointer<const MATRIX,PointerMatrix<MATRIX,VECTOR> > m;
00250 };
00251
00252
00268 template<class MATRIX, class VECTOR>
00269 class PointerMatrixAux : public PointerMatrixBase<VECTOR>
00270 {
00271 public:
00286 PointerMatrixAux (VectorMemory<VECTOR>* mem = 0,
00287 const MATRIX* M=0);
00288
00295 PointerMatrixAux(VectorMemory<VECTOR>* mem,
00296 const char* name);
00297
00307 PointerMatrixAux(VectorMemory<VECTOR>* mem,
00308 const MATRIX* M,
00309 const char* name);
00310
00311
00312 virtual void clear();
00313
00318 bool empty () const;
00319
00325 void set_memory(VectorMemory<VECTOR>* mem);
00326
00333 const PointerMatrixAux& operator= (const MATRIX* M);
00334
00338 virtual void vmult (VECTOR& dst,
00339 const VECTOR& src) const;
00340
00344 virtual void Tvmult (VECTOR& dst,
00345 const VECTOR& src) const;
00346
00351 virtual void vmult_add (VECTOR& dst,
00352 const VECTOR& src) const;
00353
00358 virtual void Tvmult_add (VECTOR& dst,
00359 const VECTOR& src) const;
00360
00361 private:
00366 virtual const void* get() const;
00367
00371 mutable GrowingVectorMemory<VECTOR> my_memory;
00372
00377 mutable SmartPointer<VectorMemory<VECTOR>,PointerMatrixAux<MATRIX,VECTOR> > mem;
00378
00382 SmartPointer<const MATRIX,PointerMatrixAux<MATRIX,VECTOR> > m;
00383 };
00384
00385
00386
00398 template <typename number>
00399 class PointerMatrixVector : public PointerMatrixBase<Vector<number> >
00400 {
00401 public:
00412 PointerMatrixVector (const Vector<number>* M=0);
00413
00419 PointerMatrixVector (const char* name);
00420
00430 PointerMatrixVector (const Vector<number>* M,
00431 const char* name);
00432
00433
00434 virtual void clear();
00435
00440 bool empty () const;
00441
00448 const PointerMatrixVector& operator= (const Vector<number>* M);
00449
00462 virtual void vmult (Vector<number>& dst,
00463 const Vector<number>& src) const;
00464
00478 virtual void Tvmult (Vector<number>& dst,
00479 const Vector<number>& src) const;
00480
00491 virtual void vmult_add (Vector<number>& dst,
00492 const Vector<number>& src) const;
00493
00504 virtual void Tvmult_add (Vector<number>& dst,
00505 const Vector<number>& src) const;
00506
00507 private:
00512 virtual const void* get() const;
00513
00517 SmartPointer<const Vector<number>,PointerMatrixVector<number> > m;
00518 };
00519
00520
00521
00541 template <class VECTOR, class MATRIX>
00542 inline
00543 PointerMatrixBase<VECTOR>*
00544 new_pointer_matrix_base(MATRIX& matrix, const VECTOR&, const char* name = "PointerMatrixAux")
00545 {
00546 return new PointerMatrixAux<MATRIX, VECTOR>(0, &matrix, name);
00547 }
00548
00555 template <typename numberv>
00556 PointerMatrixBase<Vector<numberv> >*
00557 new_pointer_matrix_base(const IdentityMatrix& matrix, const Vector<numberv>&, const char* name = "PointerMatrix")
00558 {
00559 return new PointerMatrix<IdentityMatrix, Vector<numberv> >(&matrix, name);
00560 }
00561
00562
00569 template <typename numberv, typename numberm>
00570 PointerMatrixBase<Vector<numberv> >*
00571 new_pointer_matrix_base(const FullMatrix<numberm>& matrix, const Vector<numberv>&, const char* name = "PointerMatrix")
00572 {
00573 return new PointerMatrix<FullMatrix<numberm>, Vector<numberv> >(&matrix, name);
00574 }
00575
00576
00583 template <typename numberv, typename numberm>
00584 PointerMatrixBase<Vector<numberv> >*
00585 new_pointer_matrix_base(const LAPACKFullMatrix<numberm>& matrix, const Vector<numberv>&, const char* name = "PointerMatrix")
00586 {
00587 return new PointerMatrix<LAPACKFullMatrix<numberm>, Vector<numberv> >(&matrix, name);
00588 }
00589
00590
00597 template <typename numberv, typename numberm>
00598 PointerMatrixBase<Vector<numberv> >*
00599 new_pointer_matrix_base(const SparseMatrix<numberm>& matrix, const Vector<numberv>&, const char* name = "PointerMatrix")
00600 {
00601 return new PointerMatrix<SparseMatrix<numberm>, Vector<numberv> >(&matrix, name);
00602 }
00603
00604
00611 template <class VECTOR, typename numberm>
00612 PointerMatrixBase<VECTOR>*
00613 new_pointer_matrix_base(const BlockSparseMatrix<numberm>& matrix, const VECTOR&, const char* name = "PointerMatrix")
00614 {
00615 return new PointerMatrix<BlockSparseMatrix<numberm>, VECTOR>(&matrix, name);
00616 }
00617
00618
00625 template <typename numberv, typename numberm>
00626 PointerMatrixBase<Vector<numberv> >*
00627 new_pointer_matrix_base(const SparseMatrixEZ<numberm>& matrix, const Vector<numberv>&, const char* name = "PointerMatrix")
00628 {
00629 return new PointerMatrix<SparseMatrixEZ<numberm>, Vector<numberv> >(&matrix, name);
00630 }
00631
00632
00639 template <class VECTOR, typename numberm>
00640 PointerMatrixBase<VECTOR>*
00641 new_pointer_matrix_base(const BlockSparseMatrixEZ<numberm>& matrix, const VECTOR&, const char* name = "PointerMatrix")
00642 {
00643 return new PointerMatrix<BlockSparseMatrixEZ<numberm>, VECTOR>(&matrix, name);
00644 }
00645
00646
00653 template <typename numberv, typename numberm>
00654 PointerMatrixBase<BlockVector<numberv> >*
00655 new_pointer_matrix_base(const BlockMatrixArray<numberm>& matrix, const BlockVector<numberv>&, const char* name = "PointerMatrix")
00656 {
00657 return new PointerMatrix<BlockMatrixArray<numberm>, BlockVector<numberv> >(&matrix, name);
00658 }
00659
00660
00667 template <typename numberv, typename numberm>
00668 PointerMatrixBase<Vector<numberv> >*
00669 new_pointer_matrix_base(const TridiagonalMatrix<numberm>& matrix, const Vector<numberv>&, const char* name = "PointerMatrix")
00670 {
00671 return new PointerMatrix<TridiagonalMatrix<numberm>, Vector<numberv> >(&matrix, name);
00672 }
00673
00674
00675
00677
00678
00679 template<class VECTOR>
00680 inline
00681 PointerMatrixBase<VECTOR>::~PointerMatrixBase ()
00682 {}
00683
00684
00685
00686 template<class VECTOR>
00687 inline
00688 bool
00689 PointerMatrixBase<VECTOR>::operator == (const PointerMatrixBase<VECTOR>& other) const
00690 {
00691 return (get() == other.get());
00692 }
00693
00694
00695
00696 template<class VECTOR>
00697 inline
00698 bool
00699 PointerMatrixBase<VECTOR>::operator != (const PointerMatrixBase<VECTOR>& other) const
00700 {
00701 return (get() != other.get());
00702 }
00703
00704
00705
00706 template<class VECTOR>
00707 inline
00708 bool
00709 PointerMatrixBase<VECTOR>::operator < (const PointerMatrixBase<VECTOR>& other) const
00710 {
00711 return (get() < other.get());
00712 }
00713
00714
00715
00716 template<class VECTOR>
00717 inline
00718 bool
00719 PointerMatrixBase<VECTOR>::operator <= (const PointerMatrixBase<VECTOR>& other) const
00720 {
00721 return (get() <= other.get());
00722 }
00723
00724
00725
00726 template<class VECTOR>
00727 inline
00728 bool
00729 PointerMatrixBase<VECTOR>::operator > (const PointerMatrixBase<VECTOR>& other) const
00730 {
00731 return (get() > other.get());
00732 }
00733
00734
00735
00736 template<class VECTOR>
00737 inline
00738 bool
00739 PointerMatrixBase<VECTOR>::operator >= (const PointerMatrixBase<VECTOR>& other) const
00740 {
00741 return (get() >= other.get());
00742 }
00743
00744
00745
00746
00747
00748 template<class MATRIX, class VECTOR>
00749 PointerMatrix<MATRIX, VECTOR>::PointerMatrix (const MATRIX* M)
00750 : m(M, typeid(*this).name())
00751 {}
00752
00753
00754 template<class MATRIX, class VECTOR>
00755 PointerMatrix<MATRIX, VECTOR>::PointerMatrix (const char* name)
00756 : m(0, name)
00757 {}
00758
00759
00760 template<class MATRIX, class VECTOR>
00761 PointerMatrix<MATRIX, VECTOR>::PointerMatrix (
00762 const MATRIX* M,
00763 const char* name)
00764 : m(M, name)
00765 {}
00766
00767
00768 template<class MATRIX, class VECTOR>
00769 inline void
00770 PointerMatrix<MATRIX, VECTOR>::clear ()
00771 {
00772 m = 0;
00773 }
00774
00775
00776 template<class MATRIX, class VECTOR>
00777 inline const PointerMatrix<MATRIX, VECTOR>&
00778 PointerMatrix<MATRIX, VECTOR>::operator= (const MATRIX* M)
00779 {
00780 m = M;
00781 return *this;
00782 }
00783
00784
00785 template<class MATRIX, class VECTOR>
00786 inline bool
00787 PointerMatrix<MATRIX, VECTOR>::empty () const
00788 {
00789 if (m == 0)
00790 return true;
00791 return m->empty();
00792 }
00793
00794 template<class MATRIX, class VECTOR>
00795 inline void
00796 PointerMatrix<MATRIX, VECTOR>::vmult (VECTOR& dst,
00797 const VECTOR& src) const
00798 {
00799 Assert (m != 0, ExcNotInitialized());
00800 m->vmult (dst, src);
00801 }
00802
00803
00804 template<class MATRIX, class VECTOR>
00805 inline void
00806 PointerMatrix<MATRIX, VECTOR>::Tvmult (VECTOR& dst,
00807 const VECTOR& src) const
00808 {
00809 Assert (m != 0, ExcNotInitialized());
00810 m->Tvmult (dst, src);
00811 }
00812
00813
00814 template<class MATRIX, class VECTOR>
00815 inline void
00816 PointerMatrix<MATRIX, VECTOR>::vmult_add (VECTOR& dst,
00817 const VECTOR& src) const
00818 {
00819 Assert (m != 0, ExcNotInitialized());
00820 m->vmult_add (dst, src);
00821 }
00822
00823
00824 template<class MATRIX, class VECTOR>
00825 inline void
00826 PointerMatrix<MATRIX, VECTOR>::Tvmult_add (VECTOR& dst,
00827 const VECTOR& src) const
00828 {
00829 Assert (m != 0, ExcNotInitialized());
00830 m->Tvmult_add (dst, src);
00831 }
00832
00833
00834 template<class MATRIX, class VECTOR>
00835 inline const void*
00836 PointerMatrix<MATRIX, VECTOR>::get () const
00837 {
00838 return m;
00839 }
00840
00841
00842
00843
00844
00845 template<class MATRIX, class VECTOR>
00846 PointerMatrixAux<MATRIX, VECTOR>::PointerMatrixAux (
00847 VectorMemory<VECTOR>* mem,
00848 const MATRIX* M)
00849 : mem(mem, typeid(*this).name()),
00850 m(M, typeid(*this).name())
00851 {
00852 if (mem == 0) mem = &my_memory;
00853 }
00854
00855
00856 template<class MATRIX, class VECTOR>
00857 PointerMatrixAux<MATRIX, VECTOR>::PointerMatrixAux (
00858 VectorMemory<VECTOR>* mem,
00859 const char* name)
00860 : mem(mem, name),
00861 m(0, name)
00862 {
00863 if (mem == 0) mem = &my_memory;
00864 }
00865
00866
00867 template<class MATRIX, class VECTOR>
00868 PointerMatrixAux<MATRIX, VECTOR>::PointerMatrixAux (
00869 VectorMemory<VECTOR>* mem,
00870 const MATRIX* M,
00871 const char* name)
00872 : mem(mem, name),
00873 m(M, name)
00874 {
00875 if (mem == 0) mem = &my_memory;
00876 }
00877
00878
00879 template<class MATRIX, class VECTOR>
00880 inline void
00881 PointerMatrixAux<MATRIX, VECTOR>::clear ()
00882 {
00883 m = 0;
00884 }
00885
00886
00887 template<class MATRIX, class VECTOR>
00888 inline const PointerMatrixAux<MATRIX, VECTOR>&
00889 PointerMatrixAux<MATRIX, VECTOR>::operator= (const MATRIX* M)
00890 {
00891 m = M;
00892 return *this;
00893 }
00894
00895
00896 template<class MATRIX, class VECTOR>
00897 inline void
00898 PointerMatrixAux<MATRIX, VECTOR>::set_memory(VectorMemory<VECTOR>* M)
00899 {
00900 mem = M;
00901 if (mem == 0)
00902 mem = &my_memory;
00903 }
00904
00905
00906 template<class MATRIX, class VECTOR>
00907 inline bool
00908 PointerMatrixAux<MATRIX, VECTOR>::empty () const
00909 {
00910 if (m == 0)
00911 return true;
00912 return m->empty();
00913 }
00914
00915 template<class MATRIX, class VECTOR>
00916 inline void
00917 PointerMatrixAux<MATRIX, VECTOR>::vmult (VECTOR& dst,
00918 const VECTOR& src) const
00919 {
00920 if (mem == 0)
00921 mem = &my_memory;
00922 Assert (mem != 0, ExcNotInitialized());
00923 Assert (m != 0, ExcNotInitialized());
00924 m->vmult (dst, src);
00925 }
00926
00927
00928 template<class MATRIX, class VECTOR>
00929 inline void
00930 PointerMatrixAux<MATRIX, VECTOR>::Tvmult (VECTOR& dst,
00931 const VECTOR& src) const
00932 {
00933 if (mem == 0)
00934 mem = &my_memory;
00935 Assert (mem != 0, ExcNotInitialized());
00936 Assert (m != 0, ExcNotInitialized());
00937 m->Tvmult (dst, src);
00938 }
00939
00940
00941 template<class MATRIX, class VECTOR>
00942 inline void
00943 PointerMatrixAux<MATRIX, VECTOR>::vmult_add (VECTOR& dst,
00944 const VECTOR& src) const
00945 {
00946 if (mem == 0)
00947 mem = &my_memory;
00948 Assert (mem != 0, ExcNotInitialized());
00949 Assert (m != 0, ExcNotInitialized());
00950 VECTOR* v = mem->alloc();
00951 v->reinit(dst);
00952 m->vmult (*v, src);
00953 dst += *v;
00954 mem->free(v);
00955 }
00956
00957
00958 template<class MATRIX, class VECTOR>
00959 inline void
00960 PointerMatrixAux<MATRIX, VECTOR>::Tvmult_add (VECTOR& dst,
00961 const VECTOR& src) const
00962 {
00963 if (mem == 0)
00964 mem = &my_memory;
00965 Assert (mem != 0, ExcNotInitialized());
00966 Assert (m != 0, ExcNotInitialized());
00967 VECTOR* v = mem->alloc();
00968 v->reinit(dst);
00969 m->Tvmult (*v, src);
00970 dst += *v;
00971 mem->free(v);
00972 }
00973
00974
00975 template<class MATRIX, class VECTOR>
00976 inline const void*
00977 PointerMatrixAux<MATRIX, VECTOR>::get () const
00978 {
00979 return m;
00980 }
00981
00982
00983
00984
00985
00986 template<typename number>
00987 PointerMatrixVector<number>::PointerMatrixVector (const Vector<number>* M)
00988 : m(M, typeid(*this).name())
00989 {}
00990
00991
00992 template<typename number>
00993 PointerMatrixVector<number>::PointerMatrixVector (const char* name)
00994 : m(0, name)
00995 {}
00996
00997
00998 template<typename number>
00999 PointerMatrixVector<number>::PointerMatrixVector (
01000 const Vector<number>* M,
01001 const char* name)
01002 : m(M, name)
01003 {}
01004
01005
01006 template<typename number>
01007 inline void
01008 PointerMatrixVector<number>::clear ()
01009 {
01010 m = 0;
01011 }
01012
01013
01014 template<typename number>
01015 inline const PointerMatrixVector<number>&
01016 PointerMatrixVector<number>::operator= (const Vector<number>* M)
01017 {
01018 m = M;
01019 return *this;
01020 }
01021
01022
01023 template<typename number>
01024 inline bool
01025 PointerMatrixVector<number>::empty () const
01026 {
01027 if (m == 0)
01028 return true;
01029 return m->empty();
01030 }
01031
01032 template<typename number>
01033 inline void
01034 PointerMatrixVector<number>::vmult (
01035 Vector<number>& dst,
01036 const Vector<number>& src) const
01037 {
01038 Assert (m != 0, ExcNotInitialized());
01039 Assert (dst.size() == 1, ExcDimensionMismatch(dst.size(), 1));
01040
01041 dst(0) = *m * src;
01042 }
01043
01044
01045 template<typename number>
01046 inline void
01047 PointerMatrixVector<number>::Tvmult (
01048 Vector<number>& dst,
01049 const Vector<number>& src) const
01050 {
01051 Assert (m != 0, ExcNotInitialized());
01052 Assert(src.size() == 1, ExcDimensionMismatch(src.size(), 1));
01053
01054 dst.equ (src(0), *m);
01055 }
01056
01057
01058 template<typename number>
01059 inline void
01060 PointerMatrixVector<number>::vmult_add (
01061 Vector<number>& dst,
01062 const Vector<number>& src) const
01063 {
01064 Assert (m != 0, ExcNotInitialized());
01065 Assert (dst.size() == 1, ExcDimensionMismatch(dst.size(), 1));
01066
01067 dst(0) += *m * src;
01068 }
01069
01070
01071 template<typename number>
01072 inline void
01073 PointerMatrixVector<number>::Tvmult_add (
01074 Vector<number>& dst,
01075 const Vector<number>& src) const
01076 {
01077 Assert (m != 0, ExcNotInitialized());
01078 Assert(src.size() == 1, ExcDimensionMismatch(src.size(), 1));
01079
01080 dst.add (src(0), *m);
01081 }
01082
01083
01084 template<typename number>
01085 inline const void*
01086 PointerMatrixVector<number>::get () const
01087 {
01088 return m;
01089 }
01090
01091
01092
01093 DEAL_II_NAMESPACE_CLOSE
01094
01095 #endif
01096