include/deal.II/lac/pointer_matrix.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: pointer_matrix.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                      // Use doc from base class
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                                      // Use doc from base class
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                                      // Use doc from base class
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3