scalapack.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 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_scalapack_h
#define dealii_scalapack_h

#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_SCALAPACK

#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/mpi.h>
#  include <deal.II/base/process_grid.h>
#  include <deal.II/base/thread_management.h>

#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/lapack_full_matrix.h>
#  include <deal.II/lac/lapack_support.h>

#  include <mpi.h>

#  include <memory>

DEAL_II_NAMESPACE_OPEN

template <typename NumberType>
class ScaLAPACKMatrix : protected TransposeTable<NumberType>
{
public:
  using size_type = unsigned int;

  ScaLAPACKMatrix(
    const size_type                                           n_rows,
    const size_type                                           n_columns,
    const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
    const size_type               row_block_size    = 32,
    const size_type               column_block_size = 32,
    const LAPACKSupport::Property property = LAPACKSupport::Property::general);

  ScaLAPACKMatrix(
    const size_type                                           size,
    const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
    const size_type                                           block_size = 32,
    const LAPACKSupport::Property                             property =
      LAPACKSupport::Property::symmetric);

  ScaLAPACKMatrix(
    const std::string &                                       filename,
    const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
    const size_type row_block_size    = 32,
    const size_type column_block_size = 32);

  ~ScaLAPACKMatrix() override = default;

  void
  reinit(
    const size_type                                           n_rows,
    const size_type                                           n_columns,
    const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
    const size_type               row_block_size    = 32,
    const size_type               column_block_size = 32,
    const LAPACKSupport::Property property = LAPACKSupport::Property::general);

  void
  reinit(const size_type                                           size,
         const std::shared_ptr<const Utilities::MPI::ProcessGrid> &process_grid,
         const size_type               block_size = 32,
         const LAPACKSupport::Property property =
           LAPACKSupport::Property::symmetric);

  void
  set_property(const LAPACKSupport::Property property);

  LAPACKSupport::Property
  get_property() const;

  LAPACKSupport::State
  get_state() const;

  ScaLAPACKMatrix<NumberType> &
  operator=(const FullMatrix<NumberType> &);

  void
  copy_from(const LAPACKFullMatrix<NumberType> &matrix,
            const unsigned int                  rank);

  void
  copy_to(FullMatrix<NumberType> &matrix) const;

  void
  copy_to(LAPACKFullMatrix<NumberType> &matrix, const unsigned int rank) const;

  void
  copy_to(ScaLAPACKMatrix<NumberType> &dest) const;

  void
  copy_to(ScaLAPACKMatrix<NumberType> &                B,
          const std::pair<unsigned int, unsigned int> &offset_A,
          const std::pair<unsigned int, unsigned int> &offset_B,
          const std::pair<unsigned int, unsigned int> &submatrix_size) const;

  void
  copy_transposed(const ScaLAPACKMatrix<NumberType> &B);

  void
  add(const ScaLAPACKMatrix<NumberType> &B,
      const NumberType                   a           = 0.,
      const NumberType                   b           = 1.,
      const bool                         transpose_B = false);

  void
  add(const NumberType b, const ScaLAPACKMatrix<NumberType> &B);

  void
  Tadd(const NumberType b, const ScaLAPACKMatrix<NumberType> &B);

  void
  mult(const NumberType                   b,
       const ScaLAPACKMatrix<NumberType> &B,
       const NumberType                   c,
       ScaLAPACKMatrix<NumberType> &      C,
       const bool                         transpose_A = false,
       const bool                         transpose_B = false) const;

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

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

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

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

  void
  save(const std::string &                          filename,
       const std::pair<unsigned int, unsigned int> &chunk_size =
         std::make_pair(numbers::invalid_unsigned_int,
                        numbers::invalid_unsigned_int)) const;

  void
  load(const std::string &filename);

  void
  compute_cholesky_factorization();

  void
  compute_lu_factorization();

  void
  invert();

  std::vector<NumberType>
  eigenpairs_symmetric_by_index(
    const std::pair<unsigned int, unsigned int> &index_limits,
    const bool                                   compute_eigenvectors);

  std::vector<NumberType>
  eigenpairs_symmetric_by_value(
    const std::pair<NumberType, NumberType> &value_limits,
    const bool                               compute_eigenvectors);

  std::vector<NumberType>
  eigenpairs_symmetric_by_index_MRRR(
    const std::pair<unsigned int, unsigned int> &index_limits,
    const bool                                   compute_eigenvectors);

  std::vector<NumberType>
  eigenpairs_symmetric_by_value_MRRR(
    const std::pair<NumberType, NumberType> &value_limits,
    const bool                               compute_eigenvectors);

