BLRMatrix.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  */
32 #ifndef BLR_MATRIX_HPP
33 #define BLR_MATRIX_HPP
34 
35 #include <cassert>
36 #include <memory>
37 #include <functional>
38 #include <algorithm>
39 
40 #include "BLROptions.hpp"
41 #include "BLRTileBLAS.hpp" // TODO remove
42 #include "structured/StructuredMatrix.hpp"
43 
44 namespace strumpack {
45  namespace BLR {
46 
47  // forward declarations
48  template<typename scalar> class BLRTile;
49  template<typename scalar_t,typename integer_t> class BLRExtendAdd;
50 
51 
52  template<typename T>
53  using extract_t =
54  std::function<void(const std::vector<std::size_t>&,
55  const std::vector<std::size_t>&,
56  DenseMatrix<T>&)>;
57  using adm_t = DenseMatrix<bool>;
58 
73  template<typename scalar_t> class BLRMatrix
74  : public structured::StructuredMatrix<scalar_t> {
75  using DenseM_t = DenseMatrix<scalar_t>;
76  using DenseMW_t = DenseMatrixWrapper<scalar_t>;
77  using Opts_t = BLROptions<scalar_t>;
78 
79  public:
80  BLRMatrix() = default;
81 
82  BLRMatrix(DenseM_t& A,
83  const std::vector<std::size_t>& rowtiles,
84  const std::vector<std::size_t>& coltiles,
85  const Opts_t& opts);
86 
87  BLRMatrix(DenseM_t& A, const std::vector<std::size_t>& tiles,
88  const adm_t& admissible, const Opts_t& opts);
89 
90  BLRMatrix(std::size_t m, const std::vector<std::size_t>& rowtiles,
91  std::size_t n, const std::vector<std::size_t>& coltiles);
92 
93  std::size_t rows() const override { return m_; }
94  std::size_t cols() const override { return n_; }
95 
96  std::size_t memory() const override;
97  std::size_t nonzeros() const override;
98  std::size_t rank() const override;
99 
100  DenseM_t dense() const;
101  void dense(DenseM_t& A) const;
102 
103  void draw(std::ostream& of, std::size_t roff, std::size_t coff) const;
104 
105  void print(const std::string& name) const;
106 
107  void clear();
108 
109  void solve(DenseM_t& x) const override {
110  x.laswp(piv_, true);
111  trsm(Side::L, UpLo::L, Trans::N, Diag::U, scalar_t(1.), *this, x, 0);
112  trsm(Side::L, UpLo::U, Trans::N, Diag::N, scalar_t(1.), *this, x, 0);
113  }
114 
115  const std::vector<int>& piv() const { return piv_; }
116 
128  void mult(Trans op, const DenseM_t& x, DenseM_t& y) const override;
129 
130  std::size_t rg2t(std::size_t i) const;
131  std::size_t cg2t(std::size_t j) const;
132 
133  scalar_t operator()(std::size_t i, std::size_t j) const;
134  scalar_t& operator()(std::size_t i, std::size_t j);
135  DenseM_t extract(const std::vector<std::size_t>& I,
136  const std::vector<std::size_t>& J) const;
137 
138  void decompress();
139  void decompress_local_columns(int c_min, int c_max);
140  void remove_tiles_before_local_column(int c_min, int c_max);
141 
142  std::size_t rowblocks() const { return nbrows_; }
143  std::size_t colblocks() const { return nbcols_; }
144  std::size_t tilerows(std::size_t i) const { return roff_[i+1] - roff_[i]; }
145  std::size_t tilecols(std::size_t j) const { return coff_[j+1] - coff_[j]; }
146  std::size_t tileroff(std::size_t i) const { return roff_[i]; }
147  std::size_t tilecoff(std::size_t j) const { return coff_[j]; }
148 
149  BLRTile<scalar_t>& tile(std::size_t i, std::size_t j);
150  const BLRTile<scalar_t>& tile(std::size_t i, std::size_t j) const;
151  std::unique_ptr<BLRTile<scalar_t>>& block(std::size_t i, std::size_t j);
152  DenseMW_t tile(DenseM_t& A, std::size_t i, std::size_t j) const;
153  DenseTile<scalar_t>& tile_dense(std::size_t i, std::size_t j);
154  const DenseTile<scalar_t>& tile_dense(std::size_t i, std::size_t j) const;
155 
156  void compress_tile(std::size_t i, std::size_t j, const Opts_t& opts);
157  void fill(scalar_t v);
158  void fill_col(scalar_t v, std::size_t k, std::size_t CP);
159 
160  static void
161  construct_and_partial_factor(DenseM_t& A11, DenseM_t& A12,
162  DenseM_t& A21, DenseM_t& A22,
163  BLRMatrix<scalar_t>& B11,
164  BLRMatrix<scalar_t>& B12,
165  BLRMatrix<scalar_t>& B21,
166  const std::vector<std::size_t>& tiles1,
167  const std::vector<std::size_t>& tiles2,
168  const adm_t& admissible,
169  const Opts_t& opts);
170 
171  static void
172  construct_and_partial_factor(BLRMatrix<scalar_t>& B11,
173  BLRMatrix<scalar_t>& B12,
174  BLRMatrix<scalar_t>& B21,
175  BLRMatrix<scalar_t>& B22,
176  const std::vector<std::size_t>& tiles1,
177  const std::vector<std::size_t>& tiles2,
178  const adm_t& admissible,
179  const Opts_t& opts);
180 
181  static void
182  construct_and_partial_factor_col(BLRMatrix<scalar_t>& B11,
183  BLRMatrix<scalar_t>& B12,
184  BLRMatrix<scalar_t>& B21,
185  BLRMatrix<scalar_t>& B22,
186  const std::vector<std::size_t>& tiles1,
187  const std::vector<std::size_t>& tiles2,
188  const adm_t& admissible,
189  const Opts_t& opts,
190  const std::function<void
191  (int, bool, std::size_t)>& blockcol);
192 
193  static void
194  construct_and_partial_factor(std::size_t n1, std::size_t n2,
195  const extract_t<scalar_t>& A11,
196  const extract_t<scalar_t>& A12,
197  const extract_t<scalar_t>& A21,
198  const extract_t<scalar_t>& A22,
199  BLRMatrix<scalar_t>& B11,
200  BLRMatrix<scalar_t>& B12,
201  BLRMatrix<scalar_t>& B21,
202  BLRMatrix<scalar_t>& B22,
203  const std::vector<std::size_t>& tiles1,
204  const std::vector<std::size_t>& tiles2,
205  const adm_t& admissible,
206  const BLROptions<scalar_t>& opts);
207 
208  static void
209  trsmLNU_gemm(const BLRMatrix<scalar_t>& F1,
210  const BLRMatrix<scalar_t>& F2,
211  DenseM_t& B1, DenseM_t& B2, int task_depth);
212 
213  static void
214  gemm_trsmUNN(const BLRMatrix<scalar_t>& F1,
215  const BLRMatrix<scalar_t>& F2,
216  DenseM_t& B1, DenseM_t& B2, int task_depth);
217 
218  private:
219  std::size_t m_ = 0, n_ = 0, nbrows_ = 0, nbcols_ = 0;
220  std::vector<std::size_t> roff_, coff_, cl2l_, rl2l_;
221  std::vector<std::unique_ptr<BLRTile<scalar_t>>> blocks_;
222  std::vector<int> piv_;
223 
224  void create_dense_tile(std::size_t i, std::size_t j, DenseM_t& A);
225  void create_dense_tile(std::size_t i, std::size_t j,
226  const extract_t<scalar_t>& Aelem);
227  void create_dense_tile_left_looking(std::size_t i, std::size_t j,
228  const extract_t<scalar_t>& Aelem);
229  void create_dense_tile_left_looking(std::size_t i, std::size_t j,
230  std::size_t k,
231  const extract_t<scalar_t>& Aelem,
232  const BLRMatrix<scalar_t>& B21,
233  const BLRMatrix<scalar_t>& B12);
234  void create_LR_tile(std::size_t i, std::size_t j,
235  DenseM_t& A, const Opts_t& opts);
236  void create_LR_tile_left_looking(std::size_t i, std::size_t j,
237  const extract_t<scalar_t>& Aelem,
238  const Opts_t& opts);
239 
240  void create_LR_tile_left_looking(std::size_t i, std::size_t j,
241  std::size_t k,
242  const extract_t<scalar_t>& Aelem,
243  const BLRMatrix<scalar_t>& B21,
244  const BLRMatrix<scalar_t>& B12,
245  const Opts_t& opts);
246 
247  void LUAR_B11(std::size_t i, std::size_t j, std::size_t kmax,
248  DenseM_t& A11, const Opts_t& opts, int* B);
249  void LUAR_B12(std::size_t i, std::size_t j, std::size_t kmax,
250  BLRMatrix<scalar_t>& B11, DenseM_t& A12,
251  const Opts_t& opts, int* B);
252  void LUAR_B21(std::size_t i, std::size_t j, std::size_t kmax,
253  BLRMatrix<scalar_t>& B11, DenseM_t& A21,
254  const Opts_t& opts, int* B);
255 
256  template<typename T> friend
257  void draw(const BLRMatrix<T>& H, const std::string& name);
258  template<typename T,typename I> friend class BLRExtendAdd;
259  };
260 
261  template<typename scalar_t> void
262  LUAR(const std::vector<BLRTile<scalar_t>*>& Ti,
263  const std::vector<BLRTile<scalar_t>*>& Tj,
264  DenseMatrixWrapper<scalar_t>& tij,
265  const BLROptions<scalar_t>& opts, int* B);
266 
267  template<typename scalar_t> void
268  LUAR_B22(std::size_t i, std::size_t j, std::size_t kmax,
269  BLRMatrix<scalar_t>& B12, BLRMatrix<scalar_t>& B21,
270  DenseMatrix<scalar_t>& A22,
271  const BLROptions<scalar_t>& opts, int* B);
272 
273  template<typename scalar_t> void
274  trsm(Side s, UpLo ul, Trans ta, Diag d, scalar_t alpha,
275  const BLRMatrix<scalar_t>& a, DenseMatrix<scalar_t>& b,
276  int task_depth);
277 
278  template<typename scalar_t> void
279  trsv(UpLo ul, Trans ta, Diag d, const BLRMatrix<scalar_t>& a,
280  DenseMatrix<scalar_t>& b, int task_depth);
281 
282  template<typename scalar_t> void
283  gemv(Trans ta, scalar_t alpha, const BLRMatrix<scalar_t>& a,
284  const DenseMatrix<scalar_t>& x, scalar_t beta,
285  DenseMatrix<scalar_t>& y, int task_depth);
286 
287  template<typename scalar_t> void
288  gemm(Trans ta, Trans tb, scalar_t alpha, const BLRMatrix<scalar_t>& a,
289  const BLRMatrix<scalar_t>& b, scalar_t beta,
290  DenseMatrix<scalar_t>& c, int task_depth);
291 
292  template<typename scalar_t> void
293  gemm(Trans ta, Trans tb, scalar_t alpha, const BLRMatrix<scalar_t>& A,
294  const DenseMatrix<scalar_t>& B, scalar_t beta,
295  DenseMatrix<scalar_t>& C, int task_depth);
296 
297  template<typename scalar_t>
298  void draw(const BLRMatrix<scalar_t>& B, const std::string& name) {
299  std::ofstream of("plot" + name + ".gnuplot");
300  of << "set terminal pdf enhanced color size 5,4" << std::endl;
301  of << "set output '" << name << ".pdf'" << std::endl;
302  B.draw(of, 0, 0);
303  of << "set xrange [0:" << B.cols() << "]" << std::endl;
304  of << "set yrange [" << B.rows() << ":0]" << std::endl;
305  of << "plot x lt -1 notitle" << std::endl;
306  of.close();
307  }
308 
309  } // end namespace BLR
310 } // end namespace strumpack
311 
312 #endif // BLR_MATRIX_HPP
BLROptions.hpp
Contains class holding BLROptions.
StructuredMatrix.hpp
Contains the structured matrix interfaces.
strumpack::BLR::BLRExtendAdd
Definition: BLRMatrix.hpp:49
strumpack::BLR::BLRMatrix
Class to represent a block low-rank matrix.
Definition: BLRMatrix.hpp:74
strumpack::BLR::BLRMatrix::rows
std::size_t rows() const override
Definition: BLRMatrix.hpp:93
strumpack::BLR::BLRMatrix::cols
std::size_t cols() const override
Definition: BLRMatrix.hpp:94
strumpack::BLR::BLRMatrix::nonzeros
std::size_t nonzeros() const override
strumpack::BLR::BLRMatrix::memory
std::size_t memory() const override
strumpack::BLR::BLRMatrix::rank
std::size_t rank() const override
strumpack::BLR::BLRMatrix::solve
void solve(DenseM_t &x) const override
Definition: BLRMatrix.hpp:109
strumpack::BLR::BLRMatrix::mult
void mult(Trans op, const DenseM_t &x, DenseM_t &y) const override
strumpack::BLR::BLROptions
Class containing several options for the BLR code and data-structures.
Definition: BLROptions.hpp:82
strumpack::BLR::BLRTile
Definition: BLRMatrix.hpp:48
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::DenseMatrix::laswp
void laswp(const std::vector< int > &P, bool fwd)
strumpack::structured::StructuredMatrix
Class to represent a structured matrix. This is the abstract base class for several types of structur...
Definition: StructuredMatrix.hpp:209
strumpack
Definition: StrumpackOptions.hpp:43
strumpack::UpLo
UpLo
Definition: DenseMatrix.hpp:83
strumpack::trsv
void trsv(UpLo ul, Trans ta, Diag d, const DenseMatrix< scalar_t > &a, DenseMatrix< scalar_t > &b, int depth=0)
strumpack::gemv
void gemv(Trans ta, scalar_t alpha, const DenseMatrix< scalar_t > &a, const DenseMatrix< scalar_t > &x, scalar_t beta, DenseMatrix< scalar_t > &y, int depth=0)
strumpack::Trans
Trans
Definition: DenseMatrix.hpp:51
strumpack::trsm
void trsm(Side s, UpLo ul, Trans ta, Diag d, scalar_t alpha, const DenseMatrix< scalar_t > &a, DenseMatrix< scalar_t > &b, int depth=0)
strumpack::Side
Side
Definition: DenseMatrix.hpp:74
strumpack::gemm
void gemm(Trans ta, Trans tb, scalar_t alpha, const DenseMatrix< scalar_t > &a, const DenseMatrix< scalar_t > &b, scalar_t beta, DenseMatrix< scalar_t > &c, int depth=0)
strumpack::Diag
Diag
Definition: DenseMatrix.hpp:92