ScalapackWrapper.h

Go to the documentation of this file.

/******************************************************************************
 * Copyright (c) 2021.                                                        *
 * The Regents of the University of Michigan and DFT-EFE developers.          *
 *                                                                            *
 * This file is part of the DFT-EFE code.                                     *
 *                                                                            *
 * DFT-EFE is free software: you can redistribute it and/or modify            *
 *   it under the terms of the Lesser GNU General Public License as           *
 *   published by the Free Software Foundation, either version 3 of           *
 *   the License, or (at your option) any later version.                      *
 *                                                                            *
 * DFT-EFE is distributed in the hope that it will be useful, but             *
 *   WITHOUT ANY WARRANTY; without even the implied warranty                  *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                     *
 *   See the Lesser GNU General Public License for more details.              *
 *                                                                            *
 * You should have received a copy of the GNU Lesser General Public           *
 *   License at the top level of DFT-EFE distribution.  If not, see           *
 *   <https://www.gnu.org/licenses/>.                                         *
 ******************************************************************************/
 
 
 
#ifndef ScaLAPACKMatrix_h
#define ScaLAPACKMatrix_h
 
#include "ProcessGrid.h"
#include "LapackSupport.h"
#include <memory>
#include "LinearAlgebraTypes.h"
 
namespace dftefe
{
  namespace linearAlgebra
  {
    template <typename NumberType>
    class ScaLAPACKMatrix
    {
    public:
      using size_type_scalapack = unsigned int;
      ScaLAPACKMatrix(const size_type_scalapack                 n_rows,
                      const size_type_scalapack                 n_columns,
                      const std::shared_ptr<const ProcessGrid> &process_grid,
                      const size_type_scalapack     row_block_size    = 32,
                      const size_type_scalapack     column_block_size = 32,
                      const LAPACKSupport::Property property =
                        LAPACKSupport::Property::general);
 
      ScaLAPACKMatrix(const size_type_scalapack                 size,
                      const std::shared_ptr<const ProcessGrid> &process_grid,
                      const size_type_scalapack                 block_size = 32,
                      const LAPACKSupport::Property             property =
                        LAPACKSupport::Property::hermitian);
 
      ~ScaLAPACKMatrix() = default;
 
      void
      reinit(const size_type_scalapack                 n_rows,
             const size_type_scalapack                 n_columns,
             const std::shared_ptr<const ProcessGrid> &process_grid,
             const size_type_scalapack                 row_block_size    = 32,
             const size_type_scalapack                 column_block_size = 32,
             const LAPACKSupport::Property             property =
               LAPACKSupport::Property::general);
 
      void
      reinit(const size_type_scalapack                 size,
             const std::shared_ptr<const ProcessGrid> &process_grid,
             const size_type_scalapack                 block_size = 32,
             const LAPACKSupport::Property             property =
               LAPACKSupport::Property::hermitian);
 
 
      void
      set_property(const LAPACKSupport::Property property);
 
      LAPACKSupport::Property
      get_property() const;
 
      LAPACKSupport::State
      get_state() const;
 
 
      void
      copy_to(ScaLAPACKMatrix<NumberType> &dest) const;
 
 
      void
      conjugate();
 
 
      void
      add(const ScaLAPACKMatrix<NumberType> &B,
          const NumberType                   a           = 0.,
          const NumberType                   b           = 1.,
          const bool                         transpose_B = false);
 
 
      void
      zadd(const ScaLAPACKMatrix<NumberType> &B,
           const NumberType                   a                     = 0.,
           const NumberType                   b                     = 1.,
           const bool                         conjugate_transpose_B = false);
 
      void
      copy_transposed(const ScaLAPACKMatrix<NumberType> &B);
 
 
      void
      copy_conjugate_transposed(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
      zmult(const NumberType                   b,
            const ScaLAPACKMatrix<NumberType> &B,
            const NumberType                   c,
            ScaLAPACKMatrix<NumberType> &      C,
            const bool                         conjugate_transpose_A = false,
            const bool conjugate_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
      zmmult(ScaLAPACKMatrix<NumberType> &      C,
             const ScaLAPACKMatrix<NumberType> &B,
             const bool                         adding = false) const;
 
      void
      zCmmult(ScaLAPACKMatrix<NumberType> &      C,
              const ScaLAPACKMatrix<NumberType> &B,
              const bool                         adding = false) const;
 
      void
      zmCmult(ScaLAPACKMatrix<NumberType> &      C,
              const ScaLAPACKMatrix<NumberType> &B,
              const bool                         adding = false) const;
 
      void
      zCmCmult(ScaLAPACKMatrix<NumberType> &      C,
               const ScaLAPACKMatrix<NumberType> &B,
               const bool                         adding = false) const;
 
 
      size_type_scalapack
      m() const;
 
      size_type_scalapack
      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);
 
      ScalapackError
      compute_cholesky_factorization();
 
      ScalapackError
      compute_lu_factorization();
 
      ScalapackError
      invert();
 
      void
      scale_columns(const std::vector<NumberType> &factors);
 
      void
      scale_rows(const std::vector<NumberType> &factors);
 
 
      void
      scale_columns_realfactors(const std::vector<double> &factors);
 
      void
      scale_rows_realfactors(const std::vector<double> &factors);
 
      std::vector<double>
      eigenpairs_hermitian_by_index(
        const std::pair<unsigned int, unsigned int> &index_limits,
        const bool                                   compute_eigenvectors,
        ScalapackError &                             scalapackError);
 
      std::vector<double>
      eigenpairs_hermitian_by_index_MRRR(
        const std::pair<unsigned int, unsigned int> &index_limits,
        const bool                                   compute_eigenvectors,
        ScalapackError &                             scalapackError);
 
    private:
      std::vector<double>
      eigenpairs_hermitian(
        const bool                                   compute_eigenvectors,
        ScalapackError &                             scalapackError,
        const std::pair<unsigned int, unsigned int> &index_limits =
          std::make_pair(std::numeric_limits<unsigned int>::max(),
                         std::numeric_limits<unsigned int>::max()),
        const std::pair<double, double> &value_limits =
          std::make_pair(std::numeric_limits<double>::quiet_NaN(),
                         std::numeric_limits<double>::quiet_NaN()));
 
      std::vector<double>
      eigenpairs_hermitian_MRRR(
        const bool                                   compute_eigenvectors,
        ScalapackError &                             scalapackError,
        const std::pair<unsigned int, unsigned int> &index_limits =
          std::make_pair(std::numeric_limits<unsigned int>::max(),
                         std::numeric_limits<unsigned int>::max()),
        const std::pair<double, double> &value_limits =
          std::make_pair(std::numeric_limits<double>::quiet_NaN(),
                         std::numeric_limits<double>::quiet_NaN()));
 
 
      std::vector<NumberType> values;
 
      LAPACKSupport::State state;
 
      LAPACKSupport::Property property;
 
      std::shared_ptr<const 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;
    };
 
    // ----------------------- inline functions ----------------------------
 
    template <typename NumberType>
    inline NumberType
    ScaLAPACKMatrix<NumberType>::local_el(const unsigned int loc_row,
                                          const unsigned int loc_column) const
    {
      return values[loc_column * n_local_rows + loc_row];
      // 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 values[loc_column * n_local_rows + loc_row];
      // 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;
    }
 
  } // namespace linearAlgebra
} // namespace dftefe
#include <linearAlgebra/ScalapackWrapper.t.cpp>
#endif // ScaLAPACKMatrix_h