sparseMat.h

Go to the documentation of this file.

#ifndef SPARSEMAT_H
#define SPARSEMAT_H
 
#include <set>
#include <map>
#include <vector>
#include <iostream>
#include <cassert>
#include <algorithm>
#include <execution>
#include <mutex>
 
#include "config.h"
#include "algebraCore.h"
 
namespace algebra
{
 
using MatrixShape = std::vector<std::set<int>>;
 
class SparseMatrix
{
    struct SparseVector
    {
        std::mutex mutex;   
        std::vector<int> indices;  
        std::vector<double> values;  
        double& operator[](int j)
            {
            auto it = std::lower_bound(indices.begin(), indices.end(), j);
            assert(it != indices.end() && *it == j);
            int k = it - indices.begin();
            return values[k];
            }
    };
 
public:
    SparseMatrix(const MatrixShape &shape) : rows(shape.size())
        {
        for (size_t i = 0; i < shape.size(); ++i)
            {
            const std::set<int>& row_shape = shape[i];
            SparseVector& row = rows[i];
            row.indices.reserve(row_shape.size());
            for (auto it = row_shape.begin(); it != row_shape.end(); ++it)
                row.indices.push_back(*it);
            row.values.resize(row_shape.size());
            }
        }
 
    void clear()
        {
        for (SparseVector& row: rows)
            for (double& value: row.values)
                value = 0;
        }
 
    void set(int i, int j, double val)
        {
        assert(i >= 0 && i < rows.size());
        assert(j >= 0 && j < rows.size());
        rows[i][j] = val;
        }
 
    void add(int i, int j, double val)
        {
        assert(i >= 0 && i < rows.size());
        assert(j >= 0 && j < rows.size());
        SparseVector& row = rows[i];
        std::mutex& mutex = row.mutex;
        double& value = row[j];
        mutex.lock();
        value += val;
        mutex.unlock();
        }
 
    void print(std::ostream & flux = std::cout) const
        {
        flux << "[\n";
        for (const SparseVector& row: rows)
            {
            flux << "  {";
            for (size_t k = 0; k < row.indices.size(); ++k)
                {
                flux << row.indices[k] << ": " << row.values[k];
                flux << (k < row.indices.size()-1 ? ", " : "\n");
                }
            flux << "}\n";
            }
        flux << ']';
        }
 
    double operator() (const int i, const int j) const
        {
        assert(i >= 0 && i < rows.size());
        assert(j >= 0 && j < rows.size());
        const SparseVector& row = rows[i];
        auto it = std::lower_bound(row.indices.begin(), row.indices.end(), j);
        if (it == row.indices.end() || *it != j)
            return 0;
        int k = it - row.indices.begin();
        return row.values[k];
        }
 
    template <typename T>
    void mult(std::vector<T> const& X, std::vector<T> &Y)
        {
        assert(X.size() == rows.size());
        assert(Y.size() == rows.size());
        std::transform(EXEC_POL, rows.begin(), rows.end(), Y.begin(),
            [&X](const SparseVector &row)
            {
            double val = 0;
            for (size_t k = 0; k < row.indices.size(); ++k)
                val += row.values[k] * X[row.indices[k]];
            return val;
            });
        }
 
    template <typename T>
    void build_diag_precond(std::vector<T> &D) const
        {
        assert(D.size() == rows.size());
        for (size_t i = 0; i < rows.size(); ++i)
            {
            const double c = (*this)(i,i);
            assert(c != 0);
            D[i] = 1.0 / c;
            }
        }
 
private:
    std::vector<SparseVector> rows;
}; // end class SparseMatrix
 
} // end namespace algebra
 
#endif