scalapack.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 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_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/process_grid.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/lapack_support.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/thread_management.h>
#include <mpi.h>

#include <memory>

DEAL_II_NAMESPACE_OPEN

template <typename NumberType>
class ScaLAPACKMatrix : protected TransposeTable<NumberType>
{
public:

  typedef unsigned int size_type;

  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() = 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_to (FullMatrix<NumberType> &matrix) 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 char *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 char *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 char *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 char *filename);

  /*
   * Stores the distributed matrix in @p filename
   * using parallel routines
   */
  void save_parallel(const char *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 char *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