lapack_full_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__lapack_full_matrix_h
#define dealii__lapack_full_matrix_h


#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/table.h>
#include <deal.II/lac/lapack_support.h>
#include <deal.II/lac/vector_memory.h>

#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <vector>
#include <complex>

DEAL_II_NAMESPACE_OPEN

// forward declarations
template<typename number> class Vector;
template<typename number> class BlockVector;
template<typename number> class FullMatrix;
template<typename number> class SparseMatrix;


template <typename number>
class LAPACKFullMatrix : public TransposeTable<number>
{
public:

  typedef types::global_dof_index size_type;

  explicit LAPACKFullMatrix (const size_type size = 0);


  LAPACKFullMatrix (const size_type rows,
                    const size_type cols);


  LAPACKFullMatrix (const LAPACKFullMatrix &);

  LAPACKFullMatrix<number> &
  operator = (const LAPACKFullMatrix<number> &);

  template <typename number2>
  LAPACKFullMatrix<number> &
  operator = (const FullMatrix<number2> &);

  template <typename number2>
  LAPACKFullMatrix<number> &
  operator = (const SparseMatrix<number2> &);

  LAPACKFullMatrix<number> &
  operator = (const double d);

  LAPACKFullMatrix<number> &
  operator*= (const number factor);

  LAPACKFullMatrix<number> &
  operator/= (const number factor);

  template <typename MatrixType>
  void copy_from (const MatrixType &);

  void reinit (const size_type size);

  void reinit (const size_type rows,
               const size_type cols);

  unsigned int m () const;

  unsigned int n () const;

  template<typename MatrixType>
  void fill (const MatrixType &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);


  template <typename number2>
  void vmult (Vector<number2>       &w,
              const Vector<number2> &v,
              const bool             adding = false) const;

  void vmult (Vector<number>       &w,
              const Vector<number> &v,
              const bool            adding = false) const;

  template <typename number2>
  void vmult_add (Vector<number2>       &w,
                  const Vector<number2> &v) const;

  void vmult_add (Vector<number>       &w,
                  const Vector<number> &v) const;

  template <typename number2>
  void Tvmult (Vector<number2>       &w,
               const Vector<number2> &v,
               const bool             adding=false) const;

  void Tvmult (Vector<number>       &w,
               const Vector<number> &v,
               const bool            adding=false) const;

  template <typename number2>
  void Tvmult_add (Vector<number2>       &w,
                   const Vector<number2> &v) const;

  void Tvmult_add (Vector<number>       &w,
                   const Vector<number> &v) const;


  void mmult (LAPACKFullMatrix<number>       &C,
              const LAPACKFullMatrix<number> &B,
              const bool                      adding=false) const;

  void mmult (FullMatrix<number>             &C,
              const LAPACKFullMatrix<number> &B,
              const bool                      adding=false) const;

  void Tmmult (LAPACKFullMatrix<number>       &C,
               const LAPACKFullMatrix<number> &B,
               const bool                      adding=false) const;

  void Tmmult (FullMatrix<number>             &C,
               const LAPACKFullMatrix<number> &B,
               const bool                      adding=false) const;

  void mTmult (LAPACKFullMatrix<number>       &C,
               const LAPACKFullMatrix<number> &B,
               const bool                      adding=false) const;

  void mTmult (FullMatrix<number>             &C,
               const LAPACKFullMatrix<number> &B,
               const bool                      adding=false) const;

  void TmTmult (LAPACKFullMatrix<number>       &C,
                const LAPACKFullMatrix<number> &B,
                const bool                      adding=false) const;

  void TmTmult (FullMatrix<number>             &C,
                const LAPACKFullMatrix<number> &B,
                const bool                      adding=false) const;

  void compute_lu_factorization ();

  void invert ();

  void apply_lu_factorization (Vector<number> &v,
                               const bool      transposed) const;

  void apply_lu_factorization (LAPACKFullMatrix<number> &B,
                               const bool                transposed) const;

  void compute_eigenvalues (const bool right_eigenvectors = false,
                            const bool left_eigenvectors  = false);

  void compute_eigenvalues_symmetric (const number        lower_bound,
                                      const number        upper_bound,
                                      const number        abs_accuracy,
                                      Vector<number>     &eigenvalues,
                                      FullMatrix<number> &eigenvectors);

  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 compute_generalized_eigenvalues_symmetric (LAPACKFullMatrix<number>     &B,
                                                  std::vector<Vector<number> > &eigenvectors,
                                                  const int                     itype = 1);

  void compute_svd ();

  void compute_inverse_svd (const double threshold = 0.);

  std::complex<number>
  eigenvalue (const size_type i) const;

  number
  singular_value (const size_type i) const;

  void print_formatted (std::ostream       &out,
                        const unsigned int  precision   = 3,
                        const bool          scientific  = true,
                        const unsigned int  width       = 0,
                        const char         *zero_string = " ",
                        const double        denominator = 1.,
                        const double        threshold   = 0.) const;

private:

