include/deal.II/lac/matrix_lib.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: matrix_lib.h 25526 2012-05-21 12:08:58Z bangerth @f$
00003 //
00004 //    Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2009, 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__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                                      // Doc in PointerMatrixBase
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                                      // Doc in PointerMatrixBase
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:11 by doxygen 1.7.3