draft vector ordering policy for SparseMatrix

include/micm/util/sparse_matrix.hpp

@@ -10,12 +10,12 @@
 #include <utility>
 #include <vector>
 
-#include <micm/util/standard_sparse_matrix.hpp>
+#include <micm/util/sparse_matrix_standard_ordering.hpp>
 
 namespace micm
 {
 
-  template<class T, class A>
+  template<class T, class OrderingPolicy>
   class SparseMatrixBuilder;
 
   /// @brief A sparse block-diagonal 2D matrix class with contiguous memory
@@ -26,15 +26,15 @@ namespace micm
   ///
   /// The template parameters are the type of the matrix elements and a class that
   /// defines the sizing and ordering of the data elements
-  template<class T, class A = StandardSparseMatrix>
-  class SparseMatrix : public A
+  template<class T, class OrderingPolicy = SparseMatrixStandardOrdering>
+  class SparseMatrix : public OrderingPolicy
   {
     std::size_t number_of_blocks_;        // Number of block sub-matrices in the overall matrix
     std::vector<std::size_t> row_ids_;    // Row indices of each non-zero element in a block
     std::vector<std::size_t> row_start_;  // Index in data_ and row_ids_ of the start of each row in a block
     std::vector<T> data_;                 // Value of each non-zero matrix element
 
-    friend class SparseMatrixBuilder<T, A>;
+    friend class SparseMatrixBuilder<T, OrderingPolicy>;
     friend class ProxyRow;
     friend class ConstProxyRow;
     friend class Proxy;
@@ -137,27 +137,27 @@ namespace micm
     };
 
    public:
-    static SparseMatrixBuilder<T, A> create(std::size_t block_size)
+    static SparseMatrixBuilder<T, OrderingPolicy> create(std::size_t block_size)
     {
-      return SparseMatrixBuilder<T, A>{ block_size };
+      return SparseMatrixBuilder<T, OrderingPolicy>{ block_size };
     }
 
     SparseMatrix() = default;
 
-    SparseMatrix(SparseMatrixBuilder<T, A>& builder)
+    SparseMatrix(SparseMatrixBuilder<T, OrderingPolicy>& builder)
         : number_of_blocks_(builder.number_of_blocks_),
           row_ids_(builder.RowIdsVector()),
           row_start_(builder.RowStartVector()),
-          data_(A::VectorSize(number_of_blocks_, row_ids_, row_start_), builder.initial_value_)
+          data_(OrderingPolicy::VectorSize(number_of_blocks_, row_ids_, row_start_), builder.initial_value_)
    {
     }
 
-    SparseMatrix<T, A>& operator=(SparseMatrixBuilder<T, A>& builder)
+    SparseMatrix<T, OrderingPolicy>& operator=(SparseMatrixBuilder<T, OrderingPolicy>& builder)
     {
       number_of_blocks_ = builder.number_of_blocks_;
       row_ids_ = builder.RowIdsVector();
       row_start_ = builder.RowStartVector();
-      data_ = std::vector<T>(A::VectorSize(number_of_blocks_, row_ids_, row_start_), builder.initial_value_);
+      data_ = std::vector<T>(OrderingPolicy::VectorSize(number_of_blocks_, row_ids_, row_start_), builder.initial_value_);
 
       return *this;
     }
@@ -173,7 +173,7 @@ namespace micm
 
     std::size_t VectorIndex(std::size_t block, std::size_t row, std::size_t column) const
     {
-      return A::VectorIndex(number_of_blocks_, row_ids_, row_start_, block, row, column);
+      return OrderingPolicy::VectorIndex(number_of_blocks_, row_ids_, row_start_, block, row, column);
     }
 
     std::size_t VectorIndex(std::size_t row, std::size_t column) const
@@ -226,14 +226,14 @@ namespace micm
     }
   };
 
