↰ Return to documentation for file (include/networkit/algebraic/DynamicMatrix.hpp)
/*
* DynamicMatrix.hpp
*
* Created on: 13.03.2014
* Author: Michael Wegner (michael.wegner@student.kit.edu)
*/
#ifndef NETWORKIT_ALGEBRAIC_DYNAMIC_MATRIX_HPP_
#define NETWORKIT_ALGEBRAIC_DYNAMIC_MATRIX_HPP_
#include <iostream>
#include <networkit/algebraic/AlgebraicGlobals.hpp>
#include <networkit/algebraic/SparseAccumulator.hpp>
#include <networkit/algebraic/Vector.hpp>
#include <networkit/graph/Graph.hpp>
namespace NetworKit {
class DynamicMatrix final {
protected:
Graph graph;
count nRows;
count nCols;
double zero;
// necessary for operator() const to allow for index addressing via DynamicMatrix(i,j) = value
class IndexProxy {
public:
IndexProxy(DynamicMatrix &mat, index i, index j) : Matrix{mat}, i{i}, j{j} {}
operator double() const { return const_cast<const DynamicMatrix &>(Matrix)(i, j); }
void operator=(double rhs) { Matrix.setValue(i, j, rhs); }
private:
DynamicMatrix &Matrix;
index i;
index j;
};
public:
DynamicMatrix();
DynamicMatrix(count dimension, double zero = 0.0);
DynamicMatrix(count nRows, count nCols, double zero = 0.0);
DynamicMatrix(count dimension, const std::vector<Triplet> &triplets, double zero = 0.0);
DynamicMatrix(count nRows, count nCols, const std::vector<Triplet> &triplets,
double zero = 0.0);
bool operator==(const DynamicMatrix &other) const {
bool graphsEqual = graph.numberOfNodes() == other.graph.numberOfNodes()
&& graph.numberOfEdges() == other.graph.numberOfEdges();
if (graphsEqual) {
graph.forEdges([&](node u, node v, edgeweight w) {
if (w != other.graph.weight(u, v)) {
graphsEqual = false;
return;
}
});
}
return graphsEqual && nRows == other.nRows && nCols == other.nCols && zero == other.zero;
}
bool isApprox(const DynamicMatrix &other, const double eps = 0.01) const {
bool graphsEqual = graph.numberOfNodes() == other.graph.numberOfNodes()
&& graph.numberOfEdges() == other.graph.numberOfEdges();
if (graphsEqual) {
graph.forEdges([&](node u, node v, edgeweight w) {
if (std::abs(w - other.graph.weight(u, v)) > eps) {
graphsEqual = false;
return;
}
});
}
return graphsEqual && nRows == other.nRows && nCols == other.nCols && zero == other.zero;
}
bool operator!=(const DynamicMatrix &other) const { return !((*this) == other); }
inline count numberOfRows() const { return nRows; }
inline count numberOfColumns() const { return nCols; }
inline double getZero() const { return zero; }
count nnzInRow(index i) const;
count nnz() const;
double operator()(index i, index j) const;
IndexProxy operator()(index i, index j) { return IndexProxy(*this, i, j); }
void setValue(index i, index j, double value);
Vector row(index i) const;
Vector column(index j) const;
Vector diagonal() const;
DynamicMatrix operator+(const DynamicMatrix &other) const;
DynamicMatrix &operator+=(const DynamicMatrix &other);
DynamicMatrix operator-(const DynamicMatrix &other) const;
DynamicMatrix &operator-=(const DynamicMatrix &other);
DynamicMatrix operator*(double scalar) const;
DynamicMatrix &operator*=(double scalar);
Vector operator*(const Vector &vector) const;
DynamicMatrix operator*(const DynamicMatrix &other) const;
DynamicMatrix operator/(double divisor) const;
DynamicMatrix &operator/=(double divisor);
static DynamicMatrix mTmMultiply(const DynamicMatrix &A, const DynamicMatrix &B);
static DynamicMatrix mmTMultiply(const DynamicMatrix &A, const DynamicMatrix &B);
static Vector mTvMultiply(const DynamicMatrix &matrix, const Vector &vector);
DynamicMatrix transpose() const;
DynamicMatrix extract(const std::vector<index> &rowIndices,
const std::vector<index> &columnIndices) const;
void assign(const std::vector<index> &rowIndices, const std::vector<index> &columnIndices,
const DynamicMatrix &source);
template <typename F>
void apply(F unaryElementFunction);
static DynamicMatrix adjacencyMatrix(const Graph &graph, double zero = 0.0);
static DynamicMatrix diagonalMatrix(const Vector &diagonalElements, double zero = 0.0);
static DynamicMatrix incidenceMatrix(const Graph &graph, double zero = 0.0);
static DynamicMatrix laplacianMatrix(const Graph &graph, double zero = 0.0);
static DynamicMatrix normalizedLaplacianMatrix(const Graph &graph, double zero = 0.0);
template <typename L>
void forNonZeroElementsInRow(index row, L handle) const;
template <typename L>
void parallelForNonZeroElementsInRow(index i, L handle) const;
template <typename L>
void forElementsInRow(index i, L handle) const;
template <typename L>
void parallelForElementsInRow(index row, L handle) const;
template <typename L>
void forElementsInRowOrder(L handle) const;
template <typename L>
void parallelForElementsInRowOrder(L handle) const;
template <typename L>
void forNonZeroElementsInRowOrder(L handle) const;
template <typename L>
void parallelForNonZeroElementsInRowOrder(L handle) const;
};
template <typename F>
void NetworKit::DynamicMatrix::apply(F unaryElementFunction) {
forNonZeroElementsInRowOrder(
[&](index i, index j, double value) { setValue(i, j, unaryElementFunction(value)); });
}
template <typename L>
inline void NetworKit::DynamicMatrix::forNonZeroElementsInRow(index row, L handle) const {
graph.forEdgesOf(row, [&](index j, edgeweight weight) { handle(j, weight); });
}
template <typename L>
inline void NetworKit::DynamicMatrix::parallelForNonZeroElementsInRow(index i, L handle) const {
std::vector<std::pair<index, edgeweight>> elements(graph.weightNeighborRange(i).begin(),
graph.weightNeighborRange(i).end());
#pragma omp parallel for
for (omp_index j = 0; j < static_cast<omp_index>(elements.size()); ++j) {
handle(elements[j].first, elements[j].second);
}
}
template <typename L>
inline void NetworKit::DynamicMatrix::forElementsInRow(index i, L handle) const {
Vector rowVector = row(i);
index j = 0;
rowVector.forElements([&](double value) { handle(j++, value); });
}
template <typename L>
inline void NetworKit::DynamicMatrix::parallelForElementsInRow(index i, L handle) const {
Vector rowVector = row(i);
#pragma omp parallel for
for (omp_index j = 0; j < static_cast<omp_index>(rowVector.getDimension()); ++j) {
handle(j, rowVector[j]);
}
}
template <typename L>
inline void NetworKit::DynamicMatrix::forElementsInRowOrder(L handle) const {
for (index i = 0; i < nRows; ++i) {
index col = 0;
graph.forEdgesOf(i, [&](index j, edgeweight weight) {
while (col < j) {
handle(i, col, getZero());
++col;
}
handle(i, col, weight);
++col;
});
while (col < i) {
handle(i, col, getZero());
++col;
}
}
}
template <typename L>
inline void NetworKit::DynamicMatrix::parallelForElementsInRowOrder(L handle) const {
#pragma omp parallel for
for (omp_index i = 0; i < static_cast<omp_index>(nRows); ++i) {
index col = 0;
graph.forEdgesOf(i, [&](index j, edgeweight weight) {
while (col < j) {
handle(i, col, getZero());
++col;
}
handle(i, col, weight);
++col;
});
while (col < i) {
handle(i, col, getZero());
++col;
}
}
}
template <typename L>
inline void NetworKit::DynamicMatrix::forNonZeroElementsInRowOrder(L handle) const {
for (index i = 0; i < nRows; ++i) {
graph.forEdgesOf(i, [&](index j, edgeweight weight) { handle(i, j, weight); });
}
}
template <typename L>
inline void NetworKit::DynamicMatrix::parallelForNonZeroElementsInRowOrder(L handle) const {
#pragma omp parallel for
for (omp_index i = 0; i < static_cast<omp_index>(nRows); ++i) {
graph.forEdgesOf(i, [&](index j, edgeweight weight) { handle(i, j, weight); });
}
}
// print functions for test debugging / output
inline std::ostream &operator<<(std::ostream &os, const DynamicMatrix &M) {
for (index row = 0; row < M.numberOfRows(); ++row) {
if (row != 0)
os << std::endl;
for (index col = 0; col < M.numberOfColumns(); ++col) {
if (col != 0)
os << ", ";
os << M(row, col);
}
}
return os;
}
} // namespace NetworKit
#endif // NETWORKIT_ALGEBRAIC_DYNAMIC_MATRIX_HPP_