Program Listing for File 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;

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 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;

    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;
};

} /* namespace NetworKit */

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); });
    }
}

namespace NetworKit {
// 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_