HODLRMatrix.hpp

Go to the documentation of this file.

/*
 * STRUMPACK -- STRUctured Matrices PACKage, Copyright (c) 2014, The
 * Regents of the University of California, through Lawrence Berkeley
 * National Laboratory (subject to receipt of any required approvals
 * from the U.S. Dept. of Energy).  All rights reserved.
 *
 * If you have questions about your rights to use or distribute this
 * software, please contact Berkeley Lab's Technology Transfer
 * Department at TTD@lbl.gov.
 *
 * NOTICE. This software is owned by the U.S. Department of Energy. As
 * such, the U.S. Government has been granted for itself and others
 * acting on its behalf a paid-up, nonexclusive, irrevocable,
 * worldwide license in the Software to reproduce, prepare derivative
 * works, and perform publicly and display publicly.  Beginning five
 * (5) years after the date permission to assert copyright is obtained
 * from the U.S. Department of Energy, and subject to any subsequent
 * five (5) year renewals, the U.S. Government is granted for itself
 * and others acting on its behalf a paid-up, nonexclusive,
 * irrevocable, worldwide license in the Software to reproduce,
 * prepare derivative works, distribute copies to the public, perform
 * publicly and display publicly, and to permit others to do so.
 *
 * Developers: Pieter Ghysels, Francois-Henry Rouet, Xiaoye S. Li.
 *             (Lawrence Berkeley National Lab, Computational Research
 *             Division).
 *
 */
#ifndef STRUMPACK_HODLR_MATRIX_HPP
#define STRUMPACK_HODLR_MATRIX_HPP
 
#include <cassert>
#include <functional>
 
#include "HSS/HSSPartitionTree.hpp"
#include "kernel/Kernel.hpp"
#include "dense/DistributedMatrix.hpp"
#include "HODLROptions.hpp"
#include "sparse/CSRGraph.hpp"
 
namespace strumpack {
 
  namespace HODLR {
 
    struct ExtractionMeta {
      std::unique_ptr<int[]> iwork;
      int Ninter, Nallrows, Nallcols, Nalldat_loc,
        *allrows, *allcols, *rowids, *colids, *pgids, Npmap, *pmaps;
    };
 
