Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
matrix_lib.h
1 // ---------------------------------------------------------------------
2 // @f$Id: matrix_lib.h 30036 2013-07-18 16:55:32Z maier @f$
3 //
4 // Copyright (C) 2002 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__matrix_lib_h
18 #define __deal2__matrix_lib_h
19 
20 #include <deal.II/base/subscriptor.h>
21 #include <deal.II/lac/vector_memory.h>
22 #include <deal.II/lac/pointer_matrix.h>
23 #include <deal.II/lac/solver_richardson.h>
24 
26 
27 template<typename number> class Vector;
28 template<typename number> class BlockVector;
29 template<typename number> class SparseMatrix;
30 
48 template<class VECTOR>
49 class ProductMatrix : public PointerMatrixBase<VECTOR>
50 {
51 public:
58  ProductMatrix();
59 
66 
73  template <class MATRIX1, class MATRIX2>
74  ProductMatrix(const MATRIX1 &m1,
75  const MATRIX2 &m2,
77 
81  template <class MATRIX1, class MATRIX2>
82  void reinit(const MATRIX1 &m1, const MATRIX2 &m2);
83 
87  template <class MATRIX1, class MATRIX2>
88  void initialize(const MATRIX1 &m1, const MATRIX2 &m2,
90 
95 
96  // Doc in PointerMatrixBase
97  void clear();
98 
103  virtual void vmult (VECTOR &w,
104  const VECTOR &v) const;
105 
111  virtual void Tvmult (VECTOR &w,
112  const VECTOR &v) const;
113 
118  virtual void vmult_add (VECTOR &w,
119  const VECTOR &v) const;
120 
127  virtual void Tvmult_add (VECTOR &w,
128  const VECTOR &v) const;
129 
130 private:
135 
140 
149  virtual const void *get() const;
150 };
151 
152 
161 template<class VECTOR>
162 class ScaledMatrix : public Subscriptor
163 {
164 public:
169  ScaledMatrix ();
173  template <class MATRIX>
174  ScaledMatrix (const MATRIX &M, const double factor);
175 
179  ~ScaledMatrix ();
184  template <class MATRIX>
185  void initialize (const MATRIX &M, const double factor);
186 
190  void clear ();
191 
195  void vmult (VECTOR &w, const VECTOR &v) const;
196 
201  void Tvmult (VECTOR &w, const VECTOR &v) const;
202 
203 private:
211  double factor;
212 };
213 
214 
215 
228 template<typename number, typename vector_number>
229 class ProductSparseMatrix : public PointerMatrixBase<Vector<vector_number> >
230 {
231 public:
236 
242 
249  ProductSparseMatrix(const MatrixType &m1,
250  const MatrixType &m2,
252 
261 
262  void initialize(const MatrixType &m1,
263  const MatrixType &m2,
265 
266  // Doc in PointerMatrixBase
267  void clear();
268 
273  virtual void vmult (VectorType &w,
274  const VectorType &v) const;
275 
281  virtual void Tvmult (VectorType &w,
282  const VectorType &v) const;
283 
288  virtual void vmult_add (VectorType &w,
289  const VectorType &v) const;
290 
297  virtual void Tvmult_add (VectorType &w,
298  const VectorType &v) const;
299 
300 private:
305 
310 
319  virtual const void *get() const;
320 };
321 
322 
332 {
333 public:
338 
344 
348  template <typename number>
349  void filter (Vector<number> &v) const;
350 
354  template <typename number>
355  void filter (BlockVector<number> &v) const;
356 
361  template <typename number>
362  void vmult (Vector<number> &dst,
363  const Vector<number> &src) const;
364 
369  template <typename number>
370  void vmult_add (Vector<number> &dst,
371  const Vector<number> &src) const;
372 
378  template <typename number>
379  void vmult (BlockVector<number> &dst,
380  const BlockVector<number> &src) const;
381 
387  template <typename number>
388  void vmult_add (BlockVector<number> &dst,
389  const BlockVector<number> &src) const;
390 
391 
395  template <typename VECTOR>
396  void Tvmult(VECTOR &, const VECTOR &) const;
397 
401  template <typename VECTOR>
402  void Tvmult_add(VECTOR &, const VECTOR &) const;
403 
404 private:
408  size_type component;
409 };
410 
411 
412 
430 template<class VECTOR>
432 {
433 public:
448 
455  template <class MATRIX, class PRECONDITION>
456  void initialize (const MATRIX &,
457  const PRECONDITION &);
458 
463  SolverControl &control() const;
467  void vmult (VECTOR &, const VECTOR &) const;
468 
472  void vmult_add (VECTOR &, const VECTOR &) const;
473 
477  void Tvmult (VECTOR &, const VECTOR &) const;
478 
482  void Tvmult_add (VECTOR &, const VECTOR &) const;
483 
484 private:
490 
495 
500 
505 };
506 
507 
508 
509 
511 //---------------------------------------------------------------------------
512 
513 
514 template<class VECTOR>
515 inline
517  :
518  m(0)
519 {}
520 
521 
522 
523 template<class VECTOR>
524 template<class MATRIX>
525 inline
526 ScaledMatrix<VECTOR>::ScaledMatrix(const MATRIX &mat, const double factor)
527  :
528  m(new_pointer_matrix_base(mat, VECTOR())),
529  factor(factor)
530 {}
531 
532 
533 
534 template<class VECTOR>
535 template<class MATRIX>
536 inline
537 void
538 ScaledMatrix<VECTOR>::initialize(const MATRIX &mat, const double f)
539 {
540  if (m) delete m;
541  m = new_pointer_matrix_base(mat, VECTOR());
542  factor = f;
543 }
544 
545 
546 
547 template<class VECTOR>
548 inline
549 void
551 {
552  if (m) delete m;
553  m = 0;
554 }
555 
556 
557 
558 template<class VECTOR>
559 inline
561 {
562  clear ();
563 }
564 
565 
566 template<class VECTOR>
567 inline
568 void
569 ScaledMatrix<VECTOR>::vmult (VECTOR &w, const VECTOR &v) const
570 {
571  m->vmult(w, v);
572  w *= factor;
573 }
574 
575 
576 template<class VECTOR>
577 inline
578 void
579 ScaledMatrix<VECTOR>::Tvmult (VECTOR &w, const VECTOR &v) const
580 {
581  m->Tvmult(w, v);
582  w *= factor;
583 }
584 
585 
586 //---------------------------------------------------------------------------
587 
588 template<class VECTOR>
590  : m1(0), m2(0), mem(0)
591 {}
592 
593 
594 template<class VECTOR>
596  : m1(0), m2(0), mem(&m)
597 {}
598 
599 
600 template<class VECTOR>
601 template<class MATRIX1, class MATRIX2>
603  const MATRIX1 &mat1,
604  const MATRIX2 &mat2,
606  : mem(&m)
607 {
608  m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1, typeid(*this).name());
609  m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2, typeid(*this).name());
610 }
611 
612 
613 template<class VECTOR>
614 template<class MATRIX1, class MATRIX2>
615 void
617  const MATRIX1 &mat1,
618  const MATRIX2 &mat2)
619 {
620  if (m1) delete m1;
621  if (m2) delete m2;
622  m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1, typeid(*this).name());
623  m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2, typeid(*this).name());
624 }
625 
626 
627 template<class VECTOR>
628 template<class MATRIX1, class MATRIX2>
629 void
631  const MATRIX1 &mat1,
632  const MATRIX2 &mat2,
633  VectorMemory<VECTOR> &memory)
634 {
635  mem = &memory;
636  if (m1) delete m1;
637  if (m2) delete m2;
638  m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1, typeid(*this).name());
639  m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2, typeid(*this).name());
640 }
641 
642 
643 template<class VECTOR>
645 {
646  if (m1) delete m1;
647  if (m2) delete m2;
648 }
649 
650 
651 template<class VECTOR>
652 void
654 {
655  if (m1) delete m1;
656  m1 = 0;
657  if (m2) delete m2;
658  m2 = 0;
659 }
660 
661 
662 template<class VECTOR>
663 void
664 ProductMatrix<VECTOR>::vmult (VECTOR &dst, const VECTOR &src) const
665 {
666  Assert (mem != 0, ExcNotInitialized());
667  Assert (m1 != 0, ExcNotInitialized());
668  Assert (m2 != 0, ExcNotInitialized());
669 
670  VECTOR *v = mem->alloc();
671  v->reinit(dst);
672  m2->vmult (*v, src);
673  m1->vmult (dst, *v);
674  mem->free(v);
675 }
676 
677 
678 template<class VECTOR>
679 void
680 ProductMatrix<VECTOR>::vmult_add (VECTOR &dst, const VECTOR &src) const
681 {
682  Assert (mem != 0, ExcNotInitialized());
683  Assert (m1 != 0, ExcNotInitialized());
684  Assert (m2 != 0, ExcNotInitialized());
685 
686  VECTOR *v = mem->alloc();
687  v->reinit(dst);
688  m2->vmult (*v, src);
689  m1->vmult_add (dst, *v);
690  mem->free(v);
691 }
692 
693 
694 template<class VECTOR>
695 void
696 ProductMatrix<VECTOR>::Tvmult (VECTOR &dst, const VECTOR &src) const
697 {
698  Assert (mem != 0, ExcNotInitialized());
699  Assert (m1 != 0, ExcNotInitialized());
700  Assert (m2 != 0, ExcNotInitialized());
701 
702  VECTOR *v = mem->alloc();
703  v->reinit(dst);
704  m1->Tvmult (*v, src);
705  m2->Tvmult (dst, *v);
706  mem->free(v);
707 }
708 
709 
710 template<class VECTOR>
711 void
712 ProductMatrix<VECTOR>::Tvmult_add (VECTOR &dst, const VECTOR &src) const
713 {
714  Assert (mem != 0, ExcNotInitialized());
715  Assert (m1 != 0, ExcNotInitialized());
716  Assert (m2 != 0, ExcNotInitialized());
717 
718  VECTOR *v = mem->alloc();
719  v->reinit(dst);
720  m1->Tvmult (*v, src);
721  m2->Tvmult_add (dst, *v);
722  mem->free(v);
723 }
724 
725 
726 template<class VECTOR>
727 const void *
729 {
730  return (void *) m1;
731 }
732 
733 
734 //---------------------------------------------------------------------------
735 
736 template <class VECTOR>
737 inline void
738 MeanValueFilter::Tvmult(VECTOR &, const VECTOR &) const
739 {
740  Assert(false, ExcNotImplemented());
741 }
742 
743 
744 template <class VECTOR>
745 inline void
746 MeanValueFilter::Tvmult_add(VECTOR &, const VECTOR &) const
747 {
748  Assert(false, ExcNotImplemented());
749 }
750 
751 //-----------------------------------------------------------------------//
752 
753 template <class VECTOR>
754 template <class MATRIX, class PRECONDITION>
755 inline void
756 InverseMatrixRichardson<VECTOR>::initialize (const MATRIX &m, const PRECONDITION &p)
757 {
758  if (matrix != 0)
759  delete matrix;
760  matrix = new PointerMatrix<MATRIX, VECTOR>(&m);
761  if (precondition != 0)
762  delete precondition;
763  precondition = new PointerMatrix<PRECONDITION, VECTOR>(&p);
764 }
765 
766 
767 DEAL_II_NAMESPACE_CLOSE
768 
769 #endif
void vmult(VECTOR &, const VECTOR &) const
const types::global_dof_index invalid_size_type
Definition: types.h:211
virtual const void * get() const
Definition: matrix_lib.h:728
PointerMatrixBase< VECTOR > * new_pointer_matrix_base(MATRIX &matrix, const VECTOR &, const char *name="PointerMatrixAux")
virtual void vmult(VectorType &w, const VectorType &v) const
void vmult(Vector< number > &dst, const Vector< number > &src) const
virtual void vmult(VECTOR &w, const VECTOR &v) const
Definition: matrix_lib.h:664
PointerMatrixBase< VECTOR > * m
Definition: matrix_lib.h:207
void filter(Vector< number > &v) const
virtual void Tvmult_add(VectorType &w, const VectorType &v) const
void reinit(const MATRIX1 &m1, const MATRIX2 &m2)
Definition: matrix_lib.h:616
VectorMemory< VECTOR > & mem
Definition: matrix_lib.h:489
void vmult(VECTOR &w, const VECTOR &v) const
Definition: matrix_lib.h:569
SmartPointer< VectorMemory< VectorType >, ProductSparseMatrix< number, vector_number > > mem
Definition: matrix_lib.h:314
MeanValueFilter(size_type component=numbers::invalid_size_type)
void Tvmult(VECTOR &w, const VECTOR &v) const
Definition: matrix_lib.h:579
void initialize(const MATRIX &M, const double factor)
Definition: matrix_lib.h:538
virtual void vmult_add(VECTOR &w, const VECTOR &v) const
Definition: matrix_lib.h:680
virtual void Tvmult_add(VECTOR &w, const VECTOR &v) const
Definition: matrix_lib.h:712
void initialize(const MATRIX1 &m1, const MATRIX2 &m2, VectorMemory< VECTOR > &mem)
Definition: matrix_lib.h:630
PointerMatrixBase< VECTOR > * precondition
Definition: matrix_lib.h:504
unsigned int global_dof_index
Definition: types.h:100
Vector< vector_number > VectorType
Definition: matrix_lib.h:241
SmartPointer< const MatrixType, ProductSparseMatrix< number, vector_number > > m1
Definition: matrix_lib.h:304
#define Assert(cond, exc)
Definition: exceptions.h:299
void Tvmult_add(VECTOR &, const VECTOR &) const
PointerMatrixBase< VECTOR > * m2
Definition: matrix_lib.h:139
SolverRichardson< VECTOR > solver
Definition: matrix_lib.h:494
SparseMatrix< number > MatrixType
Definition: matrix_lib.h:235
SmartPointer< VectorMemory< VECTOR >, ProductMatrix< VECTOR > > mem
Definition: matrix_lib.h:144
SmartPointer< const MatrixType, ProductSparseMatrix< number, vector_number > > m2
Definition: matrix_lib.h:309
virtual void Tvmult(VectorType &w, const VectorType &v) const
InverseMatrixRichardson(SolverControl &control, VectorMemory< VECTOR > &mem)
types::global_dof_index size_type
Definition: matrix_lib.h:337
virtual void vmult_add(VectorType &w, const VectorType &v) const
void Tvmult(VECTOR &, const VECTOR &) const
Definition: matrix_lib.h:738
void Tvmult_add(VECTOR &, const VECTOR &) const
Definition: matrix_lib.h:746
void vmult_add(VECTOR &, const VECTOR &) const
void initialize(const MATRIX &, const PRECONDITION &)
Definition: matrix_lib.h:756
void clear()
Definition: matrix_lib.h:653
SolverControl & control() const
PointerMatrixBase< VECTOR > * matrix
Definition: matrix_lib.h:499
void clear()
Definition: matrix_lib.h:550
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
PointerMatrixBase< VECTOR > * m1
Definition: matrix_lib.h:134
void vmult_add(Vector< number > &dst, const Vector< number > &src) const
size_type component
Definition: matrix_lib.h:408
double factor
Definition: matrix_lib.h:211
virtual void Tvmult(VECTOR &w, const VECTOR &v) const
Definition: matrix_lib.h:696
void Tvmult(VECTOR &, const VECTOR &) const