Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
lapack_full_matrix.h
1 // ---------------------------------------------------------------------
2 // @f$Id: lapack_full_matrix.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 2005 - 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__lapack_full_matrix_h
18 #define __deal2__lapack_full_matrix_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/base/table.h>
24 #include <deal.II/lac/lapack_support.h>
25 #include <deal.II/lac/vector_memory.h>
26 
27 #include <deal.II/base/std_cxx1x/shared_ptr.h>
28 #include <vector>
29 #include <complex>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 // forward declarations
34 template<typename number> class Vector;
35 template<typename number> class BlockVector;
36 template<typename number> class FullMatrix;
37 
38 
53 template <typename number>
54 class LAPACKFullMatrix : public TransposeTable<number>
55 {
56 public:
61 
76  explicit LAPACKFullMatrix (const size_type n = 0);
77 
83  LAPACKFullMatrix (const size_type rows,
84  const size_type cols);
85 
106 
112 
120  template <typename number2>
123 
132  operator = (const double d);
133 
144  template <class MATRIX>
145  void copy_from (const MATRIX &);
146 
169  template<class MATRIX>
170  void fill (const MATRIX &src,
171  const size_type dst_offset_i = 0,
172  const size_type dst_offset_j = 0,
173  const size_type src_offset_i = 0,
174  const size_type src_offset_j = 0,
175  const number factor = 1.,
176  const bool transpose = false);
177 
178 
225  template <class VECTOR>
226  void vmult(VECTOR &dst, const VECTOR &src, const bool adding = false) const;
235  template <class VECTOR>
236  void vmult_add (VECTOR &w, const VECTOR &v) const;
237 
257  template <class VECTOR>
258  void Tvmult (VECTOR &w, const VECTOR &v,
259  const bool adding=false) const;
260 
270  template <class VECTOR>
271  void Tvmult_add (VECTOR &w, const VECTOR &v) const;
272 
273  void vmult (Vector<number> &w,
274  const Vector<number> &v,
275  const bool adding=false) const;
276  void vmult_add (Vector<number> &w,
277  const Vector<number> &v) const;
278  void Tvmult (Vector<number> &w,
279  const Vector<number> &v,
280  const bool adding=false) const;
281  void Tvmult_add (Vector<number> &w,
282  const Vector<number> &v) const;
288  void compute_lu_factorization ();
289 
296  void invert ();
297 
308  const bool transposed) const;
309 
346  void compute_eigenvalues (const bool right_eigenvectors = false,
347  const bool left_eigenvectors = false);
348 
382  const number lower_bound,
383  const number upper_bound,
384  const number abs_accuracy,
385  Vector<number> &eigenvalues,
386  FullMatrix<number> &eigenvectors);
387 
431  const number lower_bound,
432  const number upper_bound,
433  const number abs_accuracy,
434  Vector<number> &eigenvalues,
435  std::vector<Vector<number> > &eigenvectors,
436  const int itype = 1);
437 
476  std::vector<Vector<number> > &eigenvectors,
477  const int itype = 1);
478 
492  void compute_svd();
531  void compute_inverse_svd (const double threshold = 0.);
532 
538  std::complex<number>
539  eigenvalue (const size_type i) const;
540 
547  number
548  singular_value (const size_type i) const;
549 
591  void print_formatted (std::ostream &out,
592  const unsigned int presicion=3,
593  const bool scientific = true,
594  const unsigned int width = 0,
595  const char *zero_string = " ",
596  const double denominator = 1.,
597  const double threshold = 0.) const;
598 
599 private:
607  LAPACKSupport::State state;
614  LAPACKSupport::Properties properties;
615 
620  mutable std::vector<number> work;
621 
632  std::vector<int> ipiv;
633 
639  std::vector<number> inv_work;
640 
646  std::vector<number> wr;
647 
653  std::vector<number> wi;
654 
659  std::vector<number> vl;
660 
665  std::vector<number> vr;
666 
672  std_cxx1x::shared_ptr<LAPACKFullMatrix<number> > svd_u;
679  std_cxx1x::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
680 };
681 
682 
683 
690 template <typename number>
692  :
693  public Subscriptor
694 {
695 public:
696  void initialize(const LAPACKFullMatrix<number> &);
697  void initialize(const LAPACKFullMatrix<number> &,
699  void vmult(Vector<number> &, const Vector<number> &) const;
700  void Tvmult(Vector<number> &, const Vector<number> &) const;
701  void vmult(BlockVector<number> &,
702  const BlockVector<number> &) const;
703  void Tvmult(BlockVector<number> &,
704  const BlockVector<number> &) const;
705 private:
708 };
709 
710 
711 
712 template <typename number>
713 template <class MATRIX>
714 inline void
716 {
717  this->reinit (M.m(), M.n());
718 
719  // loop over the elements of the argument matrix row by row, as suggested
720  // in the documentation of the sparse matrix iterator class, and
721  // copy them into the current object
722  for (size_type row = 0; row < M.m(); ++row)
723  {
724  const typename MATRIX::const_iterator end_row = M.end(row);
725  for (typename MATRIX::const_iterator entry = M.begin(row);
726  entry != end_row; ++entry)
727  this->el(row, entry->column()) = entry->value();
728  }
729 
730  state = LAPACKSupport::matrix;
731 }
732 
733 
734 
735 template <typename number>
736 template <class MATRIX>
737 inline void
739  const MATRIX &M,
740  const size_type dst_offset_i,
741  const size_type dst_offset_j,
742  const size_type src_offset_i,
743  const size_type src_offset_j,
744  const number factor,
745  const bool transpose)
746 {
747  // loop over the elements of the argument matrix row by row, as suggested
748  // in the documentation of the sparse matrix iterator class
749  for (size_type row = src_offset_i; row < M.m(); ++row)
750  {
751  const typename MATRIX::const_iterator end_row = M.end(row);
752  for (typename MATRIX::const_iterator entry = M.begin(row);
753  entry != end_row; ++entry)
754  {
755  const size_type i = transpose ? entry->column() : row;
756  const size_type j = transpose ? row : entry->column();
757 
758  const size_type dst_i=dst_offset_i+i-src_offset_i;
759  const size_type dst_j=dst_offset_j+j-src_offset_j;
760  if (dst_i<this->n_rows() && dst_j<this->n_cols())
761  (*this)(dst_i, dst_j) = factor * entry->value();
762  }
763  }
764 
765  state = LAPACKSupport::matrix;
766 }
767 
768 
769 template <typename number>
770 template <class VECTOR>
771 inline void
772 LAPACKFullMatrix<number>::vmult(VECTOR &, const VECTOR &, bool) const
773 {
774  Assert(false, ExcNotImplemented());
775 }
776 
777 
778 template <typename number>
779 template <class VECTOR>
780 inline void
781 LAPACKFullMatrix<number>::vmult_add(VECTOR &, const VECTOR &) const
782 {
783  Assert(false, ExcNotImplemented());
784 }
785 
786 
787 template <typename number>
788 template <class VECTOR>
789 inline void
790 LAPACKFullMatrix<number>::Tvmult(VECTOR &, const VECTOR &, bool) const
791 {
792  Assert(false, ExcNotImplemented());
793 }
794 
795 
796 template <typename number>
797 template <class VECTOR>
798 inline void
799 LAPACKFullMatrix<number>::Tvmult_add(VECTOR &, const VECTOR &) const
800 {
801  Assert(false, ExcNotImplemented());
802 }
803 
804 
805 template <typename number>
806 inline std::complex<number>
808 {
809  Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState());
810  Assert (wr.size() == this->n_rows(), ExcInternalError());
811  Assert (wi.size() == this->n_rows(), ExcInternalError());
812  Assert (i<this->n_rows(), ExcIndexRange(i,0,this->n_rows()));
813 
814  return std::complex<number>(wr[i], wi[i]);
815 }
816 
817 
818 template <typename number>
819 inline number
821 {
822  Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state));
823  AssertIndexRange(i,wr.size());
824 
825  return wr[i];
826 }
827 
828 
829 
830 DEAL_II_NAMESPACE_CLOSE
831 
832 #endif
std::vector< number > work
number singular_value(const size_type i) const
void Tvmult_add(VECTOR &w, const VECTOR &v) const
std_cxx1x::shared_ptr< LAPACKFullMatrix< number > > svd_u
LAPACKSupport::State state
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
LAPACKFullMatrix(const size_type n=0)
LAPACKSupport::Properties properties
std_cxx1x::shared_ptr< LAPACKFullMatrix< number > > svd_vt
std::vector< number > vr
void compute_eigenvalues_symmetric(const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, FullMatrix< number > &eigenvectors)
void apply_lu_factorization(Vector< number > &v, const bool transposed) const
::ExceptionBase & ExcInvalidState()
types::global_dof_index size_type
std::vector< int > ipiv
void vmult_add(VECTOR &w, const VECTOR &v) const
std::vector< number > wr
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
void print_formatted(std::ostream &out, const unsigned int presicion=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1., const double threshold=0.) const
LAPACKFullMatrix< number > & operator=(const LAPACKFullMatrix< number > &)
SymmetricTensor< rank, dim, Number > transpose(const SymmetricTensor< rank, dim, Number > &t)
void vmult(VECTOR &dst, const VECTOR &src, const bool adding=false) const
std::vector< number > inv_work
void fill(const MATRIX &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0, const number factor=1., const bool transpose=false)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::vector< number > vl
void compute_lu_factorization()
std::vector< number > wi
void compute_generalized_eigenvalues_symmetric(LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number > > &eigenvectors, const int itype=1)
void copy_from(const MATRIX &)
void compute_inverse_svd(const double threshold=0.)
void compute_eigenvalues(const bool right_eigenvectors=false, const bool left_eigenvectors=false)
::ExceptionBase & ExcNotImplemented()
std::complex< number > eigenvalue(const size_type i) const
void Tvmult(VECTOR &w, const VECTOR &v, const bool adding=false) const
::ExceptionBase & ExcInternalError()