CsrMatrix.h

Go to the documentation of this file.

#ifndef CSR_MATRIX_H
#define CSR_MATRIX_H
 
#include <algorithm>
#include <iostream>
#include <tuple>
#include <valarray>
#include <vector>
 
#include "CoreEnums.h"
#include "GenericMatrix.h"
 
class StaticSparseMatrixAccessError : public std::runtime_error
{
    public:
        StaticSparseMatrixAccessError(const int row, const int col)
            : std::runtime_error("Accessing zero element at (" + std::to_string(row) + ","
                                 + std::to_string(col) + ") in non-const context.")
        {
        }
};
 
template <typename T>
class CsrSparseMatrix : public GenericMatrix<T>
{
    private:
        std::valarray<int> m_column_indices;
        std::valarray<int> m_row_pointers;
        int m_number_non_zeros{};
        bool m_dynamic_creation{};
        std::vector<std::tuple<T, int, int>> m_temporary_dynamic_create_storage{};
 
        int findIndex(int rows, int cols) const;
 
    public:
        CsrSparseMatrix();
        CsrSparseMatrix(const std::valarray<T>& values, const std::valarray<int>& column_indices,
                        const std::valarray<int>& row_pointers, int rows, int cols);
 
        void registerElement(T value, int rows, int cols) override;
        void buildAtRuntime();
        void buildAtRuntime(int rows, int cols) override;
 
        T& operator()(int rows, int cols) override;
        T operator()(int rows, int cols) const override;
 
        std::valarray<T> scalarAdd(T value) const override;
        std::valarray<T> scalarSubtract(T value) const override;
        std::valarray<T> scalarMultiply(T value) const override;
        std::valarray<T> scalarDivide(T value) const override;
 
        std::valarray<T> vectorElemAdd(const std::valarray<T>& value) const override;
        std::valarray<T> vectorElemSubtract(const std::valarray<T>& value) const override;
        std::valarray<T> vectorElemMultiply(const std::valarray<T>& value) const override;
        std::valarray<T> vectorElemDivide(const std::valarray<T>& value) const override;
 
        std::valarray<T> matrixVectorMultiplication(const std::valarray<T>& value) const override;
 
        const MatrixType getMatrixType() override { return MatrixType::Csr; }
 
        void printMatrixInDenseFormat();
 
        const std::vector<std::pair<int, int>>& getNonZeroElements() override
        {
            return this->m_non_zero_elements;
        };
 
        int getNonZeroElementNumber() const { return this->m_matrix.size(); }
        const std::valarray<int>& getColumnIndices() { return this->m_column_indices; }
        const std::valarray<int>& getRowPointers() { return this->m_row_pointers; }
        const std::valarray<T>& getMatrixElements() { return this->m_matrix; }
        void updateNonZeroElements() override;
 
        ~CsrSparseMatrix() override = default;
};
 
template <typename T>
CsrSparseMatrix<T>::CsrSparseMatrix() : m_dynamic_creation{true}
{
}
 
template <typename T>
CsrSparseMatrix<T>::CsrSparseMatrix(const std::valarray<T>& values,
                                    const std::valarray<int>& column_indices,
                                    const std::valarray<int>& row_pointers, int rows, int cols)
    : GenericMatrix<T>{rows, cols, values, true},
      m_column_indices{column_indices},
      m_row_pointers{row_pointers}
{
}
 
template <typename T>
void CsrSparseMatrix<T>::registerElement(const T value, const int rows, const int cols)
{
    if (m_dynamic_creation)
    {
        m_temporary_dynamic_create_storage.emplace_back(value, rows, cols);
    }
}
 
