lapack_full_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2018 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/thread_management.h>

#include <memory>
#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 std::make_unsigned<types::blas_int>::type 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 number d);

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

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

  void set (const size_type i,
            const size_type j,
            const number value);

  void add (const number                    a,
            const LAPACKFullMatrix<number> &B);

  void rank1_update(const number a,
                    const Vector<number> &v);

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

  void reinit (const size_type size);

  void grow_or_shrink (const size_type size);

  void remove_row_and_column (const size_type row, const size_type col);

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

  void set_property(const LAPACKSupport::Property property);

  size_type m () const;

  size_type 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 Tmmult (LAPACKFullMatrix<number>       &C,
               const LAPACKFullMatrix<number> &B,
               const Vector<number>           &V,
               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 scale_rows(const Vector<number> &V);

  void compute_lu_factorization ();

  void compute_cholesky_factorization ();

  number reciprocal_condition_number(const number l1_norm) const;

  number reciprocal_condition_number() const;

  number determinant () const;

  number l1_norm() const;

  number linfty_norm() const;

  number frobenius_norm() const;

  number trace() const;

  void invert ();

  void solve(Vector<number> &v,
             const bool transposed = false) const;

  void solve(LAPACKFullMatrix<number> &B,
             const bool transposed = false) const;

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

  DEAL_II_DEPRECATED
  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 types::blas_int             itype = 1);

  void compute_generalized_eigenvalues_symmetric (LAPACKFullMatrix<number>     &B,
                                                  std::vector<Vector<number> > &eigenvectors,
                                                  const types::blas_int             itype = 1);

  void compute_svd ();

  void compute_inverse_svd (const double threshold = 0.);

  void compute_inverse_svd_with_kernel (const unsigned int kernel_size);

  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:
  number norm(const char type) const;

  LAPACKSupport::State state;

  LAPACKSupport::Property property;

  mutable std::vector<number> work;

  mutable std::vector<types::blas_int> iwork;

  std::vector<types::blas_int> ipiv;

  std::vector<number> inv_work;

  std::vector<number> wr;

  std::vector<number> wi;

  std::vector<number> vl;

  std::vector<number> vr;

  std::unique_ptr<LAPACKFullMatrix<number> > svd_u;

  std::unique_ptr<LAPACKFullMatrix<number> > svd_vt;

  mutable Threads::Mutex mutex;
};



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
void LAPACKFullMatrix<number>::set (const size_type i,
                                    const size_type j,
                                    const number value)
{
  (*this)(i,j) = value;
}


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

template <typename number>
inline
typename LAPACKFullMatrix<number>::size_type
LAPACKFullMatrix<number>::n () const
{
  return static_cast<size_type>(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());
  AssertIndexRange (i,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