    template<typename scalar_t> class HODLRMatrix {
      using DenseM_t = DenseMatrix<scalar_t>;
      using DenseMW_t = DenseMatrixWrapper<scalar_t>;
      using DistM_t = DistributedMatrix<scalar_t>;
      using DistMW_t = DistributedMatrixWrapper<scalar_t>;
      using opts_t = HODLROptions<scalar_t>;
      using Vec_t = std::vector<std::size_t>;
      using VecVec_t = std::vector<std::vector<std::size_t>>;
      using F2Cptr = void*;
 
    public:
      using real_t = typename RealType<scalar_t>::value_type;
      using mult_t = typename std::function
        <void(Trans, const DenseM_t&, DenseM_t&)>;
      using elem_t = typename std::function
        <scalar_t(std::size_t i, std::size_t j)>;
      using delem_blocks_t = typename std::function
        <void(VecVec_t& I, VecVec_t& J, std::vector<DistMW_t>& B,
              ExtractionMeta&)>;
      using elem_blocks_t = typename std::function
        <void(VecVec_t& I, VecVec_t& J, std::vector<DenseMW_t>& B,
              ExtractionMeta&)>;
 
      HODLRMatrix() {}
 
      HODLRMatrix(const MPIComm& c, kernel::Kernel<real_t>& K,
                  const opts_t& opts);
 
      HODLRMatrix(const MPIComm& c, const HSS::HSSPartitionTree& tree,
                  const std::function<scalar_t(int i, int j)>& Aelem,
                  const opts_t& opts);
 
 
      HODLRMatrix(const MPIComm& c, const HSS::HSSPartitionTree& tree,
                  const elem_blocks_t& Aelem, const opts_t& opts);
 
 
      HODLRMatrix(const MPIComm& c, const HSS::HSSPartitionTree& tree,
                  const std::function<
                  void(Trans op,const DenseM_t& R,DenseM_t& S)>& Amult,
                  const opts_t& opts);
 
      HODLRMatrix(const MPIComm& c, const HSS::HSSPartitionTree& tree,
                  const opts_t& opts);
 
      template<typename integer_t>
      HODLRMatrix(const MPIComm& c, const HSS::HSSPartitionTree& tree,
                  const CSRGraph<integer_t>& graph, const opts_t& opts);
 
      HODLRMatrix(const HODLRMatrix<scalar_t>& h) = delete;
 
      HODLRMatrix(HODLRMatrix<scalar_t>&& h) { *this = std::move(h); }
 
      virtual ~HODLRMatrix();
 
      HODLRMatrix<scalar_t>& operator=(const HODLRMatrix<scalar_t>& h) = delete;
 
      HODLRMatrix<scalar_t>& operator=(HODLRMatrix<scalar_t>&& h);
 
      std::size_t rows() const { return rows_; }
      std::size_t cols() const { return cols_; }
      std::size_t lrows() const { return lrows_; }
      std::size_t begin_row() const { return dist_[c_->rank()]; }
      std::size_t end_row() const { return dist_[c_->rank()+1]; }
      const MPIComm& Comm() const { return *c_; }
 
      std::size_t memory() const {
        return get_stat("Mem_Fill") * 1024 * 1024;
      }
 
      std::size_t total_memory() const {
        return c_->all_reduce(memory(), MPI_SUM);
      }
 
      std::size_t factor_memory() const {
        return get_stat("Mem_Factor") * 1.e6;
      }
 
      std::size_t total_factor_memory() const {
        return c_->all_reduce(factor_memory(), MPI_SUM);
      }
 
      std::size_t rank() const { return get_stat("Rank_max"); }
 
      std::size_t max_rank() const {
        return c_->all_reduce(rank(), MPI_MAX);
      }
 
 
      double get_stat(const std::string& name) const;
 
      void print_stats();
 
      void compress
      (const std::function<void(Trans op,const DenseM_t& R,DenseM_t& S)>& Amult);
 
      void compress
      (const std::function<void(Trans op,const DenseM_t& R,DenseM_t& S)>& Amult,
       int rank_guess);
 
      void compress(const delem_blocks_t& Aelem);
 
      void compress(const elem_blocks_t& Aelem);
 
      void mult(Trans op, const DenseM_t& X, DenseM_t& Y) const;
 
      void mult(Trans op, const DistM_t& X, DistM_t& Y) const;
 
      void factor();
 
      long long int solve(const DenseM_t& B, DenseM_t& X) const;
 
      long long int solve(const DistM_t& B, DistM_t& X) const;
 
      long long int inv_mult(Trans op, const DenseM_t& B, DenseM_t& X) const;
 
      void extract_elements(const VecVec_t& I, const VecVec_t& J,
                            std::vector<DistM_t>& B);
      void extract_elements(const VecVec_t& I, const VecVec_t& J,
                            std::vector<DenseM_t>& B);
 
      void extract_elements(const Vec_t& I, const Vec_t& J, DenseM_t& B);
 
      DistM_t dense(const BLACSGrid* g) const;
 
      void set_sampling_parameter(double sample_param);
 
      DenseM_t redistribute_2D_to_1D(const DistM_t& R) const;
      void redistribute_2D_to_1D(const DistM_t& R2D, DenseM_t& R1D) const;
      void redistribute_1D_to_2D(const DenseM_t& S1D, DistM_t& S2D) const;
 
      DenseMatrix<scalar_t> gather_from_1D(const DenseM_t& A) const;
      DenseMatrix<scalar_t> all_gather_from_1D(const DenseM_t& A) const;
      DenseMatrix<scalar_t> scatter_to_1D(const DenseM_t& A) const;
 
      const std::vector<int>& perm() const { return perm_; }
      const std::vector<int>& iperm() const { return iperm_; }
 
    private:
      F2Cptr ho_bf_ = nullptr;     // HODLR handle returned by Fortran code
      F2Cptr options_ = nullptr;   // options structure returned by Fortran code
      F2Cptr stats_ = nullptr;     // statistics structure returned by Fortran code
      F2Cptr msh_ = nullptr;       // mesh structure returned by Fortran code
      F2Cptr kerquant_ = nullptr;  // kernel quantities structure returned by Fortran code
      F2Cptr ptree_ = nullptr;     // process tree returned by Fortran code
      MPI_Fint Fcomm_;             // the fortran MPI communicator
      const MPIComm* c_;
      int rows_ = 0, cols_ = 0, lrows_ = 0, lvls_ = 0;
      std::vector<int> perm_, iperm_; // permutation used by the HODLR code
      std::vector<int> dist_;         // begin rows of each rank
      std::vector<int> leafs_;        // leaf sizes of the tree
 
      void options_init(const opts_t& opts);
      void perm_init();
      void dist_init();
 
      template<typename S> friend class ButterflyMatrix;
    };
 
 
    template<typename scalar_t> struct AelemCommPtrs {
      const typename HODLRMatrix<scalar_t>::delem_blocks_t* Aelem;
      const MPIComm* c;
      const BLACSGrid* gl;
      const BLACSGrid* g0;
    };
 
    template<typename scalar_t> void HODLR_block_evaluation
    (int* Ninter, int* Nallrows, int* Nallcols, int* Nalldat_loc,
     int* allrows, int* allcols, scalar_t* alldat_loc,
     int* rowids, int* colids, int* pgids, int* Npmap, int* pmaps,
     void* AC);
 
    template<typename scalar_t> void HODLR_block_evaluation_seq
    (int* Ninter, int* Nallrows, int* Nallcols, int* Nalldat_loc,
     int* allrows, int* allcols, scalar_t* alldat_loc,
     int* rowids, int* colids, int* pgids, int* Npmap, int* pmaps,
     void* f);
 
  } // end namespace HODLR
} // end namespace strumpack
 
#endif // STRUMPACK_HODLR_MATRIX_HPP