lapack_full_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2019 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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/base/thread_management.h>

#include <deal.II/lac/lapack_support.h>
#include <deal.II/lac/vector_memory.h>

#include <complex>
#include <memory>
#include <vector>

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:
  using size_type = std::make_unsigned<types::blas_int>::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);

  void
  apply_givens_rotation(const std::array<number, 3> &csr,
                        const size_type              i,
                        const size_type              k,
                        const bool                   left = true);

  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
  transpose(LAPACKFullMatrix<number> &B) 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<typename numbers::NumberTraits<number>::real_type> 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());

  if (numbers::NumberTraits<number>::is_complex)
    return std::complex<number>(wi[i]);
  else
    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