scalapack.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2023 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/mpi_stub.h>
#  include <deal.II/base/mutex.h>
#  include <deal.II/base/process_grid.h>
 
#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/lapack_full_matrix.h>
#  include <deal.II/lac/lapack_support.h>
 
#  include <limits>
#  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