  std::vector<NumberType>
  compute_SVD(ScaLAPACKMatrix<NumberType> *U  = nullptr,
              ScaLAPACKMatrix<NumberType> *VT = nullptr);

  void
  least_squares(ScaLAPACKMatrix<NumberType> &B, const bool transpose = false);

  unsigned int
  pseudoinverse(const NumberType ratio);

  NumberType
  reciprocal_condition_number(const NumberType a_norm) const;

  NumberType
  l1_norm() const;

  NumberType
  linfty_norm() const;

  NumberType
  frobenius_norm() const;

  size_type
  m() const;

  size_type
  n() const;

  unsigned int
  local_m() const;

  unsigned int
  local_n() const;

  unsigned int
  global_row(const unsigned int loc_row) const;

  unsigned int
  global_column(const unsigned int loc_column) const;

  NumberType
  local_el(const unsigned int loc_row, const unsigned int loc_column) const;

  NumberType &
  local_el(const unsigned int loc_row, const unsigned int loc_column);

  template <class InputVector>
  void
  scale_columns(const InputVector &factors);

  template <class InputVector>
  void
  scale_rows(const InputVector &factors);

private:
  NumberType
  norm_symmetric(const char type) const;

  NumberType
  norm_general(const char type) const;

  std::vector<NumberType>
  eigenpairs_symmetric(
    const bool                                   compute_eigenvectors,
    const std::pair<unsigned int, unsigned int> &index_limits =
      std::make_pair(numbers::invalid_unsigned_int,
                     numbers::invalid_unsigned_int),
    const std::pair<NumberType, NumberType> &value_limits =
      std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),
                     std::numeric_limits<NumberType>::quiet_NaN()));

  std::vector<NumberType>
  eigenpairs_symmetric_MRRR(
    const bool                                   compute_eigenvectors,
    const std::pair<unsigned int, unsigned int> &index_limits =
      std::make_pair(numbers::invalid_unsigned_int,
                     numbers::invalid_unsigned_int),
    const std::pair<NumberType, NumberType> &value_limits =
      std::make_pair(std::numeric_limits<NumberType>::quiet_NaN(),
                     std::numeric_limits<NumberType>::quiet_NaN()));

  /*
   * Stores the distributed matrix in @p filename
   * using serial routines
   */
  void
  save_serial(const std::string &                          filename,
              const std::pair<unsigned int, unsigned int> &chunk_size) const;

  /*
   * Loads the distributed matrix from file @p filename
   * using serial routines
   */
  void
  load_serial(const std::string &filename);

  /*
   * Stores the distributed matrix in @p filename
   * using parallel routines
   */
  void
  save_parallel(const std::string &                          filename,
                const std::pair<unsigned int, unsigned int> &chunk_size) const;

  /*
   * Loads the distributed matrix from file @p filename
   * using parallel routines
   */
  void
  load_parallel(const std::string &filename);

  LAPACKSupport::State state;

  LAPACKSupport::Property property;

  std::shared_ptr<const Utilities::MPI::ProcessGrid> grid;

  int n_rows;

  int n_columns;

  int row_block_size;

  int column_block_size;

  int n_local_rows;

  int n_local_columns;

  int descriptor[9];

  mutable std::vector<NumberType> work;

  mutable std::vector<int> iwork;

  std::vector<int> ipiv;

  const char uplo;

  const int first_process_row;

  const int first_process_column;

  const int submatrix_row;

  const int submatrix_column;

  mutable Threads::Mutex mutex;
};

// ----------------------- inline functions ----------------------------

#  ifndef DOXYGEN

template <typename NumberType>
inline NumberType
ScaLAPACKMatrix<NumberType>::local_el(const unsigned int loc_row,
                                      const unsigned int loc_column) const
{
  return (*this)(loc_row, loc_column);
}



template <typename NumberType>
inline NumberType &
ScaLAPACKMatrix<NumberType>::local_el(const unsigned int loc_row,
                                      const unsigned int loc_column)
{
  return (*this)(loc_row, loc_column);
}


template <typename NumberType>
inline unsigned int
ScaLAPACKMatrix<NumberType>::m() const
{
  return n_rows;
}



template <typename NumberType>
inline unsigned int
ScaLAPACKMatrix<NumberType>::n() const
{
  return n_columns;
}



template <typename NumberType>
unsigned int
ScaLAPACKMatrix<NumberType>::local_m() const
{
  return n_local_rows;
}



template <typename NumberType>
unsigned int
ScaLAPACKMatrix<NumberType>::local_n() const
{
  return n_local_columns;
}


#  endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_SCALAPACK

#endif