-  template<class T, class A = StandardSparseMatrix>
+  template<class T, class OrderingPolicy = SparseMatrixStandardOrdering>
   class SparseMatrixBuilder
   {
     std::size_t number_of_blocks_{ 1 };
     std::size_t block_size_;
     std::set<std::pair<std::size_t, std::size_t>> non_zero_elements_{};
     T initial_value_{};
-    friend class SparseMatrix<T>;
+    friend class SparseMatrix<T, OrderingPolicy>;
 
    public:
     SparseMatrixBuilder() = delete;
@@ -243,26 +243,26 @@ namespace micm
     {
     }
 
-    operator SparseMatrix<T, A>() const
+    operator SparseMatrix<T, OrderingPolicy>() const
     {
-      return SparseMatrix<T, A>(*this);
+      return SparseMatrix<T, OrderingPolicy>(*this);
     }
 
-    SparseMatrixBuilder<T>& number_of_blocks(std::size_t n)
+    SparseMatrixBuilder<T, OrderingPolicy>& number_of_blocks(std::size_t n)
     {
       number_of_blocks_ = n;
       return *this;
     }
 
-    SparseMatrixBuilder<T>& with_element(std::size_t x, std::size_t y)
+    SparseMatrixBuilder<T, OrderingPolicy>& with_element(std::size_t x, std::size_t y)
     {
       if (x >= block_size_ || y >= block_size_)
         throw std::invalid_argument("SparseMatrix element out of range");
       non_zero_elements_.insert(std::make_pair(x, y));
       return *this;
     }
 
-    SparseMatrixBuilder<T>& initial_value(T inital_value)
+    SparseMatrixBuilder<T, OrderingPolicy>& initial_value(T inital_value)
     {
       initial_value_ = inital_value;
       return *this;

include/micm/util/standard_sparse_matrix.hpp → include/micm/util/sparse_matrix_standard_ordering.hpp

@@ -10,7 +10,7 @@ namespace micm
   ///
   /// Data is stored with blocks in the block diagonal matrix as the highest
   /// level structure, then by row, then by non-zero columns in each row.
-  class StandardSparseMatrix
+  class SparseMatrixStandardOrdering
   {
    protected:
     static std::size_t VectorSize(

include/micm/util/sparse_matrix_vector_ordering.hpp

@@ -0,0 +1,49 @@
+// Copyright (C) 2023 National Center for Atmospheric Research,
+//
+// SPDX-License-Identifier: Apache-2.0
+#pragma once
+
+#include <cmath>
+
+namespace micm
+{
+
+  /// @brief Defines the ordering of SparseMatrix object data to encourage vectorization
+  ///
+  /// Data is stored with sets of blocks in the block diagonal matrix as the highest
+  /// level structure, then by row, then by non-zero columns in each row, then by
+  /// individual blocks in the set of blocks.
+  ///
+  /// The template argument is the number of blocks per set of blocks and should be
+  /// approximately the size of the vector register.
+  template<std::size_t L>
+  class SparseMatrixVectorOrdering
+  {
+   protected:
+    static std::size_t VectorSize(
+        std::size_t number_of_blocks,
+        const std::vector<std::size_t>& row_ids,
+        const std::vector<std::size_t>& row_start)
+    {
+      return std::ceil((double)number_of_blocks / (double)L) * row_ids.size();
+    };
+
+    std::size_t VectorIndex(
+        std::size_t number_of_blocks,
+        const std::vector<std::size_t>& row_ids,
+        const std::vector<std::size_t>& row_start,
+        std::size_t block,
+        std::size_t row,
+        std::size_t column) const
+    {
+      if (row >= row_start.size() - 1 || column >= row_start.size() - 1 || block >= number_of_blocks)
+        throw std::invalid_argument("SparseMatrix element out of range");
+      auto begin = std::next(row_ids.begin(), row_start[row]);
+      auto end = std::next(row_ids.begin(), row_start[row + 1]);
+      auto elem = std::find(begin, end, column);
+      if (elem == end)
+        throw std::invalid_argument("SparseMatrix zero element access not allowed");
+      return std::size_t{ (elem - row_ids.begin()) * L + block % L + (block / L) * row_ids.size() };
+    };
+  };
+}  // namespace micm

test/unit/util/CMakeLists.txt

@@ -7,5 +7,6 @@ include(test_util)
 # Tests
 
 create_standard_test(NAME matrix SOURCES test_matrix.cpp)
-create_standard_test(NAME sparse_matrix SOURCES test_sparse_matrix.cpp)
+create_standard_test(NAME sparse_matrix_standard_ordering SOURCES test_sparse_matrix_standard_ordering.cpp)
+create_standard_test(NAME sparse_matrix_vector_ordering SOURCES test_sparse_matrix_vector_ordering.cpp)
 create_standard_test(NAME vector_matrix SOURCES test_vector_matrix.cpp)

test/unit/util/test_sparse_matrix_policy.hpp

@@ -0,0 +1,91 @@
+#include <gtest/gtest.h>
+
+#include <micm/util/sparse_matrix.hpp>
+
+template<class OrderingPolicy>
+micm::SparseMatrix<double, OrderingPolicy> testZeroMatrix()
+{
+  auto builder = micm::SparseMatrix<double, OrderingPolicy>::create(3);
+  auto row_ids = builder.RowIdsVector();
+  auto row_starts = builder.RowStartVector();
+
+  EXPECT_EQ(builder.NumberOfElements(), 0);
+  EXPECT_EQ(row_ids.size(), 0);
+  EXPECT_EQ(row_starts.size(), 4);
+  EXPECT_EQ(row_starts[0], 0);
+  EXPECT_EQ(row_starts[1], 0);
+  EXPECT_EQ(row_starts[2], 0);
+  EXPECT_EQ(row_starts[3], 0);
+
+  micm::SparseMatrix<double, OrderingPolicy> matrix{ builder };
+
+  EXPECT_EQ(matrix.FlatBlockSize(), 0);
+
+  EXPECT_THROW(
+      try { std::size_t elem = matrix.VectorIndex(0, 0); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix zero element access not allowed");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { std::size_t elem = matrix.VectorIndex(6, 0); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { std::size_t elem = matrix.VectorIndex(1, 3); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { std::size_t elem = matrix.VectorIndex(6, 3); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { bool isZero = matrix.IsZero(6, 0); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { bool isZero = matrix.IsZero(1, 3); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { bool isZero = matrix.IsZero(6, 3); } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { matrix[0][0][4] = 2.0; } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { matrix[1][0][0] = 2.0; } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { matrix[0][3][0] = 2.0; } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
+        throw;
+      },
+      std::invalid_argument);
+  EXPECT_THROW(
+      try { matrix[0][1][1] = 2.0; } catch (const std::invalid_argument& e) {
+        EXPECT_STREQ(e.what(), "SparseMatrix zero element access not allowed");
+        throw;
+      },
+      std::invalid_argument);
+  return matrix;
+}

test/unit/util/test_sparse_matrix.cpp → test/unit/util/test_sparse_matrix_standard_ordering.cpp

@@ -1,91 +1,12 @@
 #include <gtest/gtest.h>
 
 #include <micm/util/sparse_matrix.hpp>
+#include <micm/util/sparse_matrix_standard_ordering.hpp>
+#include "test_sparse_matrix_policy.hpp"
 
 TEST(SparseMatrix, ZeroMatrix)
 {
-  auto builder = micm::SparseMatrix<double>::create(3);
-  auto row_ids = builder.RowIdsVector();
-  auto row_starts = builder.RowStartVector();
-
-  EXPECT_EQ(builder.NumberOfElements(), 0);
-  EXPECT_EQ(row_ids.size(), 0);
-  EXPECT_EQ(row_starts.size(), 4);
-  EXPECT_EQ(row_starts[0], 0);
-  EXPECT_EQ(row_starts[1], 0);
-  EXPECT_EQ(row_starts[2], 0);
-  EXPECT_EQ(row_starts[3], 0);
-
-  micm::SparseMatrix<double> matrix{ builder };
-
-  EXPECT_EQ(matrix.FlatBlockSize(), 0);
-
-  EXPECT_THROW(
-      try { std::size_t elem = matrix.VectorIndex(0, 0); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix zero element access not allowed");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { std::size_t elem = matrix.VectorIndex(6, 0); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { std::size_t elem = matrix.VectorIndex(1, 3); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { std::size_t elem = matrix.VectorIndex(6, 3); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { bool isZero = matrix.IsZero(6, 0); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { bool isZero = matrix.IsZero(1, 3); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { bool isZero = matrix.IsZero(6, 3); } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { matrix[0][0][4] = 2.0; } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { matrix[1][0][0] = 2.0; } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { matrix[0][3][0] = 2.0; } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix element out of range");
-        throw;
-      },
-      std::invalid_argument);
-  EXPECT_THROW(
-      try { matrix[0][1][1] = 2.0; } catch (const std::invalid_argument& e) {
-        EXPECT_STREQ(e.what(), "SparseMatrix zero element access not allowed");
-        throw;
-      },
-      std::invalid_argument);
+  testZeroMatrix<micm::SparseMatrixStandardOrdering>();
 }
 
 TEST(SparseMatrix, ConstZeroMatrix)

test/unit/util/test_sparse_matrix_vector_ordering.cpp

@@ -0,0 +1,10 @@
+#include <gtest/gtest.h>
+
+#include <micm/util/sparse_matrix.hpp>
+#include <micm/util/sparse_matrix_vector_ordering.hpp>
+#include "test_sparse_matrix_policy.hpp"
+
+TEST(SparseMatrix, ZeroMatrix)
+{
+  testZeroMatrix<micm::SparseMatrixVectorOrdering<2>>();
+}