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 
292  virtual real_t norm1() const = 0; //{ assert(false); return -1.; };
293 
298  bool symm_sparse() const { return symm_sparse_; }
299 
307  void set_symm_sparse(bool symm_sparse=true) { symm_sparse_ = symm_sparse; }
308 
309 
323  virtual void spmv(const DenseM_t& x, DenseM_t& y) const = 0;
324 
338  virtual void spmv(const scalar_t* x, scalar_t* y) const = 0;
339 
344  virtual void permute(const integer_t* iorder, const integer_t* order);
345 
346  virtual void permute(const std::vector<integer_t>& iorder,
347  const std::vector<integer_t>& order) {
348  permute(iorder.data(), order.data());
349  }
350 
351  virtual void permute_columns(const std::vector<integer_t>& perm) = 0;
352 
353  virtual Equil_t equilibration() const { return Equil_t(this->size()); }
354 
355  virtual void equilibrate(const Equil_t&) {}
356 
357  virtual Match_t matching(MatchingJob, bool apply=true);
358 
359  virtual void apply_matching(const Match_t&);
360 
361  virtual void symmetrize_sparsity();
362 
363  virtual void print() const;
364  virtual void print_dense(const std::string& name) const {
365  std::cerr << "print_dense not implemented for this matrix type"
366  << std::endl;
367  }
368  virtual void print_matrix_market(const std::string& filename) const {
369  std::cerr << "print_matrix_market not implemented for this matrix type"
370  << std::endl;
371  }
372 
373  virtual int read_matrix_market(const std::string& filename) = 0;
374 
375  virtual real_t max_scaled_residual(const scalar_t* x,
376  const scalar_t* b) const = 0;
377  virtual real_t max_scaled_residual(const DenseM_t& x,
378  const DenseM_t& b) const = 0;
379 
380 
381 #ifndef DOXYGEN_SHOULD_SKIP_THIS
382  virtual CSRGraph<integer_t> extract_graph
383  (int ordering_level, integer_t lo, integer_t hi) const = 0;
384  virtual CSRGraph<integer_t> extract_graph_sep_CB
385  (int ordering_level, integer_t lo, integer_t hi,
386  const std::vector<integer_t>& upd) const = 0;
387  virtual CSRGraph<integer_t> extract_graph_CB_sep
388  (int ordering_level, integer_t lo, integer_t hi,
389  const std::vector<integer_t>& upd) const = 0;
390  virtual CSRGraph<integer_t> extract_graph_CB
391  (int ordering_level, const std::vector<integer_t>& upd) const = 0;
392 
393  virtual void extract_separator
394  (integer_t sep_end, const std::vector<std::size_t>& I,
395  const std::vector<std::size_t>& J, DenseM_t& B, int depth) const = 0;
396  virtual void extract_front
397  (DenseM_t& F11, DenseM_t& F12, DenseM_t& F21,
398  integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
399  int depth) const = 0;
400  virtual void push_front_elements
401  (integer_t, integer_t, const std::vector<integer_t>&,
402  std::vector<Triplet<scalar_t>>&, std::vector<Triplet<scalar_t>>&,
403  std::vector<Triplet<scalar_t>>&) const = 0;
404 
405  virtual void front_multiply
406  (integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
407  const DenseM_t& R, DenseM_t& Sr, DenseM_t& Sc, int depth) const = 0;
408 
409  virtual void front_multiply_F11
410  (Trans op, integer_t slo, integer_t shi,
411  const DenseM_t& R, DenseM_t& S, int depth) const = 0;
412  virtual void front_multiply_F12
413  (Trans op, integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
414  const DenseM_t& R, DenseM_t& S, int depth) const = 0;
415  virtual void front_multiply_F21
416  (Trans op, integer_t slo, integer_t shi, const std::vector<integer_t>& upd,
417  const DenseM_t& R, DenseM_t& S, int depth) const = 0;
418 
419 #if defined(STRUMPACK_USE_MPI)
420  virtual void extract_F11_block
421  (scalar_t* F, integer_t ldF, integer_t row, integer_t nr_rows,
422  integer_t col, integer_t nr_cols) const = 0;
423  virtual void extract_F12_block
424  (scalar_t* F, integer_t ldF, integer_t row,
425  integer_t nr_rows, integer_t col, integer_t nr_cols,
426  const integer_t* upd) const = 0;
427  virtual void extract_F21_block
428  (scalar_t* F, integer_t ldF, integer_t row,
429  integer_t nr_rows, integer_t col, integer_t nr_cols,
430  const integer_t* upd) const = 0;
431  virtual void extract_separator_2d
432  (integer_t sep_end, const std::vector<std::size_t>& I,
433  const std::vector<std::size_t>& J, DistM_t& B) const = 0;
434  virtual void front_multiply_2d
435  (integer_t sep_begin, integer_t sep_end,
436  const std::vector<integer_t>& upd, const DistM_t& R,
437  DistM_t& Srow, DistM_t& Scol, int depth) const = 0;
438  virtual void front_multiply_2d
439  (Trans op, integer_t sep_begin, integer_t sep_end,
440  const std::vector<integer_t>& upd, const DistM_t& R,
441  DistM_t& S, int depth) const = 0;
442 #endif //STRUMPACK_USE_MPI
443 #endif //DOXYGEN_SHOULD_SKIP_THIS
444 
445  protected:
446  integer_t n_, nnz_;
447  std::vector<integer_t> ptr_, ind_;
448  std::vector<scalar_t> val_;
449  bool symm_sparse_;
450 
451  enum MMsym {GENERAL, SYMMETRIC, SKEWSYMMETRIC, HERMITIAN};
452 
453  CompressedSparseMatrix();
454  CompressedSparseMatrix
455  (integer_t n, integer_t nnz, bool symm_sparse=false);
456  CompressedSparseMatrix
457  (integer_t n, const integer_t* row_ptr, const integer_t* col_ind,
458  const scalar_t* values, bool symm_sparsity);
459 
460  std::vector<std::tuple<integer_t,integer_t,scalar_t>>
461  read_matrix_market_entries(const std::string& filename);
462 
463  virtual int strumpack_mc64(MatchingJob, Match_t&) { return 0; }
464 
465  virtual void scale(const std::vector<scalar_t>&,
466  const std::vector<scalar_t>&) = 0;
467  virtual void scale_real(const std::vector<real_t>&,
468  const std::vector<real_t>&) = 0;
469 
470  long long spmv_flops() const;
471  long long spmv_bytes() const;
472  };
473 
474 } //end namespace strumpack
475 
476 #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:307
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:298
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:98