All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GraphBLAS.h
Go to the documentation of this file.
1 /*
2  * GraphBLAS.h
3  *
4  * Created on: May 31, 2016
5  * Author: Michael Wegner (michael.wegner@student.kit.edu)
6  */
7 
8 #ifndef NETWORKIT_CPP_ALGEBRAIC_GRAPHBLAS_H_
9 #define NETWORKIT_CPP_ALGEBRAIC_GRAPHBLAS_H_
10 
11 #include <limits>
12 #include "Semirings.h"
13 #include "SparseAccumulator.h"
14 #include "AlgebraicGlobals.h"
15 #include "Vector.h"
16 
21 namespace GraphBLAS {
22 
23 // ****************************************************
24 // Operations
25 // ****************************************************
26 
35 template<class SemiRing, class Matrix, typename L>
36 Matrix eWiseBinOp(const Matrix& A, const Matrix& B, L binOp) {
37  assert(A.numberOfRows() == B.numberOfRows() && A.numberOfColumns() == B.numberOfColumns());
38  assert(A.getZero() == B.getZero() && A.getZero() == SemiRing::zero());
39 
40  std::vector<int64_t> columnPointer(A.numberOfColumns(), -1);
41  std::vector<double> Arow(A.numberOfColumns(), SemiRing::zero());
42  std::vector<double> Brow(A.numberOfColumns(), SemiRing::zero());
43 
44  std::vector<NetworKit::Triplet> triplets;
45 
46  for (NetworKit::index i = 0; i < A.numberOfRows(); ++i) {
47  NetworKit::index listHead = 0;
48  NetworKit::count nnz = 0;
49 
50  // search for nonZeros in matrix A
51  A.forNonZeroElementsInRow(i, [&](NetworKit::index j, double value) {
52  Arow[j] = value;
53 
54  columnPointer[j] = listHead;
55  listHead = j;
56  nnz++;
57  });
58 
59  // search for nonZeros in matrix B
60  B.forNonZeroElementsInRow(i, [&](NetworKit::index j, double value) {
61  Brow[j] = value;
62 
63  if (columnPointer[j] == -1) { // matrix A does not have a nonZero entry in column j
64  columnPointer[j] = listHead;
65  listHead = j;
66  nnz++;
67  }
68  });
69 
70  // apply operator on the found nonZeros in A and B
71  for (NetworKit::count k = 0; k < nnz; ++k) {
72  double value = binOp(Arow[listHead], Brow[listHead]);
73  if (value != SemiRing::zero()) {
74  triplets.push_back({i,listHead,value});
75  }
76 
77  NetworKit::index temp = listHead;
78  listHead = columnPointer[listHead];
79 
80  // reset for next row
81  columnPointer[temp] = -1;
82  Arow[temp] = SemiRing::zero();
83  Brow[temp] = SemiRing::zero();
84  }
85 
86  nnz = 0;
87  }
88 
89  return Matrix(A.numberOfRows(), A.numberOfColumns(), triplets, A.getZero());
90 }
91 
100 template<class SemiRing = ArithmeticSemiring, class Matrix>
101 Matrix MxM(const Matrix& A, const Matrix& B) {
102  assert(A.numberOfColumns() == B.numberOfRows());
103  assert(A.getZero() == SemiRing::zero() && B.getZero() == SemiRing::zero());
104 
105  std::vector<NetworKit::Triplet> triplets;
106  NetworKit::SparseAccumulator spa(B.numberOfRows());
107  for (NetworKit::index i = 0; i < A.numberOfRows(); ++i) {
108  A.forNonZeroElementsInRow(i, [&](NetworKit::index k, double w1) {
109  B.forNonZeroElementsInRow(k, [&](NetworKit::index j, double w2) {
110  spa.scatter(SemiRing::mult(w1,w2), j, *SemiRing::add);
111  });
112  });
113 
114  spa.gather([&](NetworKit::index i, NetworKit::index j, double value){
115  triplets.push_back({i,j,value});
116  });
117 
118  spa.increaseRow();
119  }
120 
121  return Matrix(A.numberOfRows(), B.numberOfColumns(), triplets, A.getZero());
122 }
123 
132 template<class SemiRing = ArithmeticSemiring, class Matrix>
133 void MxM(const Matrix& A, const Matrix& B, Matrix& C) {
134  assert(A.numberOfColumns() == B.numberOfRows() && A.numberOfRows() == C.numberOfRows() && B.numberOfColumns() == C.numberOfColumns());
135  assert(A.getZero() == SemiRing::zero() && B.getZero() == SemiRing::zero() && C.getZero() == SemiRing::zero());
136 
137  std::vector<NetworKit::Triplet> triplets;
138  NetworKit::SparseAccumulator spa(B.numberOfRows());
139  for (NetworKit::index i = 0; i < A.numberOfRows(); ++i) {
140  A.forNonZeroElementsInRow(i, [&](NetworKit::index k, double w1) {
141  B.forNonZeroElementsInRow(k, [&](NetworKit::index j, double w2) {
142  spa.scatter(SemiRing::mult(w1,w2), j, *SemiRing::add);
143  });
144  });
145 
146  spa.gather([&](NetworKit::index i, NetworKit::index j, double value){
147  triplets.push_back({i,j,value});
148  });
149 
150  spa.increaseRow();
151  }
152 
153  Matrix temp(A.numberOfRows(), B.numberOfRows(), triplets, A.getZero());
154  C = eWiseBinOp<SemiRing, Matrix>(C, temp, *SemiRing::add);
155 }
156 
166 template<class SemiRing = ArithmeticSemiring, typename F, class Matrix>
167 void MxM(const Matrix& A, const Matrix& B, Matrix& C, F accum) {
168  assert(A.numberOfColumns() == B.numberOfRows() && A.numberOfRows() == C.numberOfRows() && B.numberOfColumns() == C.numberOfColumns());
169  assert(A.getZero() == SemiRing::zero() && B.getZero() == SemiRing::zero() && C.getZero() == SemiRing::zero());
170 
171  std::vector<NetworKit::Triplet> triplets;
172  NetworKit::SparseAccumulator spa(B.numberOfRows());
173  for (NetworKit::index i = 0; i < A.numberOfRows(); ++i) {
174  A.forNonZeroElementsInRow(i, [&](NetworKit::index k, double w1) {
175  B.forNonZeroElementsInRow(k, [&](NetworKit::index j, double w2) {
176  spa.scatter(SemiRing::mult(w1,w2), j, *SemiRing::add);
177  });
178  });
179 
180  spa.gather([&](NetworKit::index i, NetworKit::index j, double value){
181  triplets.push_back({i,j,value});
182  });
183 
184  spa.increaseRow();
185  }
186 
187  Matrix temp(A.numberOfRows(), B.numberOfRows(), triplets, A.getZero());
188  C = eWiseBinOp<SemiRing, Matrix>(C, temp, accum);
189 }
190 
196 template<class SemiRing = ArithmeticSemiring, class Matrix>
198  assert(!v.isTransposed());
199  assert(A.numberOfColumns() == v.getDimension());
200  assert(A.getZero() == SemiRing::zero());
201  NetworKit::Vector result(A.numberOfRows(), A.getZero());
202 
203  A.parallelForNonZeroElementsInRowOrder([&](NetworKit::index i, NetworKit::index j, double value) {
204  result[i] = SemiRing::add(result[i], SemiRing::mult(value, v[j]));
205  });
206 
207  return result;
208 }
209 
218 template<class SemiRing = ArithmeticSemiring, class Matrix>
219 void MxV(const Matrix& A, const NetworKit::Vector& v, NetworKit::Vector& c) {
220  assert(!v.isTransposed());
221  assert(A.numberOfColumns() == v.getDimension());
222  assert(A.getZero() == SemiRing::zero());
223 
224  A.parallelForNonZeroElementsInRowOrder([&](NetworKit::index i, NetworKit::index j, double value) {
225  c[i] = SemiRing::add(c[i], SemiRing::mult(value, v[j]));
226  });
227 }
228 
237 template<class SemiRing = ArithmeticSemiring, typename F, class Matrix>
238 void MxV(const Matrix& A, const NetworKit::Vector& v, NetworKit::Vector& c, F accum) {
239  assert(!v.isTransposed());
240  assert(A.numberOfColumns() == v.getDimension());
241  assert(A.getZero() == SemiRing::zero());
242 
243  A.parallelForNonZeroElementsInRowOrder([&](NetworKit::index i, NetworKit::index j, double value) {
244  c[i] = accum(c[i], SemiRing::mult(value, v[j]));
245  });
246 }
247 
254 template<class SemiRing = ArithmeticSemiring, class Matrix>
255 Matrix eWiseAdd(const Matrix& A, const Matrix& B) {
256  return eWiseBinOp<SemiRing, Matrix>(A, B, [](const double a, const double b) {return SemiRing::add(a,b);});
257 }
258 
266 template<class SemiRing = ArithmeticSemiring, class Matrix>
267 Matrix eWiseMult(const Matrix& A, const Matrix& B) {
268  return eWiseBinOp<SemiRing, Matrix>(A, B, [](const double a, const double b) {return SemiRing::mult(a,b);});
269 }
270 
277 template<class SemiRing = ArithmeticSemiring, class Matrix>
279  assert(matrix.getZero() == SemiRing::zero());
280  NetworKit::Vector rowReduction(matrix.numberOfRows(), 0.0);
281 
282 #pragma omp parallel for
283  for (NetworKit::index i = 0; i < matrix.numberOfRows(); ++i) {
284  matrix.forNonZeroElementsInRow(i, [&](NetworKit::index j, double value) {
285  rowReduction[i] = SemiRing::add(rowReduction[i], value);
286  });
287  }
288 
289  return rowReduction;
290 }
291 
298 template<class SemiRing = ArithmeticSemiring, class Matrix>
300  assert(matrix.getZero() == SemiRing::zero());
301  NetworKit::Vector columnReduction(matrix.numberOfColumns(), 0.0);
302 
303  matrix.forNonZeroElementsInRowOrder([&](NetworKit::index i, NetworKit::index j, double value) {
304  columnReduction[j] = SemiRing::add(columnReduction[j], value);
305  });
306 
307  return columnReduction;
308 }
309 
310 }
311 
312 
313 
314 #endif /* NETWORKIT_CPP_ALGEBRAIC_GRAPHBLAS_H_ */
NetworKit::Vector rowReduce(const Matrix &matrix)
Computes the row-reduction of the matrix and returns the result as a vector.
Definition: GraphBLAS.h:278
Matrix eWiseMult(const Matrix &A, const Matrix &B)
Computes SemiRing::mult(A(i,j), B(i,j)) for all i,j element-wise and returns the resulting matrix...
Definition: GraphBLAS.h:267
count getDimension() const
Definition: Vector.h:74
Matrix MxM(const Matrix &A, const Matrix &B)
Computes the matrix-matrix multiplication of A and B.
Definition: GraphBLAS.h:101
uint64_t index
Typedefs.
Definition: Globals.h:20
The SparseAccumulator class represents the sparse accumulator datastructure as described in Kepner...
Definition: SparseAccumulator.h:22
Matrix eWiseBinOp(const Matrix &A, const Matrix &B, L binOp)
Computes binOp(A(i,j), B(i,j)) for all i,j element-wise.
Definition: GraphBLAS.h:36
count gather(L handle)
Calls handle for each non zero value of the current row.
Definition: SparseAccumulator.h:85
NetworKit::Vector MxV(const Matrix &A, const NetworKit::Vector &v)
Computes the matrix-vector product of matrix A and Vector v.
Definition: GraphBLAS.h:197
NetworKit::Vector columnReduce(const Matrix &matrix)
Computes the column-reduction of the matrix and returns the result as a Vector.
Definition: GraphBLAS.h:299
uint64_t count
Definition: Globals.h:21
std::vector< std::vector< count > > Matrix
Definition: DynamicNMIDistance.h:16
Matrix eWiseAdd(const Matrix &A, const Matrix &B)
Computes SemiRing::add(A(i,j), B(i,j)) for all i,j element-wise and returns the resulting matrix...
Definition: GraphBLAS.h:255
The Vector class represents a basic vector with double coefficients.
Definition: Vector.h:25
bool isTransposed() const
A transposed vector is a row vector.
Definition: Vector.cpp:24