template <typename T>
void CsrSparseMatrix<T>::buildAtRuntime()
{
    if (m_dynamic_creation)
    {
        const auto& max_row_ptr = std::max_element(
            m_temporary_dynamic_create_storage.begin(), m_temporary_dynamic_create_storage.end(),
            [](std::tuple<T, int, int> lhs, std::tuple<T, int, int> rhs)
            { return std::get<1>(lhs) < std::get<1>(rhs); });
        const auto& max_col_ptr = std::max_element(
            m_temporary_dynamic_create_storage.begin(), m_temporary_dynamic_create_storage.end(),
            [](std::tuple<T, int, int> lhs, std::tuple<T, int, int> rhs)
            { return std::get<2>(lhs) < std::get<2>(rhs); });
        this->m_rows = std::get<1>(*max_row_ptr) + 1;
        this->m_cols = std::get<2>(*max_col_ptr) + 1;
        // Do not use the max_ptrs afterwards, as we sort the elements and the pointers will lose
        // their max property. Sort the temporary storage by rows (ascending) and within a row by
        // columns (ascending).
        std::sort(m_temporary_dynamic_create_storage.begin(),
                  m_temporary_dynamic_create_storage.end(),
                  [](std::tuple<T, int, int> lhs, std::tuple<T, int, int> rhs)
                  {
                      if (std::get<1>(lhs) == std::get<1>(rhs))
                      {
                          return std::get<2>(lhs) < std::get<2>(rhs);
                      }
                      return std::get<1>(lhs) < std::get<1>(rhs);
                  });
        // Clear out duplicate entries. They can be produced when building a global matrix from
        // smaller ones that overlap in some entries, like when using finite elements.
        for (auto it = m_temporary_dynamic_create_storage.begin();;)
        {
            if (std::get<1>(*it) == std::get<1>(*(it + 1))
                && std::get<2>(*it) == std::get<2>(*(it + 1)))
            {
                it = m_temporary_dynamic_create_storage.erase(it);
            }
            else
            {
                ++it;
            }
            if (it + 1 == m_temporary_dynamic_create_storage.end())
            {
                break;
            }
        }
        // All non-zero entries are known, so the valarray sizes can be finalized
        this->m_matrix.resize(m_temporary_dynamic_create_storage.size());
        m_column_indices.resize(m_temporary_dynamic_create_storage.size());
        // The row pointer always has one element more that points to n+1, as to always define
        // row_pointer[i+1].
        m_row_pointers.resize(this->m_rows + 1);
 
        auto current_row_start{m_temporary_dynamic_create_storage.begin()};
        bool is_last{false};
        // We need the +1 on the end iterator here, because we need the case it = .end() within the
        // if condition
        for (auto it{m_temporary_dynamic_create_storage.begin()};; it++)
        {
            is_last = (it == m_temporary_dynamic_create_storage.end());
            if (is_last || std::get<1>(*it) != std::get<1>(*current_row_start))
            {
                // Set the row pointer for the current row to the element, which marks the start of
                // that row.
                m_row_pointers[std::get<1>(*current_row_start)] = static_cast<int>(
                    current_row_start - m_temporary_dynamic_create_storage.begin());
                // If there are empty row(s) between the entries, we need to propagate the next
                // valid row_pointer value, so that row[i+1]-row[i] = 0 for these empty rows. As
                // this is not calculated, we mark it with a placeholder and fill it after going
                // through it once. Not is_last is checked first in order to not deref the end() 
                // pointer.
                if (!is_last && std::get<1>(*it) > std::get<1>(*current_row_start) + 1)
                {
                    for (int i{std::get<1>(*current_row_start) + 1}; i < std::get<1>(*it); i++)
                    {
                        m_row_pointers[i] = -1;
                    }
                }
                current_row_start = it;
            }
            if (is_last)
            {
                break;
            }
        }
        // Now that all values are calculated, we can fill in the placeholders.
        for (int i{static_cast<int>(m_row_pointers.size()) - 1}; i > 0; i--)
        {
            if (m_row_pointers[i - 1] == -1)
            {
                m_row_pointers[i - 1] = m_row_pointers[i];
            }
        }
        // This element is not a part of the matrix, but it makes the expresion m_row_pointers[rows]
        // and m_row_pointers[rows + 1] valid for all rows <= nr_of_rows!
        // m_row_pointers[current_row_int] = static_cast<int>( current_row_start -
        // m_temporary_dynamic_create_storage.begin() ); current_row_int++;
        m_row_pointers[this->m_rows] = this->m_matrix.size();
 
        // After this block, the elements are all sorted by rows (ascending) and within one row by
        // columns (ascending) and the beginnings of the rows are already marked, so we can simply
        // write them directly to the arrays.
        for (int i{0}; i < m_temporary_dynamic_create_storage.size(); i++)
        {
            this->m_matrix[i] = std::get<0>(m_temporary_dynamic_create_storage[i]);
            m_column_indices[i] = std::get<2>(m_temporary_dynamic_create_storage[i]);
        }
    }
}
 
