Program Listing for File ConjugateGradient.hpp¶

/*
 * ConjugateGradient.hpp
 *
 *  Created on: 15.06.2014
 *      Author: Daniel Hoske and Michael Wegner
 */

#ifndef NETWORKIT_NUMERICS_CONJUGATE_GRADIENT_HPP_
#define NETWORKIT_NUMERICS_CONJUGATE_GRADIENT_HPP_

#include <cstdint>
#include <utility>

#include <networkit/algebraic/CSRMatrix.hpp>
#include <networkit/algebraic/Vector.hpp>
#include <networkit/numerics/LinearSolver.hpp>

namespace NetworKit {

template <class Matrix, class Preconditioner>
class ConjugateGradient : public LinearSolver<Matrix> {
public:
    ConjugateGradient(double tolerance = 1e-5)
        : LinearSolver<Matrix>(tolerance), matrix(Matrix()) {}

    void setup(const Matrix &matrix) override {
        this->matrix = matrix;
        precond = Preconditioner(matrix);
    }

    void setupConnected(const Matrix &matrix) override {
        this->matrix = matrix;
        precond = Preconditioner(matrix);
    }

    SolverStatus solve(const Vector &rhs, Vector &result, count maxConvergenceTime = 5 * 60 * 1000,
                       count maxIterations = std::numeric_limits<count>::max()) const override;

    std::vector<SolverStatus>
    parallelSolve(const std::vector<Vector> &rhs, std::vector<Vector> &results,
                  count maxConvergenceTime = 5 * 60 * 1000,
                  count maxIterations = std::numeric_limits<count>::max()) const override;

    template <typename RHSLoader, typename ResultProcessor>
    void parallelSolve(const RHSLoader &rhsLoader, const ResultProcessor &resultProcessor,
                       std::pair<count, count> rhsSize, count maxConvergenceTime = 5 * 60 * 1000,
                       count maxIterations = std::numeric_limits<count>::max()) {
        const index numThreads = omp_get_max_threads();
        count n = rhsSize.first;
        count m = rhsSize.second;
        std::vector<Vector> results(numThreads, Vector(m));
        std::vector<Vector> RHSs(numThreads, Vector(m));

#pragma omp parallel for
        for (omp_index i = 0; i < static_cast<omp_index>(n); ++i) {
            index threadId = omp_get_thread_num();
            const Vector &rhs = rhsLoader(i, RHSs[threadId]);
            Vector &result = results[threadId];

            this->solve(rhs, result, maxConvergenceTime, maxIterations);
            resultProcessor(i, result);
        }
    }

private:
    Matrix matrix;
    Preconditioner precond;
};

template <class Matrix, class Preconditioner>
SolverStatus ConjugateGradient<Matrix, Preconditioner>::solve(const Vector &rhs, Vector &result,
                                                              count, count maxIterations) const {
    assert(matrix.numberOfRows() == rhs.getDimension());

    // Absolute residual to achieve
    double sqr_desired_residual = this->tolerance * this->tolerance * (rhs.length() * rhs.length());

    // Main loop. See:
    // http://en.wikipedia.org/wiki/Conjugate_gradient_method#The_resulting_algorithm
    Vector residual_dir = rhs - matrix * result;
    Vector conjugate_dir = precond.rhs(residual_dir);
    double sqr_residual = Vector::innerProduct(residual_dir, residual_dir);
    double sqr_residual_precond = Vector::innerProduct(residual_dir, conjugate_dir);

    count niters = 0;
    Vector tmp, residual_precond;
    while (sqr_residual > sqr_desired_residual) {
        niters++;
        if (niters > maxIterations) {
            break;
        }

        tmp = matrix * conjugate_dir;
        double step = sqr_residual_precond / Vector::innerProduct(conjugate_dir, tmp);
        result += step * conjugate_dir;
        residual_dir -= step * tmp;
        sqr_residual = Vector::innerProduct(residual_dir, residual_dir);

        residual_precond = precond.rhs(residual_dir);
        double new_sqr_residual_precond = Vector::innerProduct(residual_dir, residual_precond);
        conjugate_dir =
            (new_sqr_residual_precond / sqr_residual_precond) * conjugate_dir + residual_precond;
        sqr_residual_precond = new_sqr_residual_precond;
    }

    SolverStatus status;
    status.numIters = niters;
    status.residual = (rhs - matrix * result).length();
    status.converged = status.residual / rhs.length() <= this->tolerance;

    return status;
}

template <class Matrix, class Preconditioner>
std::vector<SolverStatus> ConjugateGradient<Matrix, Preconditioner>::parallelSolve(
    const std::vector<Vector> &rhs, std::vector<Vector> &results, count maxConvergenceTime,
    count maxIterations) const {
    std::vector<SolverStatus> stati(rhs.size());
#pragma omp parallel for
    for (omp_index i = 0; i < static_cast<omp_index>(rhs.size()); ++i) {
        stati[i] = this->solve(rhs[i], results[i], maxConvergenceTime, maxIterations);
    }
    return stati;
}

} /* namespace NetworKit */

#endif // NETWORKIT_NUMERICS_CONJUGATE_GRADIENT_HPP_