00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__lapack_full_matrix_h
00013 #define __deal2__lapack_full_matrix_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/smartpointer.h>
00018 #include <deal.II/base/table.h>
00019 #include <deal.II/lac/lapack_support.h>
00020 #include <deal.II/lac/vector_memory.h>
00021
00022 #include <deal.II/base/std_cxx1x/shared_ptr.h>
00023 #include <vector>
00024 #include <complex>
00025
00026 DEAL_II_NAMESPACE_OPEN
00027
00028
00029 template<typename number> class Vector;
00030 template<typename number> class BlockVector;
00031 template<typename number> class FullMatrix;
00032
00033
00048 template <typename number>
00049 class LAPACKFullMatrix : public TransposeTable<number>
00050 {
00051 public:
00066 explicit LAPACKFullMatrix (const unsigned int n = 0);
00067
00073 LAPACKFullMatrix (const unsigned int rows,
00074 const unsigned int cols);
00075
00095 LAPACKFullMatrix (const LAPACKFullMatrix&);
00096
00100 LAPACKFullMatrix<number> &
00101 operator = (const LAPACKFullMatrix<number>&);
00102
00110 template <typename number2>
00111 LAPACKFullMatrix<number> &
00112 operator = (const FullMatrix<number2>&);
00113
00121 LAPACKFullMatrix<number> &
00122 operator = (const double d);
00123
00134 template <class MATRIX>
00135 void copy_from (const MATRIX&);
00136
00159 template<class MATRIX>
00160 void fill (const MATRIX &src,
00161 const unsigned int dst_offset_i = 0,
00162 const unsigned int dst_offset_j = 0,
00163 const unsigned int src_offset_i = 0,
00164 const unsigned int src_offset_j = 0,
00165 const number factor = 1.,
00166 const bool transpose = false);
00167
00168
00215 template <class VECTOR>
00216 void vmult(VECTOR& dst, const VECTOR& src, const bool adding = false) const;
00225 template <class VECTOR>
00226 void vmult_add (VECTOR& w, const VECTOR& v) const;
00227
00247 template <class VECTOR>
00248 void Tvmult (VECTOR& w, const VECTOR& v,
00249 const bool adding=false) const;
00250
00260 template <class VECTOR>
00261 void Tvmult_add (VECTOR& w, const VECTOR& v) const;
00262
00263 void vmult (Vector<number> &w,
00264 const Vector<number> &v,
00265 const bool adding=false) const;
00266 void vmult_add (Vector<number> &w,
00267 const Vector<number> &v) const;
00268 void Tvmult (Vector<number> &w,
00269 const Vector<number> &v,
00270 const bool adding=false) const;
00271 void Tvmult_add (Vector<number> &w,
00272 const Vector<number> &v) const;
00278 void compute_lu_factorization ();
00279
00286 void invert ();
00287
00297 void apply_lu_factorization (Vector<number>& v,
00298 const bool transposed) const;
00299
00336 void compute_eigenvalues (const bool right_eigenvectors = false,
00337 const bool left_eigenvectors = false);
00338
00371 void compute_eigenvalues_symmetric(
00372 const number lower_bound,
00373 const number upper_bound,
00374 const number abs_accuracy,
00375 Vector<number> & eigenvalues,
00376 FullMatrix<number> & eigenvectors);
00377
00419 void compute_generalized_eigenvalues_symmetric(
00420 LAPACKFullMatrix<number> & B,
00421 const number lower_bound,
00422 const number upper_bound,
00423 const number abs_accuracy,
00424 Vector<number> & eigenvalues,
00425 std::vector<Vector<number> > & eigenvectors,
00426 const int itype = 1);
00427
00464 void compute_generalized_eigenvalues_symmetric (
00465 LAPACKFullMatrix<number> & B,
00466 std::vector<Vector<number> > & eigenvectors,
00467 const int itype = 1);
00468
00482 void compute_svd();
00521 void compute_inverse_svd (const double threshold = 0.);
00522
00528 std::complex<number>
00529 eigenvalue (const unsigned int i) const;
00530
00537 number
00538 singular_value (const unsigned int i) const;
00539
00581 void print_formatted (std::ostream &out,
00582 const unsigned int presicion=3,
00583 const bool scientific = true,
00584 const unsigned int width = 0,
00585 const char *zero_string = " ",
00586 const double denominator = 1.,
00587 const double threshold = 0.) const;
00588
00589 private:
00597 LAPACKSupport::State state;
00604 LAPACKSupport::Properties properties;
00605
00610 mutable std::vector<number> work;
00611
00622 std::vector<int> ipiv;
00623
00629 std::vector<number> inv_work;
00630
00636 std::vector<number> wr;
00637
00643 std::vector<number> wi;
00644
00649 std::vector<number> vl;
00650
00655 std::vector<number> vr;
00656
00662 std_cxx1x::shared_ptr<LAPACKFullMatrix<number> > svd_u;
00669 std_cxx1x::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
00670 };
00671
00672
00673
00680 template <typename number>
00681 class PreconditionLU
00682 :
00683 public Subscriptor
00684 {
00685 public:
00686 void initialize(const LAPACKFullMatrix<number>&);
00687 void initialize(const LAPACKFullMatrix<number>&,
00688 VectorMemory<Vector<number> >&);
00689 void vmult(Vector<number>&, const Vector<number>&) const;
00690 void Tvmult(Vector<number>&, const Vector<number>&) const;
00691 void vmult(BlockVector<number>&,
00692 const BlockVector<number>&) const;
00693 void Tvmult(BlockVector<number>&,
00694 const BlockVector<number>&) const;
00695 private:
00696 SmartPointer<const LAPACKFullMatrix<number>,PreconditionLU<number> > matrix;
00697 SmartPointer<VectorMemory<Vector<number> >,PreconditionLU<number> > mem;
00698 };
00699
00700
00701
00702 template <typename number>
00703 template <class MATRIX>
00704 inline void
00705 LAPACKFullMatrix<number>::copy_from (const MATRIX& M)
00706 {
00707 this->reinit (M.m(), M.n());
00708 const typename MATRIX::const_iterator end = M.end();
00709 for (typename MATRIX::const_iterator entry = M.begin();
00710 entry != end; ++entry)
00711 this->el(entry->row(), entry->column()) = entry->value();
00712
00713 state = LAPACKSupport::matrix;
00714 }
00715
00716
00717
00718 template <typename number>
00719 template <class MATRIX>
00720 inline void
00721 LAPACKFullMatrix<number>::fill (
00722 const MATRIX& M,
00723 const unsigned int dst_offset_i,
00724 const unsigned int dst_offset_j,
00725 const unsigned int src_offset_i,
00726 const unsigned int src_offset_j,
00727 const number factor,
00728 const bool transpose)
00729 {
00730 const typename MATRIX::const_iterator end = M.end();
00731 for (typename MATRIX::const_iterator entry = M.begin(src_offset_i);
00732 entry != end; ++entry)
00733 {
00734 const unsigned int i = transpose ? entry->column() : entry->row();
00735 const unsigned int j = transpose ? entry->row() : entry->column();
00736
00737 const unsigned int dst_i=dst_offset_i+i-src_offset_i;
00738 const unsigned int dst_j=dst_offset_j+j-src_offset_j;
00739 if (dst_i<this->n_rows() && dst_j<this->n_cols())
00740 (*this)(dst_i, dst_j) = factor * entry->value();
00741 }
00742
00743 state = LAPACKSupport::matrix;
00744 }
00745
00746
00747 template <typename number>
00748 template <class VECTOR>
00749 inline void
00750 LAPACKFullMatrix<number>::vmult(VECTOR&, const VECTOR&, bool) const
00751 {
00752 Assert(false, ExcNotImplemented());
00753 }
00754
00755
00756 template <typename number>
00757 template <class VECTOR>
00758 inline void
00759 LAPACKFullMatrix<number>::vmult_add(VECTOR&, const VECTOR&) const
00760 {
00761 Assert(false, ExcNotImplemented());
00762 }
00763
00764
00765 template <typename number>
00766 template <class VECTOR>
00767 inline void
00768 LAPACKFullMatrix<number>::Tvmult(VECTOR&, const VECTOR&, bool) const
00769 {
00770 Assert(false, ExcNotImplemented());
00771 }
00772
00773
00774 template <typename number>
00775 template <class VECTOR>
00776 inline void
00777 LAPACKFullMatrix<number>::Tvmult_add(VECTOR&, const VECTOR&) const
00778 {
00779 Assert(false, ExcNotImplemented());
00780 }
00781
00782
00783 template <typename number>
00784 inline std::complex<number>
00785 LAPACKFullMatrix<number>::eigenvalue (const unsigned int i) const
00786 {
00787 Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState());
00788 Assert (wr.size() == this->n_rows(), ExcInternalError());
00789 Assert (wi.size() == this->n_rows(), ExcInternalError());
00790 Assert (i<this->n_rows(), ExcIndexRange(i,0,this->n_rows()));
00791
00792 return std::complex<number>(wr[i], wi[i]);
00793 }
00794
00795
00796 template <typename number>
00797 inline number
00798 LAPACKFullMatrix<number>::singular_value (const unsigned int i) const
00799 {
00800 Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state));
00801 AssertIndexRange(i,wr.size());
00802
00803 return wr[i];
00804 }
00805
00806
00807
00808 DEAL_II_NAMESPACE_CLOSE
00809
00810 #endif
00811