template <typename T>
void CsrSparseMatrix<T>::buildAtRuntime(const int rows, const int cols)
{
    if (m_dynamic_creation)
    {
        this->m_rows = rows;
        this->m_cols = cols;
        // Do not use the max_ptrs afterwards, as we sort the elements and the pointers will lose
        // their max property. Sort the temporary storage by rows (ascending) and within a row by
        // columns (ascending).
        std::sort(m_temporary_dynamic_create_storage.begin(),
                  m_temporary_dynamic_create_storage.end(),
                  [](std::tuple<T, int, int> lhs, std::tuple<T, int, int> rhs)
                  {
                      if (std::get<1>(lhs) == std::get<1>(rhs))
                      {
                          return std::get<2>(lhs) < std::get<2>(rhs);
                      }
                      return std::get<1>(lhs) < std::get<1>(rhs);
                  });
        // Clear out duplicate entries. They can be produced when building a global matrix from
        // smaller ones that overlap in some entries, like when using finite elements.
        for (auto it = m_temporary_dynamic_create_storage.begin();;)
        {
            if (std::get<1>(*it) == std::get<1>(*(it + 1))
                && std::get<2>(*it) == std::get<2>(*(it + 1)))
            {
                it = m_temporary_dynamic_create_storage.erase(it);
            }
            else
            {
                ++it;
            }
            if (it + 1 == m_temporary_dynamic_create_storage.end())
            {
                break;
            }
        }
        // All non-zero entries are known, so the valarray sizes can be finalized
        this->m_matrix.resize(m_temporary_dynamic_create_storage.size());
        m_column_indices.resize(m_temporary_dynamic_create_storage.size());
        // The last row always has entries for all FEM/FV/FD discretizations, so we can set the row
        // pointer array size to this value + 2
        m_row_pointers.resize(rows + 1);
 
        auto current_row_start{m_temporary_dynamic_create_storage.begin()};
        bool is_last{false};
        bool last_is_empty{false};
        // We need the +1 on the end iterator here, because we need the case it = .end() within the
        // if condition
        for (auto it{m_temporary_dynamic_create_storage.begin()};; it++)
        {
            is_last = (it == m_temporary_dynamic_create_storage.end());
            if (is_last || std::get<1>(*it) != std::get<1>(*current_row_start))
            {
                // Set the row pointer for the current row to the element, which marks the start of
                // that row.
                m_row_pointers[std::get<1>(*current_row_start)] = static_cast<int>(
                    current_row_start - m_temporary_dynamic_create_storage.begin());
                // If there are empty row(s) between the entries, we need to propagate the next
                // valid row_pointer value, so that row[i+1]-row[i] = 0 for these empty rows. As
                // this is not calculated, we mark it with a placeholder and fill it after going
                // through it once. Not is_last is checked first in order to not deref the end() 
                // pointer.
                if (!is_last && std::get<1>(*it) > std::get<1>(*current_row_start) + 1)
                {
                    for (int i{std::get<1>(*current_row_start) + 1}; i < std::get<1>(*it); i++)
                    {
                        m_row_pointers[i] = -1;
                    }
                }
                // If the last row was not the nomially last one, it is empty and needs to be 
                // treated differently.
                if(is_last && std::get<1>(*current_row_start) != rows-1)
                {
                    last_is_empty = true;
                }
                current_row_start = it;
            }
            if (is_last)
            {
                break;
            }
        }
        // Now that all values are calculated, we can fill in the placeholders.
        if(last_is_empty)
        {
            m_row_pointers[rows-1] = this->m_matrix.size();
        }
        for (int i{static_cast<int>(m_row_pointers.size()) - 1}; i > 0; i--)
        {
            if (m_row_pointers[i - 1] == -1)
            {
                m_row_pointers[i - 1] = m_row_pointers[i];
            }
        }
        // This element is not a part of the matrix, but it makes the expresion m_row_pointers[rows]
        // and m_row_pointers[rows + 1] valid for all rows <= nr_of_rows!
        // m_row_pointers[current_row_int] = static_cast<int>( current_row_start -
        // m_temporary_dynamic_create_storage.begin() ); current_row_int++;
        m_row_pointers[rows] = this->m_matrix.size();
 
        // After this block, the elements are all sorted by rows (ascending) and within one row by
        // columns (ascending) and the beginnings of the rows are already marked, so we can simply
        // write them directly to the arrays.
        for (int i{0}; i < m_temporary_dynamic_create_storage.size(); i++)
        {
            this->m_matrix[i] = std::get<0>(m_temporary_dynamic_create_storage[i]);
            m_column_indices[i] = std::get<2>(m_temporary_dynamic_create_storage[i]);
        }
    }
}
 
