include/deal.II/lac/lapack_full_matrix.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: lapack_full_matrix.h 25374 2012-04-03 00:33:32Z kanschat @f$
00003 //
00004 //    Copyright (C) 2005, 2006, 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__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 // forward declarations
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 
 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