CompressedSparseMatrix.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  */
34 #ifndef COMPRESSED_SPARSE_MATRIX_HPP
35 #define COMPRESSED_SPARSE_MATRIX_HPP
36 
37 #include <vector>
38 #include <string>
39 #include <tuple>
40 
41 #include "misc/Tools.hpp"
42 #include "misc/Triplet.hpp"
43 #include "dense/DenseMatrix.hpp"
44 #include "StrumpackOptions.hpp"
45 
46 // where is this used?? in MC64?
47 #ifdef _LONGINT
48  typedef long long int int_t;
49 #else // Default
50  typedef int int_t;
51 #endif
52 
53 namespace strumpack {
54 
55  // forward declarations
56  template<typename integer_t> class CSRGraph;
57  template<typename scalar_t> class DenseMatrix;
58  template<typename scalar_t> class DistributedMatrix;
59 
60  extern "C" {
61  int_t strumpack_mc64id_(int_t*);
62  int_t strumpack_mc64ad_
63  (int_t*, int_t*, int_t*, int_t*, int_t*, double*, int_t*,
64  int_t*, int_t*, int_t*, int_t*, double*, int_t*, int_t*);
65  }
66 
67  template<typename scalar_t, typename integer_t,
68  typename real_t = typename RealType<scalar_t>::value_type>
69  class MatchingData {
70  public:
71  MatchingData() {}
72  MatchingData(MatchingJob j, std::size_t n) : job(j) {
73  if (job != MatchingJob::NONE)
74  Q.resize(n);
75  if (job == MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING) {
76  R.resize(n);
77  C.resize(n);
78  }
79  }
80 
81  MatchingJob job = MatchingJob::NONE;
82  std::vector<integer_t> Q;
83  std::vector<real_t> R, C;
84 
85  int_t mc64_work_int(std::size_t n, std::size_t nnz) const {
86  switch (job) {
87  case MatchingJob::MAX_SMALLEST_DIAGONAL: return 4*n; break;
88  case MatchingJob::MAX_SMALLEST_DIAGONAL_2: return 10*n + nnz; break;
89  case MatchingJob::MAX_CARDINALITY:
90  case MatchingJob::MAX_DIAGONAL_SUM:
91  case MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING:
92  default: return 5*n;
93  }
94  }
95 
96  int_t mc64_work_double(std::size_t n, std::size_t nnz) const {
97  switch (job) {
98  case MatchingJob::MAX_CARDINALITY: return 0; break;
99  case MatchingJob::MAX_SMALLEST_DIAGONAL: return n; break;
100  case MatchingJob::MAX_SMALLEST_DIAGONAL_2: return nnz; break;
101  case MatchingJob::MAX_DIAGONAL_SUM: return 2*n + nnz; break;
102  case MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING:
103  default: return 3*n + nnz; break;
104  }
105  }
106  };
107 
108  template<typename scalar_t,
109  typename real_t = typename RealType<scalar_t>::value_type>
110  class Equilibration {
111  public:
112  Equilibration() {}
113  Equilibration(std::size_t N) : R(N), C(N) {}
114  // use this one for block-row distributed sparse matrix, with rows
115  // local rows, and cols global columns
116  Equilibration(std::size_t rows, std::size_t cols) : R(rows), C(cols) {}
117 
118  int info = 0;
119  EquilibrationType type = EquilibrationType::NONE;
120  real_t rcond = 1, ccond = 1, Amax = 0;
121  std::vector<real_t> R, C;
122 
123  void set_type() {
124  const real_t thres = 0.1;
125  const real_t small = blas::lamch<real_t>('S') / blas::lamch<real_t>('P');
126  const real_t large = 1. / small;
127  if (rcond >= thres && Amax >= small && Amax <= large) {
128  R.clear();
129  if (ccond >= thres) {
130  type = EquilibrationType::NONE;
131  C.clear();
132  } else type = EquilibrationType::COLUMN;
133  } else {
134  if (ccond >= thres) {
135  type = EquilibrationType::ROW;
136  C.clear();
137  } else type = EquilibrationType::BOTH;
138  }
139  }
140  };
141 
162  template<typename scalar_t,typename integer_t>
163  class CompressedSparseMatrix {
164  using DenseM_t = DenseMatrix<scalar_t>;
165  using real_t = typename RealType<scalar_t>::value_type;
166 #if defined(STRUMPACK_USE_MPI)
167  using DistM_t = DistributedMatrix<scalar_t>;
168 #endif
169  using Match_t = MatchingData<scalar_t,integer_t>;
170  using Equil_t = Equilibration<scalar_t>;
171 
172  public:
176  virtual ~CompressedSparseMatrix() {}
177 
183  integer_t size() const { return n_; }
184 
190  integer_t nnz() const { return nnz_; }
191 
201  const integer_t* ptr() const { return ptr_.data(); }
202 
210  const integer_t* ind() const { return ind_.data(); }
211 
219  const scalar_t* val() const { return val_.data(); }
220 
229  integer_t* ptr() { return ptr_.data(); }
230 
238  integer_t* ind() { return ind_.data(); }
239 
246  scalar_t* val() { return val_.data(); }
247 
253  const integer_t& ptr(integer_t i) const { assert(i >= 0 && i <= size()); return ptr_[i]; }
254 
261  const integer_t& ind(integer_t i) const { assert(i >= 0 && i < nnz()); return ind_[i]; }
262 
268  const scalar_t& val(integer_t i) const { assert(i >= 0 && i < nnz()); return val_[i]; }
269 
275  integer_t& ptr(integer_t i) { assert(i <= size()); return ptr_[i]; }
276 
283  integer_t& ind(integer_t i) { assert(i < nnz()); return ind_[i]; }
284 
290  scalar_t& val(integer_t i) { assert(i < nnz()); return val_[i]; }
291 
296  bool symm_sparse() const { return symm_sparse_; }
297 
305  void set_symm_sparse(bool symm_sparse=true) { symm_sparse_ = symm_sparse; }
306 
307 
321  virtual void spmv(const DenseM_t& x, DenseM_t& y) const = 0;
322 
336  virtual void spmv(const scalar_t* x, scalar_t* y) const = 0;
337 
342  virtual void permute(const integer_t* iorder, const integer_t* order);
343 
344  virtual void permute(const std::vector<integer_t>& iorder,
345  const std::vector<integer_t>& order) {
346  permute(iorder.data(), order.data());
347  }
348 
349  virtual void permute_columns(const std::vector<integer_t>& perm) = 0;
350 
351  virtual Equil_t equilibration() const { return Equil_t(this->size()); }
352 
353  virtual void equilibrate(const Equil_t&) {}
354 
355  virtual Match_t matching(MatchingJob, bool apply=true);
356 
357  virtual void apply_matching(const Match_t&);
358 
359  virtual void symmetrize_sparsity();
360 
361  virtual void print() const;
362  virtual void print_dense(const std::string& name) const {
363  std::cerr << "print_dense not implemented for this matrix type"
364  << std::endl;
365  }
366  virtual void print_matrix_market(const std::string& filename) const {
367  std::cerr << "print_matrix_market not implemented for this matrix type"
368  << std::endl;
369  }
370 
371  virtual int read_matrix_market(const std::string& filename) = 0;
372 
373  virtual real_t max_scaled_residual
374  (const scalar_t* x, const scalar_t* b) const = 0;
375  virtual real_t max_scaled_residual
376  (const DenseM_t& x, const DenseM_t& b) const = 0;
377 
378 
379 #ifndef DOXYGEN_SHOULD_SKIP_THIS
380  virtual CSRGraph<integer_t> extract_graph
381  (int ordering_level, integer_t lo, integer_t hi) const = 0;
382  virtual CSRGraph<integer_t> extract_graph_sep_CB
383  (int ordering_level, integer_t lo, integer_t hi,
384  const std::vector<integer_t>& upd) const = 0;
385  virtual CSRGraph<integer_t> extract_graph_CB_sep
386  (int ordering_level, integer_t lo, integer_t hi,
387  const std::vector<integer_t>& upd) const = 0;
388  virtual CSRGraph<integer_t> extract_graph_CB
389  (int ordering_level, const std::vector<integer_t>& upd) const = 0;
390 
391  virtual void extract_separator
392  (integer_t sep_end, const std::vector<std::size_t>& I,
393  const std::vector<std::size_t>& J, DenseM_t& B, int depth) const = 0;
394  virtual void extract_front
395  (DenseM_t& F11, DenseM_t& F12, DenseM_t& F21,
396  integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
397  int depth) const = 0;
398  virtual void push_front_elements
399  (integer_t, integer_t, const std::vector<integer_t>&,
400  std::vector<Triplet<scalar_t>>&, std::vector<Triplet<scalar_t>>&,
401  std::vector<Triplet<scalar_t>>&) const = 0;
402 
403  virtual void front_multiply
404  (integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
405  const DenseM_t& R, DenseM_t& Sr, DenseM_t& Sc, int depth) const = 0;
406 
407  virtual void front_multiply_F11
408  (Trans op, integer_t slo, integer_t shi,
409  const DenseM_t& R, DenseM_t& S, int depth) const = 0;
410  virtual void front_multiply_F12
411  (Trans op, integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
412  const DenseM_t& R, DenseM_t& S, int depth) const = 0;
413  virtual void front_multiply_F21
414  (Trans op, integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
415  const DenseM_t& R, DenseM_t& S, int depth) const = 0;
416 
417 #if defined(STRUMPACK_USE_MPI)
418  virtual void extract_F11_block
419  (scalar_t* F, integer_t ldF, integer_t row, integer_t nr_rows,
420  integer_t col, integer_t nr_cols) const = 0;
421  virtual void extract_F12_block
422  (scalar_t* F, integer_t ldF, integer_t row,
423  integer_t nr_rows, integer_t col, integer_t nr_cols,
424  const integer_t* upd) const = 0;
425  virtual void extract_F21_block
426  (scalar_t* F, integer_t ldF, integer_t row,
427  integer_t nr_rows, integer_t col, integer_t nr_cols,
428  const integer_t* upd) const = 0;
429  virtual void extract_separator_2d
430  (integer_t sep_end, const std::vector<std::size_t>& I,
431  const std::vector<std::size_t>& J, DistM_t& B) const = 0;
432  virtual void front_multiply_2d
433  (integer_t sep_begin, integer_t sep_end,
434  const std::vector<integer_t>& upd, const DistM_t& R,
435  DistM_t& Srow, DistM_t& Scol, int depth) const = 0;
436  virtual void front_multiply_2d
437  (Trans op, integer_t sep_begin, integer_t sep_end,
438  const std::vector<integer_t>& upd, const DistM_t& R,
439  DistM_t& S, int depth) const = 0;
440 #endif //STRUMPACK_USE_MPI
441 #endif //DOXYGEN_SHOULD_SKIP_THIS
442 
443  protected:
444  integer_t n_, nnz_;
445  std::vector<integer_t> ptr_, ind_;
446  std::vector<scalar_t> val_;
447  bool symm_sparse_;
448 
449  enum MMsym {GENERAL, SYMMETRIC, SKEWSYMMETRIC, HERMITIAN};
450 
451  CompressedSparseMatrix();
452  CompressedSparseMatrix
453  (integer_t n, integer_t nnz, bool symm_sparse=false);
454  CompressedSparseMatrix
455  (integer_t n, const integer_t* row_ptr, const integer_t* col_ind,
456  const scalar_t* values, bool symm_sparsity);
457 
458  std::vector<std::tuple<integer_t,integer_t,scalar_t>>
459  read_matrix_market_entries(const std::string& filename);
460 
461  virtual int strumpack_mc64(MatchingJob, Match_t&) { return 0; }
462 
463  virtual void scale(const std::vector<scalar_t>&,
464  const std::vector<scalar_t>&) = 0;
465  virtual void scale_real(const std::vector<real_t>&,
466  const std::vector<real_t>&) = 0;
467 
468  long long spmv_flops() const;
469  long long spmv_bytes() const;
470  };
471 
472 } //end namespace strumpack
473 
474 #endif // COMPRESSED_SPARSE_MATRIX_HPP
strumpack::Equilibration
Definition: CompressedSparseMatrix.hpp:110
strumpack::MatchingJob::MAX_SMALLEST_DIAGONAL_2
@ MAX_SMALLEST_DIAGONAL_2
strumpack::Trans::N
@ N
strumpack::CompressedSparseMatrix::ind
const integer_t & ind(integer_t i) const
Definition: CompressedSparseMatrix.hpp:261
strumpack::Trans::C
@ C
strumpack::MatchingJob::MAX_DIAGONAL_SUM
@ MAX_DIAGONAL_SUM
DenseMatrix.hpp
Contains the DenseMatrix and DenseMatrixWrapper classes, simple wrappers around BLAS/LAPACK style den...
strumpack::CompressedSparseMatrix::spmv
virtual void spmv(const DenseM_t &x, DenseM_t &y) const =0
strumpack::Side::R
@ R
strumpack::MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING
@ MAX_DIAGONAL_PRODUCT_SCALING
strumpack::CompressedSparseMatrix::val
scalar_t & val(integer_t i)
Definition: CompressedSparseMatrix.hpp:290
strumpack
Definition: StrumpackOptions.hpp:42
strumpack::CompressedSparseMatrix::permute
virtual void permute(const integer_t *iorder, const integer_t *order)
strumpack::DistributedMatrix
Definition: CompressedSparseMatrix.hpp:58
strumpack::MatchingJob::MAX_CARDINALITY
@ MAX_CARDINALITY
strumpack::CompressedSparseMatrix::size
integer_t size() const
Definition: CompressedSparseMatrix.hpp:183
strumpack::CompressedSparseMatrix::set_symm_sparse
void set_symm_sparse(bool symm_sparse=true)
Definition: CompressedSparseMatrix.hpp:305
strumpack::DenseMatrix
This class represents a matrix, stored in column major format, to allow direct use of BLAS/LAPACK rou...
Definition: CompressedSparseMatrix.hpp:57
strumpack::CompressedSparseMatrix::ind
integer_t & ind(integer_t i)
Definition: CompressedSparseMatrix.hpp:283
strumpack::CompressedSparseMatrix::ind
integer_t * ind()
Definition: CompressedSparseMatrix.hpp:238
StrumpackOptions.hpp
Holds options for the sparse solver.
strumpack::CompressedSparseMatrix::val
const scalar_t * val() const
Definition: CompressedSparseMatrix.hpp:219
strumpack::CompressedSparseMatrix::val
const scalar_t & val(integer_t i) const
Definition: CompressedSparseMatrix.hpp:268
strumpack::CompressedSparseMatrix::ind
const integer_t * ind() const
Definition: CompressedSparseMatrix.hpp:210
strumpack::CompressedSparseMatrix::symm_sparse
bool symm_sparse() const
Definition: CompressedSparseMatrix.hpp:296
strumpack::CompressedSparseMatrix::ptr
integer_t & ptr(integer_t i)
Definition: CompressedSparseMatrix.hpp:275
strumpack::MatchingJob::NONE
@ NONE
strumpack::CompressedSparseMatrix::val
scalar_t * val()
Definition: CompressedSparseMatrix.hpp:246
strumpack::CompressedSparseMatrix::ptr
const integer_t * ptr() const
Definition: CompressedSparseMatrix.hpp:201
strumpack::CompressedSparseMatrix::nnz
integer_t nnz() const
Definition: CompressedSparseMatrix.hpp:190
strumpack::CompressedSparseMatrix::ptr
integer_t * ptr()
Definition: CompressedSparseMatrix.hpp:229
strumpack::MatchingJob::MAX_SMALLEST_DIAGONAL
@ MAX_SMALLEST_DIAGONAL
strumpack::CompressedSparseMatrix
Abstract base class for compressed sparse matrix storage.
Definition: CompressedSparseMatrix.hpp:163
strumpack::CompressedSparseMatrix::~CompressedSparseMatrix
virtual ~CompressedSparseMatrix()
Definition: CompressedSparseMatrix.hpp:176
strumpack::Trans
Trans
Definition: DenseMatrix.hpp:50
strumpack::CompressedSparseMatrix::ptr
const integer_t & ptr(integer_t i) const
Definition: CompressedSparseMatrix.hpp:253
strumpack::MatchingData
Definition: CompressedSparseMatrix.hpp:69
strumpack::CSRGraph
Definition: CompressedSparseMatrix.hpp:56
strumpack::MatchingJob
MatchingJob
Definition: StrumpackOptions.hpp:95