  LAPACKSupport::State state;

  LAPACKSupport::Properties properties;

  mutable std::vector<number> work;

  std::vector<int> ipiv;

  std::vector<number> inv_work;

  std::vector<number> wr;

  std::vector<number> wi;

  std::vector<number> vl;

  std::vector<number> vr;

  std_cxx11::shared_ptr<LAPACKFullMatrix<number> > svd_u;

  std_cxx11::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
};



template <typename number>
class PreconditionLU
  :
  public Subscriptor
{
public:
  void initialize(const LAPACKFullMatrix<number> &);
  void initialize(const LAPACKFullMatrix<number> &,
                  VectorMemory<Vector<number> > &);
  void vmult(Vector<number> &, const Vector<number> &) const;
  void Tvmult(Vector<number> &, const Vector<number> &) const;
  void vmult(BlockVector<number> &,
             const BlockVector<number> &) const;
  void Tvmult(BlockVector<number> &,
              const BlockVector<number> &) const;
private:
  SmartPointer<const LAPACKFullMatrix<number>,PreconditionLU<number> > matrix;
  SmartPointer<VectorMemory<Vector<number> >,PreconditionLU<number> > mem;
};

/*---------------------- Inline functions -----------------------------------*/

template <typename number>
inline
unsigned int
LAPACKFullMatrix<number>::m () const
{
  return this->n_rows ();
}

template <typename number>
inline
unsigned int
LAPACKFullMatrix<number>::n () const
{
  return this->n_cols ();
}

template <typename number>
template <typename MatrixType>
inline void
LAPACKFullMatrix<number>::copy_from (const MatrixType &M)
{
  this->reinit (M.m(), M.n());

  // loop over the elements of the argument matrix row by row, as suggested
  // in the documentation of the sparse matrix iterator class, and
  // copy them into the current object
  for (size_type row = 0; row < M.m(); ++row)
    {
      const typename MatrixType::const_iterator end_row = M.end(row);
      for (typename MatrixType::const_iterator entry = M.begin(row);
           entry != end_row; ++entry)
        this->el(row, entry->column()) = entry->value();
    }

  state = LAPACKSupport::matrix;
}



template <typename number>
template <typename MatrixType>
inline void
LAPACKFullMatrix<number>::fill (const MatrixType &M,
                                const size_type   dst_offset_i,
                                const size_type   dst_offset_j,
                                const size_type   src_offset_i,
                                const size_type   src_offset_j,
                                const number      factor,
                                const bool        transpose)
{
  // loop over the elements of the argument matrix row by row, as suggested
  // in the documentation of the sparse matrix iterator class
  for (size_type row = src_offset_i; row < M.m(); ++row)
    {
      const typename MatrixType::const_iterator end_row = M.end(row);
      for (typename MatrixType::const_iterator entry = M.begin(row);
           entry != end_row; ++entry)
        {
          const size_type i = transpose ? entry->column() : row;
          const size_type j = transpose ? row : entry->column();

          const size_type dst_i=dst_offset_i+i-src_offset_i;
          const size_type dst_j=dst_offset_j+j-src_offset_j;
          if (dst_i<this->n_rows() && dst_j<this->n_cols())
            (*this)(dst_i, dst_j) = factor * entry->value();
        }
    }

  state = LAPACKSupport::matrix;
}


template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::vmult (Vector<number2> &,
                                 const Vector<number2> &,
                                 const bool) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::vmult must be called with a "
                    "matching Vector<double> vector type."));
}


template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::vmult_add (Vector<number2> &,
                                     const Vector<number2> &) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::vmult_add must be called with a "
                    "matching Vector<double> vector type."));
}


template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::Tvmult (Vector<number2> &,
                                  const Vector<number2> &,
                                  const bool) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::Tvmult must be called with a "
                    "matching Vector<double> vector type."));
}


template <typename number>
template <typename number2>
void
LAPACKFullMatrix<number>::Tvmult_add (Vector<number2> &,
                                      const Vector<number2> &) const
{
  Assert(false,
         ExcMessage("LAPACKFullMatrix<number>::Tvmult_add must be called with a "
                    "matching Vector<double> vector type."));
}


template <typename number>
inline std::complex<number>
LAPACKFullMatrix<number>::eigenvalue (const size_type i) const
{
  Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState());
  Assert (wr.size() == this->n_rows(), ExcInternalError());
  Assert (wi.size() == this->n_rows(), ExcInternalError());
  Assert (i<this->n_rows(), ExcIndexRange(i,0,this->n_rows()));

  return std::complex<number>(wr[i], wi[i]);
}


template <typename number>
inline number
LAPACKFullMatrix<number>::singular_value (const size_type i) const
{
  Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state));
  AssertIndexRange(i,wr.size());

  return wr[i];
}



DEAL_II_NAMESPACE_CLOSE

#endif