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  // implement and test, separate from factorization
88  // BLRMatrix(std::size_t n, const std::vector<std::size_t>& tiles,
89  // const adm_t& admissible, const elem_t& Aelem,
90  // std::vector<int>& piv, const Opts_t& opts);
91 
92  BLRMatrix(DenseM_t& A, const std::vector<std::size_t>& tiles,
93  const adm_t& admissible, std::vector<int>& piv,
94  const Opts_t& opts);
95 
96  BLRMatrix(std::size_t m, const std::vector<std::size_t>& rowtiles,
97  std::size_t n, const std::vector<std::size_t>& coltiles);
98 
99  std::size_t rows() const override { return m_; }
100  std::size_t cols() const override { return n_; }
101 
102  std::size_t memory() const override;
103  std::size_t nonzeros() const override;
104  std::size_t rank() const override;
105 
106  DenseM_t dense() const;
107  void dense(DenseM_t& A) const;
108 
109  void draw(std::ostream& of, std::size_t roff, std::size_t coff) const;
110 
111  void print(const std::string& name) const;
112 
113  void clear();
114 
115  void solve(const std::vector<int>& P, DenseM_t& x) const {
116  x.laswp(P, true);
117  trsm(Side::L, UpLo::L, Trans::N, Diag::U, scalar_t(1.), *this, x, 0);
118  trsm(Side::L, UpLo::U, Trans::N, Diag::N, scalar_t(1.), *this, x, 0);
119  }
120 
132  void mult(Trans op, const DenseM_t& x, DenseM_t& y) const override;
133 
134  std::size_t rg2t(std::size_t i) const;
135  std::size_t cg2t(std::size_t j) const;
136 
137  scalar_t operator()(std::size_t i, std::size_t j) const;
138  scalar_t& operator()(std::size_t i, std::size_t j);
139  DenseM_t extract(const std::vector<std::size_t>& I,
140  const std::vector<std::size_t>& J) const;
141 
142  void decompress();
143  void decompress_local_columns(int c_min, int c_max);
144  void remove_tiles_before_local_column(int c_min, int c_max);
145 
146  std::size_t rowblocks() const { return nbrows_; }
147  std::size_t colblocks() const { return nbcols_; }
148  std::size_t tilerows(std::size_t i) const { return roff_[i+1] - roff_[i]; }
149  std::size_t tilecols(std::size_t j) const { return coff_[j+1] - coff_[j]; }
150  std::size_t tileroff(std::size_t i) const { return roff_[i]; }
151  std::size_t tilecoff(std::size_t j) const { return coff_[j]; }
152 
153  BLRTile<scalar_t>& tile(std::size_t i, std::size_t j);
154  const BLRTile<scalar_t>& tile(std::size_t i, std::size_t j) const;
155  std::unique_ptr<BLRTile<scalar_t>>& block(std::size_t i, std::size_t j);
156  DenseMW_t tile(DenseM_t& A, std::size_t i, std::size_t j) const;
157  DenseTile<scalar_t>& tile_dense(std::size_t i, std::size_t j);
158  const DenseTile<scalar_t>& tile_dense(std::size_t i, std::size_t j) const;
159 
160  void compress_tile(std::size_t i, std::size_t j, const Opts_t& opts);
161  void fill(scalar_t v);
162  void fill_col(scalar_t v, std::size_t k, std::size_t CP);
163 
164  static void
165  construct_and_partial_factor(DenseM_t& A11, DenseM_t& A12,
166  DenseM_t& A21, DenseM_t& A22,
167  BLRMatrix<scalar_t>& B11,
168  std::vector<int>& piv,
169  BLRMatrix<scalar_t>& B12,
170  BLRMatrix<scalar_t>& B21,
171  const std::vector<std::size_t>& tiles1,
172  const std::vector<std::size_t>& tiles2,
173  const adm_t& admissible,
174  const Opts_t& opts);
175 
176  static void
177  construct_and_partial_factor(BLRMatrix<scalar_t>& B11,
178  BLRMatrix<scalar_t>& B12,
179  BLRMatrix<scalar_t>& B21,
180  BLRMatrix<scalar_t>& B22,
181  std::vector<int>& piv,
182  const std::vector<std::size_t>& tiles1,
183  const std::vector<std::size_t>& tiles2,
184  const adm_t& admissible,
185  const Opts_t& opts);
186 
187  static void
188  construct_and_partial_factor_col(BLRMatrix<scalar_t>& B11,
189  BLRMatrix<scalar_t>& B12,
190  BLRMatrix<scalar_t>& B21,
191  BLRMatrix<scalar_t>& B22,
192  std::vector<int>& piv,
193  const std::vector<std::size_t>& tiles1,
194  const std::vector<std::size_t>& tiles2,
195  const adm_t& admissible,
196  const Opts_t& opts,
197  const std::function<void
198  (int, bool, std::size_t)>& blockcol);
199 
200  static void
201  construct_and_partial_factor(std::size_t n1, std::size_t n2,
202  const extract_t<scalar_t>& A11,
203  const extract_t<scalar_t>& A12,
204  const extract_t<scalar_t>& A21,
205  const extract_t<scalar_t>& A22,
206  BLRMatrix<scalar_t>& B11,
207  std::vector<int>& piv,
208  BLRMatrix<scalar_t>& B12,
209  BLRMatrix<scalar_t>& B21,
210  BLRMatrix<scalar_t>& B22,
211  const std::vector<std::size_t>& tiles1,
212  const std::vector<std::size_t>& tiles2,
213  const adm_t& admissible,
214  const BLROptions<scalar_t>& opts);
215 
216  static void
217  trsmLNU_gemm(const BLRMatrix<scalar_t>& F1,
218  const BLRMatrix<scalar_t>& F2,
219  DenseM_t& B1, DenseM_t& B2, int task_depth);
220 
221  static void
222  gemm_trsmUNN(const BLRMatrix<scalar_t>& F1,
223  const BLRMatrix<scalar_t>& F2,
224  DenseM_t& B1, DenseM_t& B2, int task_depth);
225 
226  private:
227  std::size_t m_ = 0, n_ = 0, nbrows_ = 0, nbcols_ = 0;
228  std::vector<std::size_t> roff_, coff_, cl2l_, rl2l_;
229  std::vector<std::unique_ptr<BLRTile<scalar_t>>> blocks_;
230 
231  void create_dense_tile(std::size_t i, std::size_t j, DenseM_t& A);
232  void create_dense_tile(std::size_t i, std::size_t j,
233  const extract_t<scalar_t>& Aelem);
234  void create_dense_tile_left_looking(std::size_t i, std::size_t j,
235  const extract_t<scalar_t>& Aelem);
236  void create_dense_tile_left_looking(std::size_t i, std::size_t j,
237  std::size_t k,
238  const extract_t<scalar_t>& Aelem,
239  const BLRMatrix<scalar_t>& B21,
240  const BLRMatrix<scalar_t>& B12);
241  void create_LR_tile(std::size_t i, std::size_t j,
242  DenseM_t& A, const Opts_t& opts);
243  void create_LR_tile_left_looking(std::size_t i, std::size_t j,
244  const extract_t<scalar_t>& Aelem,
245  const Opts_t& opts);
246 
247  void create_LR_tile_left_looking(std::size_t i, std::size_t j,
248  std::size_t k,
249  const extract_t<scalar_t>& Aelem,
250  const BLRMatrix<scalar_t>& B21,
251  const BLRMatrix<scalar_t>& B12,
252  const Opts_t& opts);
253 
254  void LUAR_B11(std::size_t i, std::size_t j, std::size_t kmax,
255  DenseM_t& A11, const Opts_t& opts, int* B);
256  void LUAR_B12(std::size_t i, std::size_t j, std::size_t kmax,
257  BLRMatrix<scalar_t>& B11, DenseM_t& A12,
258  const Opts_t& opts, int* B);
259  void LUAR_B21(std::size_t i, std::size_t j, std::size_t kmax,
260  BLRMatrix<scalar_t>& B11, DenseM_t& A21,
261  const Opts_t& opts, int* B);
262 
263  template<typename T> friend
264  void draw(const BLRMatrix<T>& H, const std::string& name);
265  template<typename T,typename I> friend class BLRExtendAdd;
266  };
267 
268  template<typename scalar_t> void
269  LUAR(const std::vector<BLRTile<scalar_t>*>& Ti,
270  const std::vector<BLRTile<scalar_t>*>& Tj,
271  DenseMatrixWrapper<scalar_t>& tij,
272  const BLROptions<scalar_t>& opts, int* B);
273 
274  template<typename scalar_t> void
275  LUAR_B22(std::size_t i, std::size_t j, std::size_t kmax,
276  BLRMatrix<scalar_t>& B12, BLRMatrix<scalar_t>& B21,
277  DenseMatrix<scalar_t>& A22,
278  const BLROptions<scalar_t>& opts, int* B);
279 
280  template<typename scalar_t> void
281  trsm(Side s, UpLo ul, Trans ta, Diag d, scalar_t alpha,
282  const BLRMatrix<scalar_t>& a, DenseMatrix<scalar_t>& b,
283  int task_depth);
284 
285  template<typename scalar_t> void
286  trsv(UpLo ul, Trans ta, Diag d, const BLRMatrix<scalar_t>& a,
287  DenseMatrix<scalar_t>& b, int task_depth);
288 
289  template<typename scalar_t> void
290  gemv(Trans ta, scalar_t alpha, const BLRMatrix<scalar_t>& a,
291  const DenseMatrix<scalar_t>& x, scalar_t beta,
292  DenseMatrix<scalar_t>& y, int task_depth);
293 
294  template<typename scalar_t> void
295  gemm(Trans ta, Trans tb, scalar_t alpha, const BLRMatrix<scalar_t>& a,
296  const BLRMatrix<scalar_t>& b, scalar_t beta,
297  DenseMatrix<scalar_t>& c, int task_depth);
298 
299  template<typename scalar_t> void
300  gemm(Trans ta, Trans tb, scalar_t alpha, const BLRMatrix<scalar_t>& A,
301  const DenseMatrix<scalar_t>& B, scalar_t beta,
302  DenseMatrix<scalar_t>& C, int task_depth);
303 
304  template<typename scalar_t>
305  void draw(const BLRMatrix<scalar_t>& B, const std::string& name) {
306  std::ofstream of("plot" + name + ".gnuplot");
307  of << "set terminal pdf enhanced color size 5,4" << std::endl;
308  of << "set output '" << name << ".pdf'" << std::endl;
309  B.draw(of, 0, 0);
310  of << "set xrange [0:" << B.cols() << "]" << std::endl;
311  of << "set yrange [" << B.rows() << ":0]" << std::endl;
312  of << "plot x lt -1 notitle" << std::endl;
313  of.close();
314  }
315 
316  } // end namespace BLR
317 } // end namespace strumpack
318 
319 #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:99
strumpack::BLR::BLRMatrix::cols
std::size_t cols() const override
Definition: BLRMatrix.hpp:100
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::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:42
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