HSSMatrix.hpp
Go to the documentation of this file.
1 /*
2  * STRUMPACK -- STRUctured Matrices PACKage, Copyright (c) 2014, The
3  * Regents of the University of California, through Lawrence Berkeley
4  * National Laboratory (subject to receipt of any required approvals
5  * from the U.S. Dept. of Energy). All rights reserved.
6  *
7  * If you have questions about your rights to use or distribute this
8  * software, please contact Berkeley Lab's Technology Transfer
9  * Department at TTD@lbl.gov.
10  *
11  * NOTICE. This software is owned by the U.S. Department of Energy. As
12  * such, the U.S. Government has been granted for itself and others
13  * acting on its behalf a paid-up, nonexclusive, irrevocable,
14  * worldwide license in the Software to reproduce, prepare derivative
15  * works, and perform publicly and display publicly. Beginning five
16  * (5) years after the date permission to assert copyright is obtained
17  * from the U.S. Department of Energy, and subject to any subsequent
18  * five (5) year renewals, the U.S. Government is granted for itself
19  * and others acting on its behalf a paid-up, nonexclusive,
20  * irrevocable, worldwide license in the Software to reproduce,
21  * prepare derivative works, distribute copies to the public, perform
22  * publicly and display publicly, and to permit others to do so.
23  *
24  * Developers: Pieter Ghysels, Francois-Henry Rouet, Xiaoye S. Li.
25  * (Lawrence Berkeley National Lab, Computational Research
26  * Division).
27  *
28  */
37 #ifndef HSS_MATRIX_HPP
38 #define HSS_MATRIX_HPP
39 
40 #include <cassert>
41 #include <functional>
42 #include <string>
43 
44 #include "HSSBasisID.hpp"
45 #include "HSSOptions.hpp"
46 #include "HSSExtra.hpp"
47 #include "HSSMatrixBase.hpp"
48 #include "kernel/Kernel.hpp"
49 #include "HSSMatrix.sketch.hpp"
50 
51 namespace strumpack {
52  namespace HSS {
53 
54 #ifndef DOXYGEN_SHOULD_SKIP_THIS
55  // forward declaration
56  template<typename scalar_t> class HSSMatrixMPI;
57 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
58 
59 
79  template<typename scalar_t> class HSSMatrix
80  : public HSSMatrixBase<scalar_t> {
81  using real_t = typename RealType<scalar_t>::value_type;
82  using DenseM_t = DenseMatrix<scalar_t>;
83  using DenseMW_t = DenseMatrixWrapper<scalar_t>;
84  using elem_t = typename std::function
85  <void(const std::vector<std::size_t>& I,
86  const std::vector<std::size_t>& J, DenseM_t& B)>;
87  using mult_t = typename std::function
88  <void(DenseM_t& Rr, DenseM_t& Rc, DenseM_t& Sr, DenseM_t& Sc)>;
89  using opts_t = HSSOptions<scalar_t>;
90 
91  public:
95  HSSMatrix();
96 
114  HSSMatrix(const DenseM_t& A, const opts_t& opts);
115 
131  HSSMatrix(std::size_t m, std::size_t n, const opts_t& opts);
132 
143  HSSMatrix(const structured::ClusterTree& t, const opts_t& opts);
144 
154  HSSMatrix(kernel::Kernel<real_t>& K, const opts_t& opts);
155 
161  HSSMatrix(const HSSMatrix<scalar_t>& other);
162 
168  HSSMatrix<scalar_t>& operator=(const HSSMatrix<scalar_t>& other);
169 
174  HSSMatrix(HSSMatrix<scalar_t>&& other) = default;
175 
180  HSSMatrix<scalar_t>& operator=(HSSMatrix<scalar_t>&& other) = default;
181 
186  std::unique_ptr<HSSMatrixBase<scalar_t>> clone() const override;
187 
194  const HSSMatrix<scalar_t>* child(int c) const {
195  return dynamic_cast<HSSMatrix<scalar_t>*>(this->ch_[c].get());
196  }
197 
204  HSSMatrix<scalar_t>* child(int c) {
205  return dynamic_cast<HSSMatrix<scalar_t>*>(this->ch_[c].get());
206  }
207 
224  void compress(const DenseM_t& A, const opts_t& opts);
225 
270  void compress(const std::function<void(DenseM_t& Rr,
271  DenseM_t& Rc,
272  DenseM_t& Sr,
273  DenseM_t& Sc)>& Amult,
274  const std::function<void(const std::vector<std::size_t>& I,
275  const std::vector<std::size_t>& J,
276  DenseM_t& B)>& Aelem,
277  const opts_t& opts);
278 
279 
306  void compress_with_coordinates(const DenseMatrix<real_t>& coords,
307  const std::function
308  <void(const std::vector<std::size_t>& I,
309  const std::vector<std::size_t>& J,
310  DenseM_t& B)>& Aelem,
311  const opts_t& opts);
312 
317  void reset() override;
318 
322  void factor() override;
323 
334  void partial_factor();
335 
347  void solve(DenseM_t& b) const override;
348 
364  void forward_solve(WorkSolve<scalar_t>& w, const DenseM_t& b,
365  bool partial) const override;
366 
380  void backward_solve(WorkSolve<scalar_t>& w, DenseM_t& x) const override;
381 
390  DenseM_t apply(const DenseM_t& b) const;
391 
404  void mult(Trans op, const DenseM_t& x, DenseM_t& y) const override;
405 
415  DenseM_t applyC(const DenseM_t& b) const;
416 
427  scalar_t get(std::size_t i, std::size_t j) const;
428 
439  DenseM_t extract(const std::vector<std::size_t>& I,
440  const std::vector<std::size_t>& J) const;
441 
455  void extract_add(const std::vector<std::size_t>& I,
456  const std::vector<std::size_t>& J,
457  DenseM_t& B) const;
458 
459 #ifndef DOXYGEN_SHOULD_SKIP_THIS
460  void Schur_update(DenseM_t& Theta,
461  DenseM_t& DUB01,
462  DenseM_t& Phi) const;
463  void Schur_product_direct(const DenseM_t& Theta,
464  const DenseM_t& DUB01,
465  const DenseM_t& Phi,
466  const DenseM_t&_ThetaVhatC_or_VhatCPhiC,
467  const DenseM_t& R,
468  DenseM_t& Sr, DenseM_t& Sc) const;
469  void Schur_product_indirect(const DenseM_t& DUB01,
470  const DenseM_t& R1,
471  const DenseM_t& R2, const DenseM_t& Sr2,
472  const DenseM_t& Sc2,
473  DenseM_t& Sr, DenseM_t& Sc) const;
474  void delete_trailing_block() override;
475 #endif // DOXYGEN_SHOULD_SKIP_THIS
476 
477  std::size_t rank() const override;
478  std::size_t memory() const override;
479  std::size_t nonzeros() const override;
480  std::size_t levels() const override;
481 
482  void print_info(std::ostream& out=std::cout,
483  std::size_t roff=0,
484  std::size_t coff=0) const override;
485 
492  DenseM_t dense() const;
493 
494  void shift(scalar_t sigma) override;
495 
496  void draw(std::ostream& of,
497  std::size_t rlo=0, std::size_t clo=0) const override;
498 
505  void write(const std::string& fname) const;
506 
513  static HSSMatrix<scalar_t> read(const std::string& fname);
514 
515  const HSSFactors<scalar_t>& ULV() { return this->ULV_; }
516 
517  protected:
518  HSSMatrix(std::size_t m, std::size_t n,
519  const opts_t& opts, bool active);
520  HSSMatrix(const structured::ClusterTree& t,
521  const opts_t& opts, bool active);
522  HSSMatrix(std::ifstream& is);
523 
524  HSSBasisID<scalar_t> U_, V_;
525  DenseM_t D_, B01_, B10_;
526 
527  void compress_original(const DenseM_t& A,
528  const opts_t& opts);
529  void compress_original(const mult_t& Amult,
530  const elem_t& Aelem,
531  const opts_t& opts);
532  void compress_stable(const DenseM_t& A,
533  const opts_t& opts);
534  void compress_stable(const mult_t& Amult,
535  const elem_t& Aelem,
536  const opts_t& opts);
537  void compress_hard_restart(const DenseM_t& A,
538  const opts_t& opts);
539  void compress_hard_restart(const mult_t& Amult,
540  const elem_t& Aelem,
541  const opts_t& opts);
542 
543  void compress_recursive_original(DenseM_t& Rr, DenseM_t& Rc,
544  DenseM_t& Sr, DenseM_t& Sc,
545  const elem_t& Aelem,
546  const opts_t& opts,
547  WorkCompress<scalar_t>& w,
548  int dd, int depth) override;
549  void compress_recursive_stable(DenseM_t& Rr, DenseM_t& Rc,
550  DenseM_t& Sr, DenseM_t& Sc,
551  const elem_t& Aelem,
552  const opts_t& opts,
553  WorkCompress<scalar_t>& w,
554  int d, int dd, int depth) override;
555  // SJLT_Matrix<scalar_t,int>* S=nullptr
556  void compute_local_samples(DenseM_t& Rr, DenseM_t& Rc,
557  DenseM_t& Sr, DenseM_t& Sc,
558  WorkCompress<scalar_t>& w,
559  int d0, int d, int depth,
560  SJLTMatrix<scalar_t, int>* S=nullptr);
561  bool compute_U_V_bases(DenseM_t& Sr, DenseM_t& Sc, const opts_t& opts,
562  WorkCompress<scalar_t>& w, int d, int depth);
563  void compute_U_basis_stable(DenseM_t& Sr, const opts_t& opts,
564  WorkCompress<scalar_t>& w,
565  int d, int dd, int depth);
566  void compute_V_basis_stable(DenseM_t& Sc, const opts_t& opts,
567  WorkCompress<scalar_t>& w,
568  int d, int dd, int depth);
569  void reduce_local_samples(DenseM_t& Rr, DenseM_t& Rc,
570  WorkCompress<scalar_t>& w,
571  int d0, int d, int depth);
572  bool update_orthogonal_basis(const opts_t& opts, scalar_t& r_max_0,
573  const DenseM_t& S, DenseM_t& Q,
574  int d, int dd, bool untouched,
575  int L, int depth);
576  void set_U_full_rank(WorkCompress<scalar_t>& w);
577  void set_V_full_rank(WorkCompress<scalar_t>& w);
578 
579  void compress_level_original(DenseM_t& Rr, DenseM_t& Rc,
580  DenseM_t& Sr, DenseM_t& Sc,
581  const opts_t& opts,
582  WorkCompress<scalar_t>& w,
583  int dd, int lvl, int depth) override;
584  void compress_level_stable(DenseM_t& Rr, DenseM_t& Rc,
585  DenseM_t& Sr, DenseM_t& Sc,
586  const opts_t& opts,
587  WorkCompress<scalar_t>& w,
588  int d, int dd, int lvl, int depth) override;
589  void get_extraction_indices(std::vector<std::vector<std::size_t>>& I,
590  std::vector<std::vector<std::size_t>>& J,
591  const std::pair<std::size_t,std::size_t>& off,
592  WorkCompress<scalar_t>& w,
593  int& self, int lvl) override;
594  void get_extraction_indices(std::vector<std::vector<std::size_t>>& I,
595  std::vector<std::vector<std::size_t>>& J,
596  std::vector<DenseM_t*>& B,
597  const std::pair<std::size_t,std::size_t>& off,
598  WorkCompress<scalar_t>& w,
599  int& self, int lvl) override;
600  void extract_D_B(const elem_t& Aelem, const opts_t& opts,
601  WorkCompress<scalar_t>& w, int lvl) override;
602 
603  void compress(const kernel::Kernel<real_t>& K, const opts_t& opts);
604  void compress_recursive_ann(DenseMatrix<std::uint32_t>& ann,
605  DenseMatrix<real_t>& scores,
606  const elem_t& Aelem, const opts_t& opts,
607  WorkCompressANN<scalar_t>& w,
608  int depth) override;
609  void compute_local_samples_ann(DenseMatrix<std::uint32_t>& ann,
610  DenseMatrix<real_t>& scores,
611  WorkCompressANN<scalar_t>& w,
612  const elem_t& Aelem, const opts_t& opts);
613  bool compute_U_V_bases_ann(DenseM_t& S, const opts_t& opts,
614  WorkCompressANN<scalar_t>& w, int depth);
615 
616  void factor_recursive(WorkFactor<scalar_t>& w,
617  bool isroot, bool partial,
618  int depth) override;
619 
620  void apply_fwd(const DenseM_t& b, WorkApply<scalar_t>& w,
621  bool isroot, int depth,
622  std::atomic<long long int>& flops) const override;
623  void apply_bwd(const DenseM_t& b, scalar_t beta, DenseM_t& c,
624  WorkApply<scalar_t>& w, bool isroot, int depth,
625  std::atomic<long long int>& flops) const override;
626  void applyT_fwd(const DenseM_t& b, WorkApply<scalar_t>& w, bool isroot,
627  int depth, std::atomic<long long int>& flops) const override;
628  void applyT_bwd(const DenseM_t& b, scalar_t beta, DenseM_t& c,
629  WorkApply<scalar_t>& w, bool isroot, int depth,
630  std::atomic<long long int>& flops) const override;
631 
632  void solve_fwd(const DenseM_t& b, WorkSolve<scalar_t>& w,
633  bool partial, bool isroot, int depth) const override;
634  void solve_bwd(DenseM_t& x, WorkSolve<scalar_t>& w,
635  bool isroot, int depth) const override;
636 
637  void extract_fwd(WorkExtract<scalar_t>& w,
638  bool odiag, int depth) const override;
639  void extract_bwd(DenseM_t& B, WorkExtract<scalar_t>& w,
640  int depth) const override;
641  void extract_bwd(std::vector<Triplet<scalar_t>>& triplets,
642  WorkExtract<scalar_t>& w, int depth) const override;
643  void extract_bwd_internal(WorkExtract<scalar_t>& w, int depth) const;
644 
645  void apply_UV_big(DenseM_t& Theta, DenseM_t& Uop, DenseM_t& Phi,
646  DenseM_t& Vop,
647  const std::pair<std::size_t, std::size_t>& offset,
648  int depth, std::atomic<long long int>& flops)
649  const override;
650  void apply_UtVt_big(const DenseM_t& A, DenseM_t& UtA, DenseM_t& VtA,
651  const std::pair<std::size_t, std::size_t>& offset,
652  int depth, std::atomic<long long int>& flops)
653  const override;
654 
655  void dense_recursive(DenseM_t& A, WorkDense<scalar_t>& w,
656  bool isroot, int depth) const override;
657 
661  template<typename T> friend
662  void apply_HSS(Trans ta, const HSSMatrix<T>& a, const DenseMatrix<T>& b,
663  T beta, DenseMatrix<T>& c);
664 
668  template<typename T> friend
669  void draw(const HSSMatrix<T>& H, const std::string& name);
670 
671  void read(std::ifstream& is) override;
672  void write(std::ofstream& os) const override;
673 
674  friend class HSSMatrixMPI<scalar_t>;
675 
676  using HSSMatrixBase<scalar_t>::child;
677  };
678 
686  template<typename scalar_t>
687  void draw(const HSSMatrix<scalar_t>& H, const std::string& name);
688 
700  template<typename scalar_t> void
701  apply_HSS(Trans op, const HSSMatrix<scalar_t>& A,
702  const DenseMatrix<scalar_t>& B,
703  scalar_t beta, DenseMatrix<scalar_t>& C);
704 
705  } // end namespace HSS
706 } // end namespace strumpack
707 
708 #endif // HSS_MATRIX_HPP
HSSMatrixBase.hpp
This file contains the HSSMatrixBase class definition, an abstract class for HSS matrix representatio...
HSSOptions.hpp
Contains the HSSOptions class as well as general routines for HSS options.
Kernel.hpp
Definitions of several kernel functions, and helper routines. Also provides driver routines for kerne...
strumpack::DenseMatrixWrapper
Like DenseMatrix, this class represents a matrix, stored in column major format, to allow direct use ...
Definition: DenseMatrix.hpp:1015
strumpack::DenseMatrix
This class represents a matrix, stored in column major format, to allow direct use of BLAS/LAPACK rou...
Definition: DenseMatrix.hpp:138
strumpack::HSS::HSSFactors
Contains data related to ULV factorization of an HSS matrix.
Definition: HSSExtra.hpp:161
strumpack::HSS::HSSMatrixBase
Abstract base class for Hierarchically Semi-Separable (HSS) matrices.
Definition: HSSMatrixBase.hpp:83
strumpack::HSS::HSSMatrixBase::active
bool active() const
Definition: HSSMatrixBase.hpp:239
strumpack::HSS::HSSMatrixMPI
Distributed memory implementation of the HSS (Hierarchically Semi-Separable) matrix format.
Definition: HSSMatrixMPI.hpp:68
strumpack::HSS::HSSMatrix
Class to represent a sequential/threaded Hierarchically Semi-Separable matrix.
Definition: HSSMatrix.hpp:80
strumpack::HSS::HSSMatrix::factor
void factor() override
strumpack::HSS::HSSMatrix::forward_solve
void forward_solve(WorkSolve< scalar_t > &w, const DenseM_t &b, bool partial) const override
strumpack::HSS::HSSMatrix::levels
std::size_t levels() const override
strumpack::HSS::HSSMatrix::print_info
void print_info(std::ostream &out=std::cout, std::size_t roff=0, std::size_t coff=0) const override
strumpack::HSS::HSSMatrix::child
const HSSMatrix< scalar_t > * child(int c) const
Definition: HSSMatrix.hpp:194
strumpack::HSS::HSSMatrix::HSSMatrix
HSSMatrix(const DenseM_t &A, const opts_t &opts)
strumpack::HSS::HSSMatrix::draw
friend void draw(const HSSMatrix< T > &H, const std::string &name)
strumpack::HSS::HSSMatrix::draw
void draw(std::ostream &of, std::size_t rlo=0, std::size_t clo=0) const override
strumpack::HSS::HSSMatrix::extract_add
void extract_add(const std::vector< std::size_t > &I, const std::vector< std::size_t > &J, DenseM_t &B) const
strumpack::HSS::HSSMatrix::get
scalar_t get(std::size_t i, std::size_t j) const
strumpack::HSS::HSSMatrix::shift
void shift(scalar_t sigma) override
strumpack::HSS::HSSMatrix::compress_with_coordinates
void compress_with_coordinates(const DenseMatrix< real_t > &coords, const std::function< void(const std::vector< std::size_t > &I, const std::vector< std::size_t > &J, DenseM_t &B)> &Aelem, const opts_t &opts)
strumpack::HSS::HSSMatrix::clone
std::unique_ptr< HSSMatrixBase< scalar_t > > clone() const override
strumpack::HSS::HSSMatrix::compress
void compress(const std::function< void(DenseM_t &Rr, DenseM_t &Rc, DenseM_t &Sr, DenseM_t &Sc)> &Amult, const std::function< void(const std::vector< std::size_t > &I, const std::vector< std::size_t > &J, DenseM_t &B)> &Aelem, const opts_t &opts)
strumpack::HSS::HSSMatrix::memory
std::size_t memory() const override
strumpack::HSS::HSSMatrix::applyC
DenseM_t applyC(const DenseM_t &b) const
strumpack::HSS::HSSMatrix::HSSMatrix
HSSMatrix(HSSMatrix< scalar_t > &&other)=default
strumpack::HSS::HSSMatrix::write
void write(const std::string &fname) const
strumpack::HSS::HSSMatrix::read
static HSSMatrix< scalar_t > read(const std::string &fname)
strumpack::HSS::HSSMatrix::HSSMatrix
HSSMatrix(const HSSMatrix< scalar_t > &other)
strumpack::HSS::HSSMatrix::HSSMatrix
HSSMatrix(kernel::Kernel< real_t > &K, const opts_t &opts)
strumpack::HSS::HSSMatrix::child
HSSMatrix< scalar_t > * child(int c)
Definition: HSSMatrix.hpp:204
strumpack::HSS::HSSMatrix::backward_solve
void backward_solve(WorkSolve< scalar_t > &w, DenseM_t &x) const override
strumpack::HSS::HSSMatrix::apply_HSS
friend void apply_HSS(Trans ta, const HSSMatrix< T > &a, const DenseMatrix< T > &b, T beta, DenseMatrix< T > &c)
strumpack::HSS::HSSMatrix::operator=
HSSMatrix< scalar_t > & operator=(const HSSMatrix< scalar_t > &other)
strumpack::HSS::HSSMatrix::extract
DenseM_t extract(const std::vector< std::size_t > &I, const std::vector< std::size_t > &J) const
strumpack::HSS::HSSMatrix::solve
void solve(DenseM_t &b) const override
strumpack::HSS::HSSMatrix::compress
void compress(const DenseM_t &A, const opts_t &opts)
strumpack::HSS::HSSMatrix::HSSMatrix
HSSMatrix(std::size_t m, std::size_t n, const opts_t &opts)
strumpack::HSS::HSSMatrix::rank
std::size_t rank() const override
strumpack::HSS::HSSMatrix::HSSMatrix
HSSMatrix(const structured::ClusterTree &t, const opts_t &opts)
strumpack::HSS::HSSMatrix::dense
DenseM_t dense() const
strumpack::HSS::HSSMatrix::mult
void mult(Trans op, const DenseM_t &x, DenseM_t &y) const override
strumpack::HSS::HSSMatrix::operator=
HSSMatrix< scalar_t > & operator=(HSSMatrix< scalar_t > &&other)=default
strumpack::HSS::HSSMatrix::nonzeros
std::size_t nonzeros() const override
strumpack::HSS::HSSMatrix::apply
DenseM_t apply(const DenseM_t &b) const
strumpack::HSS::HSSOptions
Class containing several options for the HSS code and data-structures.
Definition: HSSOptions.hpp:152
strumpack::kernel::Kernel
Representation of a kernel matrix.
Definition: Kernel.hpp:73
strumpack::structured::ClusterTree
The cluster tree, or partition tree that represents the partitioning of the rows or columns of a hier...
Definition: ClusterTree.hpp:67
strumpack::HSS::apply_HSS
void apply_HSS(Trans op, const HSSMatrix< scalar_t > &A, const DenseMatrix< scalar_t > &B, scalar_t beta, DenseMatrix< scalar_t > &C)
strumpack::HSS::draw
void draw(const HSSMatrix< scalar_t > &H, const std::string &name)
strumpack
Definition: StrumpackOptions.hpp:43
strumpack::Trans
Trans
Definition: DenseMatrix.hpp:51