↰ Return to documentation for file (include/networkit/algebraic/DenseMatrix.hpp)
/*
* DenseMatrix.hpp
*
* Created on: Nov 25, 2015
* Author: Michael Wegner (michael.wegner@student.kit.edu)
*/
#ifndef NETWORKIT_ALGEBRAIC_DENSE_MATRIX_HPP_
#define NETWORKIT_ALGEBRAIC_DENSE_MATRIX_HPP_
#include <cassert>
#include <iostream>
#include <vector>
#include <networkit/Globals.hpp>
#include <networkit/algebraic/AlgebraicGlobals.hpp>
#include <networkit/algebraic/Vector.hpp>
#include <networkit/graph/Graph.hpp>
namespace NetworKit {
class DenseMatrix final {
private:
count nRows;
count nCols;
std::vector<double> entries;
double zero;
public:
DenseMatrix();
DenseMatrix(count dimension, double zero = 0.0);
DenseMatrix(count nRows, count nCols, double zero = 0.0);
DenseMatrix(count dimension, const std::vector<Triplet> &triplets, double zero = 0.0);
DenseMatrix(count nRows, count nCols, const std::vector<Triplet> &triplets, double zero = 0.0);
DenseMatrix(count nRows, count nCols, const std::vector<double> &entries, double zero = 0.0);
~DenseMatrix() = default;
DenseMatrix(const DenseMatrix &other) = default;
DenseMatrix(DenseMatrix &&other) = default;
DenseMatrix &operator=(DenseMatrix &&other) = default;
DenseMatrix &operator=(const DenseMatrix &other) = default;
bool operator==(const DenseMatrix &other) const {
bool equal =
nRows == other.nRows && nCols == other.nCols && entries.size() == other.entries.size();
if (equal) {
forElementsInRowOrder([&](index i, index j, double value) {
if (other(i, j) != value) {
equal = false;
return;
}
});
}
return equal;
}
bool isApprox(const DenseMatrix &other, const double eps = 0.01) const {
bool equal =
nRows == other.nRows && nCols == other.nCols && entries.size() == other.entries.size();
if (equal) {
forElementsInRowOrder([&](index i, index j, double value) {
if (std::abs(other(i, j) - value) > eps) {
equal = false;
return;
}
});
}
return equal;
}
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;
double &operator()(index i, index j);
void setValue(index i, index j, double value);
Vector row(index i) const;
Vector column(index j) const;
Vector diagonal() const;
DenseMatrix operator+(const DenseMatrix &other) const;
DenseMatrix &operator+=(const DenseMatrix &other);
DenseMatrix operator-(const DenseMatrix &other) const;
DenseMatrix &operator-=(const DenseMatrix &other);
DenseMatrix operator*(double scalar) const;
DenseMatrix &operator*=(double scalar);
Vector operator*(const Vector &vector) const;
DenseMatrix operator*(const DenseMatrix &other) const;
DenseMatrix operator/(double divisor) const;
DenseMatrix &operator/=(double divisor);
DenseMatrix transpose() const;
DenseMatrix 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 DenseMatrix &source);
template <typename F>
void apply(F unaryElementFunction);
static DenseMatrix adjacencyMatrix(const Graph &graph, double zero = 0.0);
static DenseMatrix diagonalMatrix(const Vector &diagonalElements, double zero = 0.0);
static DenseMatrix incidenceMatrix(const Graph &graph, double zero = 0.0);
static DenseMatrix laplacianMatrix(const Graph &graph, double zero = 0.0);
static void LUDecomposition(DenseMatrix &matrix);
static Vector LUSolve(const DenseMatrix &LU, const Vector &b);
template <typename L>
static DenseMatrix binaryOperator(const DenseMatrix &A, const DenseMatrix &B, L binaryOp);
template <typename L>
void forElementsInRow(index row, 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 forNonZeroElementsInRow(index row, L handle) const;
template <typename L>
void parallelForNonZeroElementsInRow(index row, L handle) const;
template <typename L>
void forNonZeroElementsInRowOrder(L handle) const;
template <typename L>
void parallelForNonZeroElementsInRowOrder(L handle) const;
};
template <typename F>
void DenseMatrix::apply(F unaryElementFunction) {
#pragma omp parallel for
for (omp_index k = 0; k < static_cast<omp_index>(entries.size()); ++k) {
entries[k] = unaryElementFunction(entries[k]);
}
}
template <typename L>
inline DenseMatrix DenseMatrix::binaryOperator(const DenseMatrix &A, const DenseMatrix &B,
L binaryOp) {
assert(A.nRows == B.nRows && A.nCols == B.nCols);
std::vector<double> resultEntries(A.numberOfRows() * A.numberOfColumns(), 0.0);
#pragma omp parallel for
for (omp_index i = 0; i < static_cast<omp_index>(A.numberOfRows()); ++i) {
index offset = i * A.numberOfColumns();
for (index j = offset; j < offset + A.numberOfColumns(); ++j) {
resultEntries[j] = binaryOp(A.entries[j], B.entries[j]);
}
}
return DenseMatrix(A.numberOfRows(), A.numberOfColumns(), resultEntries);
}
template <typename L>
inline void DenseMatrix::forElementsInRow(index i, L handle) const {
index offset = i * numberOfColumns();
for (index k = offset, j = 0; k < offset + numberOfColumns(); ++k, ++j) {
handle(j, entries[k]);
}
}
template <typename L>
inline void DenseMatrix::parallelForElementsInRow(index i, L handle) const {
index offset = i * numberOfColumns();
#pragma omp parallel for
for (omp_index j = 0; j < static_cast<omp_index>(numberOfColumns()); ++j) {
handle(j, entries[offset + j]);
}
}
template <typename L>
inline void DenseMatrix::forElementsInRowOrder(L handle) const {
for (index i = 0; i < nRows; ++i) {
index offset = i * numberOfColumns();
for (index k = offset, j = 0; k < offset + numberOfColumns(); ++k, ++j) {
handle(i, j, entries[k]);
}
}
}
template <typename L>
inline void DenseMatrix::parallelForElementsInRowOrder(L handle) const {
#pragma omp parallel for
for (omp_index i = 0; i < static_cast<omp_index>(nRows); ++i) {
index offset = i * numberOfColumns();
for (index k = offset, j = 0; k < offset + numberOfColumns(); ++k, ++j) {
handle(i, j, entries[k]);
}
}
}
template <typename L>
inline void DenseMatrix::forNonZeroElementsInRow(index row, L handle) const {
for (index j = 0, k = row * numberOfColumns(); j < numberOfColumns(); ++j, ++k) {
if (entries[k] != getZero()) {
handle(j, entries[k]);
}
}
}
template <typename L>
inline void DenseMatrix::parallelForNonZeroElementsInRow(index row, L handle) const {
#pragma omp parallel for
for (omp_index j = 0; j < static_cast<omp_index>(numberOfColumns()); ++j) {
index k = row * numberOfColumns() + j;
if (entries[k] != getZero()) {
handle(j, entries[k]);
}
}
}
template <typename L>
inline void DenseMatrix::forNonZeroElementsInRowOrder(L handle) const {
for (index i = 0; i < numberOfRows(); ++i) {
for (index j = 0, k = i * numberOfColumns(); j < numberOfColumns(); ++j, ++k) {
if (entries[k] != getZero()) {
handle(i, j, entries[k]);
}
}
}
}
template <typename L>
inline void DenseMatrix::parallelForNonZeroElementsInRowOrder(L handle) const {
#pragma omp parallel for
for (omp_index i = 0; i < static_cast<omp_index>(numberOfRows()); ++i) {
for (index j = 0, k = i * numberOfColumns(); j < numberOfColumns(); ++j, ++k) {
if (entries[k] != getZero()) {
handle(i, j, entries[k]);
}
}
}
}
// print functions for test debugging / output
inline std::ostream &operator<<(std::ostream &os, const DenseMatrix &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_DENSE_MATRIX_HPP_