template <typename T>
int CsrSparseMatrix<T>::findIndex(const int rows, const int cols) const
{
    int lower_bounds{m_row_pointers[rows]};
    int upper_bounds{m_row_pointers[rows + 1] - 1};
    int test_index{-1};
 
    // Use devide and conquer algorithm to find result in log(n) time
    while (lower_bounds <= upper_bounds)
    {
        test_index = ((lower_bounds + upper_bounds) / 2);
        if (m_column_indices[test_index] < cols)
        {
            lower_bounds = test_index + 1;
        }
        else if (m_column_indices[test_index] > cols)
        {
            upper_bounds = test_index - 1;
        }
        else
        {
            return test_index;
        }
    }
    // If no match was found, then return a negative value to show this.
    return -1;
}
 
template <typename T>
T& CsrSparseMatrix<T>::operator()(const int rows, const int cols)
{
    if (rows >= this->m_rows || cols >= this->m_cols)
    {
        throw MatrixOutOfBoundsError(rows, cols, this->m_rows, this->m_cols);
    }
 
    int matrix_value_index{findIndex(rows, cols)};
    // If the value is non-zero, is has to be in the given row, whose values are all between
    // m_row_pointers[rows] and m_row_pointers[rows + 1]
 
    // If _find_index returned -1, it has to be a zero value.
    // This is technically not intended to access 0 elements in non-const context, as we would need
    // to write it into the matrix (which schould never carry elements that are zero) If we want to
    // allow element insertion in place (warning, this can be VERY slow), then we can replace that
    // here. Right now it's just a placeholder
    if (matrix_value_index == -1)
    {
        throw StaticSparseMatrixAccessError(rows, cols);
    }
    return this->m_matrix[matrix_value_index];
}
 
template <typename T>
T CsrSparseMatrix<T>::operator()(const int rows, const int cols) const
{
    if (rows >= this->m_rows || cols >= this->m_cols || rows < 0 || cols < 0)
    {
        throw MatrixOutOfBoundsError(rows, cols, this->m_rows, this->m_cols);
    }
 
    int matrix_value_index{findIndex(rows, cols)};
 
    // If the value is a non-defined zero value, return zero.
    if (matrix_value_index == -1)
    {
        throw static_cast<T>(0.0);
    }
    return this->m_matrix[matrix_value_index];
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::scalarAdd(const T value) const
{
    return this->m_matrix + value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::scalarSubtract(const T value) const
{
    return this->m_matrix - value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::scalarMultiply(const T value) const
{
    return this->m_matrix * value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::scalarDivide(const T value) const
{
    return this->m_matrix / value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::vectorElemAdd(const std::valarray<T>& value) const
{
    return this->m_matrix + value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::vectorElemSubtract(const std::valarray<T>& value) const
{
    return this->m_matrix - value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::vectorElemMultiply(const std::valarray<T>& value) const
{
    return this->m_matrix * value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::vectorElemDivide(const std::valarray<T>& value) const
{
    return this->m_matrix / value;
}
 
template <typename T>
std::valarray<T> CsrSparseMatrix<T>::matrixVectorMultiplication(const std::valarray<T>& value) const
{
    std::valarray<T> result(this->m_rows);
    int index_row_start;
    int index_row_end;
    for (int i{0}; i < this->m_rows; i++)
    {
        index_row_start = m_row_pointers[i];
        index_row_end = m_row_pointers[i + 1];
        for (int j{index_row_start}; j < index_row_end; j++)
        {
            result[i] += this->m_matrix[j] * value[m_column_indices[j]];
        }
    }
    return result;
}
 
template <typename T>
void CsrSparseMatrix<T>::printMatrixInDenseFormat()
{
    if (this->m_matrix.size() > 0)
    {
        for (int i{0}; i < m_row_pointers.size() - 1; i++)
        {
            std::valarray<int> column_values_in_row = m_column_indices[std::slice(
                m_row_pointers[i], m_row_pointers[i + 1] - m_row_pointers[i], 1)];
            int current_column_id{0};
            for (int j{0}; j < this->m_cols; j++)
            {
                if (current_column_id < column_values_in_row.size()
                    && j == column_values_in_row[current_column_id])
                {
                    std::cout << this->m_matrix[m_row_pointers[i] + current_column_id] << ' ';
                    current_column_id++;
                }
                else
                {
                    std::cout << "0" << ' ';
                }
            }
            std::cout << "\n";
        }
    }
    else
    {
        std::cout << "Matrix is empty\n";
    }
}
 
template <typename T>
void CsrSparseMatrix<T>::updateNonZeroElements()
{
    for (const auto& position_tuple : m_temporary_dynamic_create_storage)
    {
        this->m_non_zero_elements.emplace_back(std::get<1>(position_tuple),
                                               std::get<2>(position_tuple));
    }
}
#endif