All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CSRMatrix.h
Go to the documentation of this file.
1 /*
2  * CSRMatrix.h
3  *
4  * Created on: May 6, 2015
5  * Author: Michael Wegner (michael.wegner@student.kit.edu)
6  */
7 
8 #ifndef CSRMATRIX_H_
9 #define CSRMATRIX_H_
10 
11 #include <vector>
12 #include "../Globals.h"
13 #include "AlgebraicGlobals.h"
14 #include "Vector.h"
15 #include "../graph/Graph.h"
16 #include "../algebraic/SparseAccumulator.h"
17 #include "../auxiliary/Timer.h"
18 
19 namespace NetworKit {
20 
26 class CSRMatrix {
27 private:
28  std::vector<index> rowIdx;
29  std::vector<index> columnIdx;
30  std::vector<double> nonZeros;
31 
32  count nRows;
33  count nCols;
34  bool isSorted;
35  double zero;
36 
42  void quicksort(index left, index right);
43 
50  index partition(index left, index right);
51 
62  index binarySearchColumns(index left, index right, index j) const;
63 
64 public:
66  CSRMatrix();
67 
73  CSRMatrix(const count dimension, const double zero = 0.0);
74 
81  CSRMatrix(const count nRows, const count nCols, const double zero = 0.0);
82 
90  CSRMatrix(const count dimension, const std::vector<Triplet>& triplets, const double zero = 0.0, bool isSorted = false);
91 
100  CSRMatrix(const count nRows, const count nCols, const std::vector<Triplet>& triplets, const double zero = 0.0, bool isSorted = false);
101 
111  CSRMatrix(const count nRows, const count nCols, const std::vector<std::vector<index>> &columnIdx, const std::vector<std::vector<double>> &values, const double zero = 0.0, bool isSorted = false);
112 
123  CSRMatrix(const count nRows, const count nCols, const std::vector<index>& rowIdx, const std::vector<index>& columnIdx, const std::vector<double>& nonZeros, const double zero = 0.0, bool isSorted = false);
124 
126  CSRMatrix (const CSRMatrix &other) = default;
127 
129  CSRMatrix (CSRMatrix &&other) = default;
130 
132  virtual ~CSRMatrix() = default;
133 
135  CSRMatrix& operator=(CSRMatrix &&other) = default;
136 
138  CSRMatrix& operator=(const CSRMatrix &other) = default;
139 
145  bool operator==(const CSRMatrix& other) const {
146  bool equal = nRows == other.nRows && nCols == other.nCols && zero == other.zero;
147  if (equal) {
148  forNonZeroElementsInRowOrder([&](index i, index j, double value) {
149  if (other(i,j) != value) {
150  equal = false;
151  return;
152  }
153  });
154  }
155 
156  return equal;
157  }
158 
164  bool operator!=(const CSRMatrix& other) const {
165  return !((*this) == other);
166  }
167 
171  inline count numberOfRows() const {
172  return nRows;
173  }
174 
178  inline count numberOfColumns() const {
179  return nCols;
180  }
181 
185  inline double getZero() const {
186  return zero;
187  }
188 
193  count nnzInRow(const index i) const;
194 
198  count nnz() const;
199 
203  double operator()(const index i, const index j) const;
204 
209  void setValue(const index i, const index j, const double value);
210 
214  void sort();
215 
219  bool sorted() const;
220 
224  Vector row(const index i) const;
225 
229  Vector column(const index j) const;
230 
234  Vector diagonal() const;
235 
240  CSRMatrix operator+(const CSRMatrix &other) const;
241 
246  CSRMatrix& operator+=(const CSRMatrix &other);
247 
253  CSRMatrix operator-(const CSRMatrix &other) const;
254 
259  CSRMatrix& operator-=(const CSRMatrix &other);
260 
265  CSRMatrix operator*(const double &scalar) const;
266 
271  CSRMatrix& operator*=(const double &scalar);
272 
277  Vector operator*(const Vector &vector) const;
278 
283  CSRMatrix operator*(const CSRMatrix &other) const;
284 
289  CSRMatrix operator/(const double &divisor) const;
290 
295  CSRMatrix& operator/=(const double &divisor);
296 
305  template<typename L> static CSRMatrix binaryOperator(const CSRMatrix &A, const CSRMatrix &B, L binaryOp);
306 
314  static CSRMatrix mTmMultiply(const CSRMatrix &A, const CSRMatrix &B);
315 
323  static CSRMatrix mmTMultiply(const CSRMatrix &A, const CSRMatrix &B);
324 
332  static Vector mTvMultiply(const CSRMatrix &matrix, const Vector &vector);
333 
337  CSRMatrix transpose() const;
338 
346  CSRMatrix extract(const std::vector<index>& rowIndices, const std::vector<index>& columnIndices) const;
347 
357  void assign(const std::vector<index>& rowIndices, const std::vector<index>& columnIndices, const CSRMatrix& source);
358 
364  template<typename F>
365  void apply(const F unaryElementFunction);
366 
371  static CSRMatrix adjacencyMatrix(const Graph& graph, double zero = 0.0);
372 
378  static CSRMatrix diagonalMatrix(const Vector& diagonalElements, double zero = 0.0);
379 
384  static CSRMatrix incidenceMatrix(const Graph& graph, double zero = 0.0);
385 
390  static CSRMatrix laplacianMatrix(const Graph& graph, double zero = 0.0);
391 
396  static CSRMatrix normalizedLaplacianMatrix(const Graph& graph, double zero = 0.0);
397 
398 
399 
403  template<typename L> void forNonZeroElementsInRow(index row, L handle) const;
404 
408  template<typename L> void parallelForNonZeroElementsInRow(index row, L handle) const;
409 
413  template<typename L> void forElementsInRow(index i, L handle) const;
414 
418  template<typename L> void forNonZeroElementsInRowOrder(L handle) const;
419 
423  template<typename L> void parallelForNonZeroElementsInRowOrder(L handle) const;
424 };
425 
426 template<typename L> inline CSRMatrix NetworKit::CSRMatrix::binaryOperator(const CSRMatrix &A, const CSRMatrix &B, L binaryOp) {
427  assert(A.nRows == B.nRows && A.nCols == B.nCols);
428 
429  if (A.sorted() && B.sorted()) {
430  std::vector<index> rowIdx(A.nRows+1);
431  std::vector<std::vector<index>> columns(A.nRows);
432 
433  rowIdx[0] = 0;
434 #pragma omp parallel for
435  for (index i = 0; i < A.nRows; ++i) {
436  index k = A.rowIdx[i];
437  index l = B.rowIdx[i];
438  while (k < A.rowIdx[i+1] && l < B.rowIdx[i+1]) {
439  if (A.columnIdx[k] < B.columnIdx[l]) {
440  columns[i].push_back(A.columnIdx[k]);
441  ++k;
442  } else if (A.columnIdx[k] > B.columnIdx[l]) {
443  columns[i].push_back(B.columnIdx[l]);
444  ++l;
445  } else { // A.columnIdx[k] == B.columnIdx[l]
446  columns[i].push_back(A.columnIdx[k]);
447  ++k;
448  ++l;
449  }
450  ++rowIdx[i+1];
451  }
452 
453  while (k < A.rowIdx[i+1]) {
454  columns[i].push_back(A.columnIdx[k]);
455  ++k;
456  ++rowIdx[i+1];
457  }
458 
459  while (l < B.rowIdx[i+1]) {
460  columns[i].push_back(B.columnIdx[l]);
461  ++l;
462  ++rowIdx[i+1];
463  }
464  }
465 
466 
467  for (index i = 0; i < A.nRows; ++i) {
468  rowIdx[i+1] += rowIdx[i];
469  }
470 
471  count nnz = rowIdx[A.nRows];
472  std::vector<index> columnIdx(nnz);
473  std::vector<double> nonZeros(nnz, A.zero);
474 
475 #pragma omp parallel for
476  for (index i = 0; i < A.nRows; ++i) {
477  for (index cIdx = rowIdx[i], j = 0; cIdx < rowIdx[i+1]; ++cIdx, ++j) {
478  columnIdx[cIdx] = columns[i][j];
479  }
480  columns[i].clear();
481  columns[i].resize(0);
482  columns[i].shrink_to_fit();
483  }
484 
485 #pragma omp parallel for
486  for (index i = 0; i < A.nRows; ++i) {
487  index k = A.rowIdx[i];
488  index l = B.rowIdx[i];
489  for (index cIdx = rowIdx[i]; cIdx < rowIdx[i+1]; ++cIdx) {
490  if (k < A.rowIdx[i+1] && columnIdx[cIdx] == A.columnIdx[k]) {
491  nonZeros[cIdx] = A.nonZeros[k];
492  ++k;
493  }
494 
495  if (l < B.rowIdx[i+1] && columnIdx[cIdx] == B.columnIdx[l]) {
496  nonZeros[cIdx] = binaryOp(nonZeros[cIdx], B.nonZeros[l]);
497  ++l;
498  }
499  }
500  }
501 
502  return CSRMatrix(A.nRows, A.nCols, rowIdx, columnIdx, nonZeros, A.zero, true);
503  } else { // A or B not sorted
504  std::vector<int64_t> columnPointer(A.nCols, -1);
505  std::vector<double> Arow(A.nCols, A.zero);
506  std::vector<double> Brow(A.nCols, B.zero);
507 
508  std::vector<Triplet> triplets;
509 
510  for (index i = 0; i < A.nRows; ++i) {
511  index listHead = 0;
512  count nnz = 0;
513 
514  // search for nonZeros in our own matrix
515  for (index k = A.rowIdx[i]; k < A.rowIdx[i+1]; ++k) {
516  index j = A.columnIdx[k];
517  Arow[j] = A.nonZeros[k];
518 
519  columnPointer[j] = listHead;
520  listHead = j;
521  nnz++;
522  }
523 
524  // search for nonZeros in the other matrix
525  for (index k = B.rowIdx[i]; k < B.rowIdx[i+1]; ++k) {
526  index j = B.columnIdx[k];
527  Brow[j] = B.nonZeros[k];
528 
529  if (columnPointer[j] == -1) { // our own matrix does not have a nonZero entry in column j
530  columnPointer[j] = listHead;
531  listHead = j;
532  nnz++;
533  }
534  }
535 
536  // apply operator on the found nonZeros in A and B
537  for (count k = 0; k < nnz; ++k) {
538  double value = binaryOp(Arow[listHead], Brow[listHead]);
539  if (value != A.zero) {
540  triplets.push_back({i,listHead,value});
541  }
542 
543  index temp = listHead;
544  listHead = columnPointer[listHead];
545 
546  // reset for next row
547  columnPointer[temp] = -1;
548  Arow[temp] = A.zero;
549  Brow[temp] = B.zero;
550  }
551 
552  nnz = 0;
553  }
554 
555  return CSRMatrix(A.numberOfRows(), A.numberOfColumns(), triplets);
556  }
557 }
558 
559 template<typename F>
560 void CSRMatrix::apply(const F unaryElementFunction) {
561 #pragma omp parallel for
562  for (index k = 0; k < nonZeros.size(); ++k) {
563  nonZeros[k] = unaryElementFunction(nonZeros[k]);
564  }
565 }
566 
567 } /* namespace NetworKit */
568 
569 template<typename L>
570 inline void NetworKit::CSRMatrix::forNonZeroElementsInRow(index i, L handle) const {
571  for (index k = rowIdx[i]; k < rowIdx[i+1]; ++k) {
572  handle(columnIdx[k], nonZeros[k]);
573  }
574 }
575 
576 template<typename L>
578 #pragma omp parallel for
579  for (index k = rowIdx[i]; k < rowIdx[i+1]; ++k) {
580  handle(columnIdx[k], nonZeros[k]);
581  }
582 }
583 
584 template<typename L>
585 inline void NetworKit::CSRMatrix::forElementsInRow(index i, L handle) const {
586  Vector rowVector = row(i);
587  index j = 0;
588  rowVector.forElements([&](double val) {
589  handle(j++, val);
590  });
591 }
592 
593 template<typename L>
595  for (index i = 0; i < nRows; ++i) {
596  for (index k = rowIdx[i]; k < rowIdx[i+1]; ++k) {
597  handle(i, columnIdx[k], nonZeros[k]);
598  }
599  }
600 }
601 
602 template<typename L>
604 #pragma omp parallel for
605  for (index i = 0; i < nRows; ++i) {
606  for (index k = rowIdx[i]; k < rowIdx[i+1]; ++k) {
607  handle(i, columnIdx[k], nonZeros[k]);
608  }
609  }
610 }
611 
612 #endif /* TESTMATRIX_H_ */
void parallelForNonZeroElementsInRow(index row, L handle) const
Iterate in parallel over all non-zero elements of row row in the matrix and call handler(index column...
Definition: CSRMatrix.h:577
CSRMatrix operator-(const CSRMatrix &other) const
Subtracts other from this matrix and returns the result.
Definition: CSRMatrix.cpp:309
CSRMatrix operator+(const CSRMatrix &other) const
Adds this matrix to other and returns the result.
Definition: CSRMatrix.cpp:298
bool sorted() const
Definition: CSRMatrix.cpp:236
static Vector mTvMultiply(const CSRMatrix &matrix, const Vector &vector)
Computes matrix^T * vector.
Definition: CSRMatrix.cpp:492
CSRMatrix & operator-=(const CSRMatrix &other)
Subtracts other from this matrix.
Definition: CSRMatrix.cpp:314
void assign(const std::vector< index > &rowIndices, const std::vector< index > &columnIndices, const CSRMatrix &source)
Assign the contents of the matrix source to this matrix at rows and columns specified by rowIndices a...
Definition: CSRMatrix.cpp:561
virtual ~CSRMatrix()=default
Default destructor.
bool operator==(const CSRMatrix &other) const
Compares this matrix to other and returns true if the shape and zero element are the same as well as ...
Definition: CSRMatrix.h:145
count numberOfColumns() const
Definition: CSRMatrix.h:178
uint64_t index
Typedefs.
Definition: Globals.h:20
static CSRMatrix diagonalMatrix(const Vector &diagonalElements, double zero=0.0)
Creates a diagonal matrix with dimension equal to the dimension of the Vector diagonalElements.
Definition: CSRMatrix.cpp:586
CSRMatrix operator/(const double &divisor) const
Divides this matrix by a divisor specified in divisor and returns the result in a new matrix...
Definition: CSRMatrix.cpp:423
count numberOfRows() const
Definition: CSRMatrix.h:171
void parallelForNonZeroElementsInRowOrder(L handle) const
Iterate in parallel over all rows and call handler (lambda closure) on non-zero elements of the matri...
Definition: CSRMatrix.h:603
void apply(const F unaryElementFunction)
Applies the unary function unaryElementFunction to each value in the matrix.
Definition: CSRMatrix.h:560
static CSRMatrix normalizedLaplacianMatrix(const Graph &graph, double zero=0.0)
Returns the (weighted) normalized Laplacian matrix of the (weighted) Graph graph. ...
Definition: CSRMatrix.cpp:652
void forNonZeroElementsInRowOrder(L handle) const
Iterate over all non-zero elements of the matrix in row order and call handler (lambda closure)...
Definition: CSRMatrix.h:594
static CSRMatrix mmTMultiply(const CSRMatrix &A, const CSRMatrix &B)
Computes A * B^T.
Definition: CSRMatrix.cpp:460
CSRMatrix & operator/=(const double &divisor)
Divides this matrix by a divisor specified in divisor.
Definition: CSRMatrix.cpp:427
static CSRMatrix laplacianMatrix(const Graph &graph, double zero=0.0)
Compute the (weighted) Laplacian of the (weighted) Graph graph.
Definition: CSRMatrix.cpp:633
bool operator!=(const CSRMatrix &other) const
Compares this matrix to other and returns false if the shape and zero element are the same as well as...
Definition: CSRMatrix.h:164
static CSRMatrix adjacencyMatrix(const Graph &graph, double zero=0.0)
Compute the (weighted) adjacency matrix of the (weighted) Graph graph.
Definition: CSRMatrix.cpp:572
CSRMatrix & operator=(CSRMatrix &&other)=default
Default move assignment operator.
void setValue(const index i, const index j, const double value)
Set the matrix at position (i, j) to value.
Definition: CSRMatrix.cpp:119
uint64_t count
Definition: Globals.h:21
CSRMatrix & operator+=(const CSRMatrix &other)
Adds other to this matrix.
Definition: CSRMatrix.cpp:303
CSRMatrix()
Default constructor.
Definition: CSRMatrix.cpp:16
The Vector class represents a basic vector with double coefficients.
Definition: Vector.h:25
CSRMatrix transpose() const
Transposes this matrix and returns it.
Definition: CSRMatrix.cpp:505
bool equal(const double x, const double y, const double error)
Test doubles for equality within a given error.
Definition: NumericTools.cpp:16
The CSRMatrix class represents a sparse matrix stored in CSR-Format (i.e.
Definition: CSRMatrix.h:26
CSRMatrix operator*(const double &scalar) const
Multiplies this matrix with a scalar specified in scalar and returns the result.
Definition: CSRMatrix.cpp:320
count nnz() const
Definition: CSRMatrix.cpp:91
Vector diagonal() const
Definition: CSRMatrix.cpp:263
void forElementsInRow(index i, L handle) const
Iterate over all elements in row i in the matrix and call handle(index column, double value) ...
Definition: CSRMatrix.h:585
static CSRMatrix incidenceMatrix(const Graph &graph, double zero=0.0)
Returns the (weighted) incidence matrix of the (weighted) Graph graph.
Definition: CSRMatrix.cpp:603
CSRMatrix & operator*=(const double &scalar)
Multiplies this matrix with a scalar specified in scalar.
Definition: CSRMatrix.cpp:324
void forNonZeroElementsInRow(index row, L handle) const
Iterate over all non-zero elements of row row in the matrix and call handler(index column...
Definition: CSRMatrix.h:570
A graph (with optional weights) and parallel iterator methods.
Definition: Graph.h:79
static CSRMatrix mTmMultiply(const CSRMatrix &A, const CSRMatrix &B)
Computes A^T * B.
Definition: CSRMatrix.cpp:431
double getZero() const
Returns the zero element of the matrix.
Definition: CSRMatrix.h:185
void sort()
Sorts the column indices in each row for faster access.
Definition: CSRMatrix.cpp:210
void forElements(L handle)
Iterate over all elements of the vector and call handler (lambda closure).
Definition: Vector.h:330
Vector row(const index i) const
Definition: CSRMatrix.cpp:240
double operator()(const index i, const index j) const
Definition: CSRMatrix.cpp:95
static CSRMatrix binaryOperator(const CSRMatrix &A, const CSRMatrix &B, L binaryOp)
Computes A binaryOp B on the elements of matrix A and matrix B.
Definition: CSRMatrix.h:426
Vector column(const index j) const
Definition: CSRMatrix.cpp:251
CSRMatrix extract(const std::vector< index > &rowIndices, const std::vector< index > &columnIndices) const
Extracts a matrix with rows and columns specified by rowIndices and columnIndices from this matrix...
Definition: CSRMatrix.cpp:537
count nnzInRow(const index i) const
Definition: CSRMatrix.cpp:86