diff --git a/.gitignore b/.gitignore index f5fd378..40c0a25 100644 --- a/.gitignore +++ b/.gitignore @@ -54,4 +54,4 @@ build* *.obj *.dll *.lib -*.pdb \ No newline at end of file +*.pdb_codeql_* diff --git a/include/CSRComm.h b/include/CSRComm.h index 60e5982..855f4dc 100644 --- a/include/CSRComm.h +++ b/include/CSRComm.h @@ -3,6 +3,7 @@ #include "Global.h" #include "CSRMatrix.h" #include "CommPattern.h" +#include #include #include @@ -136,6 +137,7 @@ class cCSRComm { } // Accumulate received contributions into local CSR matrix. + // Use binary search for better performance (O(log n) instead of O(n) per lookup) for (const cTriplet& cT : vRecvBuf) { iIndex iGlobalRow = cT.m_iRow; iIndex iGlobalCol = cT.m_iCol; @@ -144,11 +146,14 @@ class cCSRComm { iIndex iStart = vRowPtr[static_cast(iLocalRow)]; iIndex iEnd = vRowPtr[static_cast(iLocalRow + 1)]; - for (iIndex i = iStart; i < iEnd; ++i) { - if (vColInd[static_cast(i)] == iGlobalCol) { - vLocalVal[static_cast(i)] += cT.m_vVal; - break; - } + // Binary search within the row's column indices + auto itBegin = vColInd.begin() + iStart; + auto itEnd = vColInd.begin() + iEnd; + auto it = std::lower_bound(itBegin, itEnd, iGlobalCol); + + if (it != itEnd && *it == iGlobalCol) { + std::size_t idx = static_cast(it - vColInd.begin()); + vLocalVal[idx] += cT.m_vVal; } } } diff --git a/include/CSRMatrixBuilder.h b/include/CSRMatrixBuilder.h new file mode 100644 index 0000000..2772697 --- /dev/null +++ b/include/CSRMatrixBuilder.h @@ -0,0 +1,245 @@ +#pragma once + +#include "CSRMatrix.h" +#include "Global.h" +#include +#include +#include + +namespace csr4mpi { + +template +struct cTriplet { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + iIndex m_iRow; + iIndex m_iCol; + Scalar m_vValue; + + cTriplet() : m_iRow(0), m_iCol(0), m_vValue(static_cast(0)) {} + cTriplet(iIndex iRow, iIndex iCol, Scalar vValue) + : m_iRow(iRow), m_iCol(iCol), m_vValue(vValue) {} +}; + +// Type aliases for common scalar types +using cTripletF = cTriplet; +using cTripletD = cTriplet; +using cTripletCF = cTriplet>; +using cTripletCD = cTriplet>; + +template +class cCSRMatrixBuilder { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + +public: + using scalar_type = Scalar; + + cCSRMatrixBuilder() = default; + + // Add a single triplet entry + void AddEntry(iIndex iRow, iIndex iCol, Scalar vValue) + { + m_vTriplets.push_back(cTriplet(iRow, iCol, vValue)); + } + + // Add multiple triplet entries + void AddEntries(const std::vector>& vTriplets) + { + m_vTriplets.insert(m_vTriplets.end(), vTriplets.begin(), vTriplets.end()); + } + + // Reserve space for efficiency + void Reserve(std::size_t size) + { + m_vTriplets.reserve(size); + } + + // Clear all entries + void Clear() + { + m_vTriplets.clear(); + } + + // Build the CSR matrix with automatic deduplication + // iGlobalRowBegin: starting global row index (typically 0 for a standalone local matrix) + // iGlobalRowEnd: ending global row index (exclusive) + // iGlobalColCount: total number of columns + // bAccumulateDuplicates: if true, accumulate duplicate entries; if false, keep last value + cCSRMatrix Build(iSize iGlobalRowBegin, iSize iGlobalRowEnd, iSize iGlobalColCount, bool bAccumulateDuplicates = true) const + { + struct cKey { + iIndex r; + iIndex c; + bool operator==(const cKey& other) const { return r == other.r && c == other.c; } + }; + struct cKeyHash { + std::size_t operator()(const cKey& k) const noexcept + { + // Use a robust hash combining function to reduce collisions + // Based on boost::hash_combine algorithm + std::size_t seed = std::hash()(k.r); + seed ^= std::hash()(k.c) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + return seed; + } + }; + + iSize iLocalRows = iGlobalRowEnd - iGlobalRowBegin; + + // Deduplicate using hash map + std::unordered_map mEntries; + mEntries.reserve(m_vTriplets.size()); + + for (const auto& triplet : m_vTriplets) { + // Validate row is within range + if (triplet.m_iRow < iGlobalRowBegin || triplet.m_iRow >= iGlobalRowEnd) { + continue; // Skip out-of-range entries + } + + cKey key { triplet.m_iRow, triplet.m_iCol }; + auto it = mEntries.find(key); + if (it == mEntries.end()) { + mEntries.emplace(key, triplet.m_vValue); + } else { + if (bAccumulateDuplicates) { + it->second += triplet.m_vValue; + } else { + it->second = triplet.m_vValue; // Keep last value + } + } + } + + // Convert to sorted triplets + struct cEntry { + iIndex r; + iIndex c; + Scalar v; + }; + std::vector vSortedEntries; + vSortedEntries.reserve(mEntries.size()); + + for (const auto& p : mEntries) { + vSortedEntries.push_back({ p.first.r, p.first.c, p.second }); + } + + // Sort by row then column + std::sort(vSortedEntries.begin(), vSortedEntries.end(), [](const cEntry& a, const cEntry& b) { + if (a.r != b.r) + return a.r < b.r; + return a.c < b.c; + }); + + // Build CSR arrays + std::vector vRowPtr(static_cast(iLocalRows + 1), 0); + std::vector vColInd; + std::vector vValues; + + vColInd.reserve(vSortedEntries.size()); + vValues.reserve(vSortedEntries.size()); + + for (const auto& e : vSortedEntries) { + iIndex iLocalRow = e.r - iGlobalRowBegin; + vRowPtr[static_cast(iLocalRow + 1)]++; + vColInd.push_back(e.c); + vValues.push_back(e.v); + } + + // Prefix sum to get row pointers + for (iSize r = 0; r < iLocalRows; ++r) { + vRowPtr[static_cast(r + 1)] += vRowPtr[static_cast(r)]; + } + + return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, + std::move(vRowPtr), std::move(vColInd), std::move(vValues)); + } + +private: + std::vector> m_vTriplets; +}; + +// Type aliases for common scalar types +using cCSRMatrixBuilderF = cCSRMatrixBuilder; +using cCSRMatrixBuilderD = cCSRMatrixBuilder; +using cCSRMatrixBuilderCF = cCSRMatrixBuilder>; +using cCSRMatrixBuilderCD = cCSRMatrixBuilder>; + +// Convenience function to build CSR matrix directly from triplets +template +cCSRMatrix BuildCSRFromTriplets( + const std::vector>& vTriplets, + iSize iGlobalRowBegin, + iSize iGlobalRowEnd, + iSize iGlobalColCount, + bool bAccumulateDuplicates = true) +{ + cCSRMatrixBuilder builder; + builder.AddEntries(vTriplets); + return builder.Build(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, bAccumulateDuplicates); +} + +// Validation function to check if CSR matrix structure is valid +template +bool ValidateCSRMatrix(const cCSRMatrix& matrix, std::string* pErrorMsg = nullptr) +{ + const auto& rowPtr = matrix.vRowPtr(); + const auto& colInd = matrix.vColInd(); + const auto& values = matrix.vValues(); + + iSize localRows = matrix.iGlobalRowEnd() - matrix.iGlobalRowBegin(); + + // Check row pointer size + if (rowPtr.size() != static_cast(localRows + 1)) { + if (pErrorMsg) *pErrorMsg = "Row pointer size mismatch"; + return false; + } + + // Check row pointers are non-decreasing + for (iSize i = 0; i < localRows; ++i) { + if (rowPtr[static_cast(i)] > rowPtr[static_cast(i + 1)]) { + if (pErrorMsg) *pErrorMsg = "Row pointers are not non-decreasing"; + return false; + } + } + + // Check first row pointer is 0 + if (rowPtr[0] != 0) { + if (pErrorMsg) *pErrorMsg = "First row pointer is not 0"; + return false; + } + + // Check last row pointer matches number of non-zeros + if (rowPtr[static_cast(localRows)] != static_cast(colInd.size())) { + if (pErrorMsg) *pErrorMsg = "Last row pointer does not match column index size"; + return false; + } + + // Check column indices and values have same size + if (colInd.size() != values.size()) { + if (pErrorMsg) *pErrorMsg = "Column index and value array size mismatch"; + return false; + } + + // Check column indices are sorted within each row and within bounds + for (iSize r = 0; r < localRows; ++r) { + iIndex start = rowPtr[static_cast(r)]; + iIndex end = rowPtr[static_cast(r + 1)]; + + for (iIndex k = start; k < end; ++k) { + iIndex col = colInd[static_cast(k)]; + + // Check column is within bounds + if (col < 0 || col >= matrix.iGlobalColCount()) { + if (pErrorMsg) *pErrorMsg = "Column index out of bounds"; + return false; + } + + // Check columns are sorted within row + if (k > start && colInd[static_cast(k - 1)] >= col) { + if (pErrorMsg) *pErrorMsg = "Column indices are not sorted within row"; + return false; + } + } + } + + return true; +} + +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ca0c31f..c9391cb 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,6 +31,8 @@ endif() add_executable(csr4mpi_tests test_csr_matrix.cpp + test_csr_builder.cpp + test_csrcomm_optimization.cpp test_distribution_pattern.cpp test_mumps_adapter.cpp test_operations.cpp diff --git a/tests/test_csr_builder.cpp b/tests/test_csr_builder.cpp new file mode 100644 index 0000000..37f9b0e --- /dev/null +++ b/tests/test_csr_builder.cpp @@ -0,0 +1,451 @@ +#include "CSRMatrixBuilder.h" +#include "CSRMatrix.h" +#include "Operations.h" +#include +#include + +using namespace csr4mpi; + +// Test with double as the default scalar type +using Scalar = double; + +TEST(CSRBuilderTest, BasicConstruction) +{ + cCSRMatrixBuilder builder; + + // Build a simple 3x3 matrix: + // [1 2 0] + // [0 3 4] + // [5 0 6] + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(0, 1, 2.0); + builder.AddEntry(1, 1, 3.0); + builder.AddEntry(1, 2, 4.0); + builder.AddEntry(2, 0, 5.0); + builder.AddEntry(2, 2, 6.0); + + auto matrix = builder.Build(0, 3, 3); + + EXPECT_EQ(matrix.iGlobalRowBegin(), 0); + EXPECT_EQ(matrix.iGlobalRowEnd(), 3); + EXPECT_EQ(matrix.iGlobalColCount(), 3); + EXPECT_EQ(matrix.vRowPtr().size(), 4u); + EXPECT_EQ(matrix.vColInd().size(), 6u); + EXPECT_EQ(matrix.vValues().size(), 6u); + + // Verify structure + const auto& rowPtr = matrix.vRowPtr(); + const auto& colInd = matrix.vColInd(); + const auto& values = matrix.vValues(); + + EXPECT_EQ(rowPtr[0], 0); + EXPECT_EQ(rowPtr[1], 2); + EXPECT_EQ(rowPtr[2], 4); + EXPECT_EQ(rowPtr[3], 6); + + // Row 0: cols 0,1 + EXPECT_EQ(colInd[0], 0); + EXPECT_EQ(colInd[1], 1); + EXPECT_DOUBLE_EQ(values[0], 1.0); + EXPECT_DOUBLE_EQ(values[1], 2.0); + + // Row 1: cols 1,2 + EXPECT_EQ(colInd[2], 1); + EXPECT_EQ(colInd[3], 2); + EXPECT_DOUBLE_EQ(values[2], 3.0); + EXPECT_DOUBLE_EQ(values[3], 4.0); + + // Row 2: cols 0,2 + EXPECT_EQ(colInd[4], 0); + EXPECT_EQ(colInd[5], 2); + EXPECT_DOUBLE_EQ(values[4], 5.0); + EXPECT_DOUBLE_EQ(values[5], 6.0); +} + +TEST(CSRBuilderTest, DuplicateAccumulation) +{ + cCSRMatrixBuilder builder; + + // Add duplicates that should accumulate + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(0, 0, 2.0); + builder.AddEntry(0, 0, 3.0); + builder.AddEntry(1, 1, 5.0); + builder.AddEntry(1, 1, 7.0); + + auto matrix = builder.Build(0, 2, 2, true); // Accumulate duplicates + + const auto& rowPtr = matrix.vRowPtr(); + const auto& colInd = matrix.vColInd(); + const auto& values = matrix.vValues(); + + EXPECT_EQ(colInd.size(), 2u); + EXPECT_EQ(values.size(), 2u); + + // Row 0: (0,0) should be 1+2+3=6 + EXPECT_EQ(colInd[0], 0); + EXPECT_DOUBLE_EQ(values[0], 6.0); + + // Row 1: (1,1) should be 5+7=12 + EXPECT_EQ(colInd[1], 1); + EXPECT_DOUBLE_EQ(values[1], 12.0); +} + +TEST(CSRBuilderTest, DuplicateLastValue) +{ + cCSRMatrixBuilder builder; + + // Add duplicates, but keep only the last value + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(0, 0, 2.0); + builder.AddEntry(0, 0, 3.0); + + auto matrix = builder.Build(0, 1, 1, false); // Keep last value + + const auto& values = matrix.vValues(); + + EXPECT_EQ(values.size(), 1u); + EXPECT_DOUBLE_EQ(values[0], 3.0); // Should be the last value +} + +TEST(CSRBuilderTest, EmptyMatrix) +{ + cCSRMatrixBuilder builder; + + // Build without adding any entries + auto matrix = builder.Build(0, 5, 5); + + EXPECT_EQ(matrix.iGlobalRowBegin(), 0); + EXPECT_EQ(matrix.iGlobalRowEnd(), 5); + EXPECT_EQ(matrix.iGlobalColCount(), 5); + EXPECT_EQ(matrix.vRowPtr().size(), 6u); + EXPECT_EQ(matrix.vColInd().size(), 0u); + EXPECT_EQ(matrix.vValues().size(), 0u); + + // All row pointers should be 0 + const auto& rowPtr = matrix.vRowPtr(); + for (std::size_t i = 0; i < rowPtr.size(); ++i) { + EXPECT_EQ(rowPtr[i], 0); + } +} + +TEST(CSRBuilderTest, SingleElement) +{ + cCSRMatrixBuilder builder; + + builder.AddEntry(0, 0, 42.0); + + auto matrix = builder.Build(0, 1, 1); + + EXPECT_EQ(matrix.vColInd().size(), 1u); + EXPECT_EQ(matrix.vValues().size(), 1u); + EXPECT_DOUBLE_EQ(matrix.vValues()[0], 42.0); +} + +TEST(CSRBuilderTest, AllDuplicates) +{ + cCSRMatrixBuilder builder; + + // All entries are duplicates of the same element + for (int i = 0; i < 100; ++i) { + builder.AddEntry(2, 3, 1.0); + } + + auto matrix = builder.Build(0, 5, 5); + + EXPECT_EQ(matrix.vColInd().size(), 1u); + EXPECT_EQ(matrix.vValues().size(), 1u); + EXPECT_DOUBLE_EQ(matrix.vValues()[0], 100.0); +} + +TEST(CSRBuilderTest, UnsortedInput) +{ + cCSRMatrixBuilder builder; + + // Add entries in random order + builder.AddEntry(2, 2, 9.0); + builder.AddEntry(0, 1, 2.0); + builder.AddEntry(1, 0, 3.0); + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(2, 0, 7.0); + builder.AddEntry(1, 2, 6.0); + builder.AddEntry(1, 1, 5.0); + builder.AddEntry(2, 1, 8.0); + + auto matrix = builder.Build(0, 3, 3); + + // Should be sorted by row then column + const auto& colInd = matrix.vColInd(); + const auto& values = matrix.vValues(); + const auto& rowPtr = matrix.vRowPtr(); + + // Row 0: cols 0,1 + EXPECT_EQ(colInd[0], 0); + EXPECT_EQ(colInd[1], 1); + EXPECT_DOUBLE_EQ(values[0], 1.0); + EXPECT_DOUBLE_EQ(values[1], 2.0); + + // Row 1: cols 0,1,2 + EXPECT_EQ(colInd[2], 0); + EXPECT_EQ(colInd[3], 1); + EXPECT_EQ(colInd[4], 2); + EXPECT_DOUBLE_EQ(values[2], 3.0); + EXPECT_DOUBLE_EQ(values[3], 5.0); + EXPECT_DOUBLE_EQ(values[4], 6.0); + + // Row 2: cols 0,1,2 + EXPECT_EQ(colInd[5], 0); + EXPECT_EQ(colInd[6], 1); + EXPECT_EQ(colInd[7], 2); + EXPECT_DOUBLE_EQ(values[5], 7.0); + EXPECT_DOUBLE_EQ(values[6], 8.0); + EXPECT_DOUBLE_EQ(values[7], 9.0); +} + +TEST(CSRBuilderTest, SparseWithGaps) +{ + cCSRMatrixBuilder builder; + + // Large matrix with only a few entries + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(50, 75, 2.0); + builder.AddEntry(99, 99, 3.0); + + auto matrix = builder.Build(0, 100, 100); + + EXPECT_EQ(matrix.iGlobalRowBegin(), 0); + EXPECT_EQ(matrix.iGlobalRowEnd(), 100); + EXPECT_EQ(matrix.vRowPtr().size(), 101u); + EXPECT_EQ(matrix.vColInd().size(), 3u); + EXPECT_EQ(matrix.vValues().size(), 3u); +} + +TEST(CSRBuilderTest, ScatteredDuplicates) +{ + cCSRMatrixBuilder builder; + + // Mix of unique and duplicate entries + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(0, 1, 2.0); + builder.AddEntry(0, 0, 1.0); // duplicate + builder.AddEntry(1, 1, 3.0); + builder.AddEntry(0, 1, 2.0); // duplicate + builder.AddEntry(1, 1, 3.0); // duplicate + builder.AddEntry(2, 2, 5.0); + + auto matrix = builder.Build(0, 3, 3); + + const auto& colInd = matrix.vColInd(); + const auto& values = matrix.vValues(); + + EXPECT_EQ(colInd.size(), 4u); // 4 unique positions: (0,0), (0,1), (1,1), (2,2) + + // (0,0) should be 2.0 + EXPECT_DOUBLE_EQ(values[0], 2.0); + // (0,1) should be 4.0 + EXPECT_DOUBLE_EQ(values[1], 4.0); + // (1,1) should be 6.0 + EXPECT_DOUBLE_EQ(values[2], 6.0); + // (2,2) should be 5.0 + EXPECT_DOUBLE_EQ(values[3], 5.0); +} + +TEST(CSRBuilderTest, SpMVWithBuiltMatrix) +{ + cCSRMatrixBuilder builder; + + // Build a matrix and verify SpMV works correctly + // [2 1] + // [1 3] + builder.AddEntry(0, 0, 2.0); + builder.AddEntry(0, 1, 1.0); + builder.AddEntry(1, 0, 1.0); + builder.AddEntry(1, 1, 3.0); + + auto matrix = builder.Build(0, 2, 2); + + std::vector x { 1.0, 2.0 }; + std::vector y; + + SpMV(matrix, x, y); + + ASSERT_EQ(y.size(), 2u); + EXPECT_DOUBLE_EQ(y[0], 4.0); // 2*1 + 1*2 + EXPECT_DOUBLE_EQ(y[1], 7.0); // 1*1 + 3*2 +} + +TEST(CSRBuilderTest, RowRangeSubset) +{ + cCSRMatrixBuilder builder; + + // Add entries for a larger matrix but build only a subset + builder.AddEntry(5, 0, 1.0); + builder.AddEntry(6, 1, 2.0); + builder.AddEntry(7, 2, 3.0); + builder.AddEntry(3, 0, 99.0); // Out of range + builder.AddEntry(10, 0, 99.0); // Out of range + + auto matrix = builder.Build(5, 8, 3); // Only rows 5-7 + + EXPECT_EQ(matrix.iGlobalRowBegin(), 5); + EXPECT_EQ(matrix.iGlobalRowEnd(), 8); + EXPECT_EQ(matrix.vColInd().size(), 3u); + EXPECT_EQ(matrix.vValues().size(), 3u); + + // Verify only in-range entries were included + const auto& values = matrix.vValues(); + EXPECT_DOUBLE_EQ(values[0], 1.0); + EXPECT_DOUBLE_EQ(values[1], 2.0); + EXPECT_DOUBLE_EQ(values[2], 3.0); +} + +TEST(CSRBuilderTest, AddEntriesBatch) +{ + cCSRMatrixBuilder builder; + + std::vector> triplets; + triplets.push_back(cTriplet(0, 0, 1.0)); + triplets.push_back(cTriplet(0, 1, 2.0)); + triplets.push_back(cTriplet(1, 1, 3.0)); + + builder.AddEntries(triplets); + + auto matrix = builder.Build(0, 2, 2); + + EXPECT_EQ(matrix.vColInd().size(), 3u); + EXPECT_EQ(matrix.vValues().size(), 3u); +} + +TEST(CSRBuilderTest, ClearAndReuse) +{ + cCSRMatrixBuilder builder; + + builder.AddEntry(0, 0, 1.0); + auto matrix1 = builder.Build(0, 1, 1); + EXPECT_EQ(matrix1.vValues().size(), 1u); + + builder.Clear(); + builder.AddEntry(0, 0, 2.0); + builder.AddEntry(0, 1, 3.0); + auto matrix2 = builder.Build(0, 1, 2); + EXPECT_EQ(matrix2.vValues().size(), 2u); +} + +// Test with different scalar types +TEST(CSRBuilderTest, FloatScalar) +{ + cCSRMatrixBuilder builder; + + builder.AddEntry(0, 0, 1.0f); + builder.AddEntry(0, 0, 2.0f); // Duplicate + + auto matrix = builder.Build(0, 1, 1); + + EXPECT_EQ(matrix.vValues().size(), 1u); + EXPECT_FLOAT_EQ(matrix.vValues()[0], 3.0f); +} + +TEST(CSRBuilderTest, ComplexScalar) +{ + using ComplexScalar = std::complex; + cCSRMatrixBuilder builder; + + builder.AddEntry(0, 0, ComplexScalar(1.0, 2.0)); + builder.AddEntry(0, 0, ComplexScalar(3.0, 4.0)); // Duplicate + + auto matrix = builder.Build(0, 1, 1); + + EXPECT_EQ(matrix.vValues().size(), 1u); + EXPECT_DOUBLE_EQ(matrix.vValues()[0].real(), 4.0); + EXPECT_DOUBLE_EQ(matrix.vValues()[0].imag(), 6.0); +} + +TEST(CSRBuilderTest, ConvenienceFunctionBuildFromTriplets) +{ + std::vector> triplets; + triplets.push_back(cTriplet(0, 0, 1.0)); + triplets.push_back(cTriplet(0, 1, 2.0)); + triplets.push_back(cTriplet(1, 1, 3.0)); + + auto matrix = BuildCSRFromTriplets(triplets, 0, 2, 2); + + EXPECT_EQ(matrix.vColInd().size(), 3u); + EXPECT_EQ(matrix.vValues().size(), 3u); + + const auto& values = matrix.vValues(); + EXPECT_DOUBLE_EQ(values[0], 1.0); + EXPECT_DOUBLE_EQ(values[1], 2.0); + EXPECT_DOUBLE_EQ(values[2], 3.0); +} + +TEST(CSRBuilderTest, ValidateCorrectMatrix) +{ + cCSRMatrixBuilder builder; + builder.AddEntry(0, 0, 1.0); + builder.AddEntry(0, 1, 2.0); + builder.AddEntry(1, 1, 3.0); + + auto matrix = builder.Build(0, 2, 2); + + std::string errorMsg; + EXPECT_TRUE(ValidateCSRMatrix(matrix, &errorMsg)) << errorMsg; +} + +TEST(CSRBuilderTest, ValidateDetectsInvalidRowPtrSize) +{ + // Manually create invalid matrix + std::vector rowPtr = { 0, 1 }; // Should be size 3 for 2 rows + std::vector colInd = { 0 }; + std::vector values = { 1.0 }; + + cCSRMatrix matrix(0, 2, 2, rowPtr, colInd, values); + + std::string errorMsg; + EXPECT_FALSE(ValidateCSRMatrix(matrix, &errorMsg)); + EXPECT_EQ(errorMsg, "Row pointer size mismatch"); +} + +TEST(CSRBuilderTest, ValidateDetectsUnsortedColumns) +{ + // Manually create matrix with unsorted columns + std::vector rowPtr = { 0, 2 }; + std::vector colInd = { 1, 0 }; // Unsorted! + std::vector values = { 1.0, 2.0 }; + + cCSRMatrix matrix(0, 1, 2, rowPtr, colInd, values); + + std::string errorMsg; + EXPECT_FALSE(ValidateCSRMatrix(matrix, &errorMsg)); + EXPECT_EQ(errorMsg, "Column indices are not sorted within row"); +} + +TEST(CSRBuilderTest, ValidateDetectsColumnOutOfBounds) +{ + // Manually create matrix with out-of-bounds column + std::vector rowPtr = { 0, 1 }; + std::vector colInd = { 10 }; // Out of bounds for colCount=2 + std::vector values = { 1.0 }; + + cCSRMatrix matrix(0, 1, 2, rowPtr, colInd, values); + + std::string errorMsg; + EXPECT_FALSE(ValidateCSRMatrix(matrix, &errorMsg)); + EXPECT_EQ(errorMsg, "Column index out of bounds"); +} + +TEST(CSRBuilderTest, BuiltMatrixAlwaysValid) +{ + // Any matrix built with the builder should be valid + cCSRMatrixBuilder builder; + + // Add random entries + builder.AddEntry(5, 3, 1.0); + builder.AddEntry(2, 7, 2.0); + builder.AddEntry(8, 1, 3.0); + builder.AddEntry(0, 9, 4.0); + builder.AddEntry(2, 7, 5.0); // duplicate + + auto matrix = builder.Build(0, 10, 10); + + EXPECT_TRUE(ValidateCSRMatrix(matrix)); +} diff --git a/tests/test_csrcomm_optimization.cpp b/tests/test_csrcomm_optimization.cpp new file mode 100644 index 0000000..3094fb3 --- /dev/null +++ b/tests/test_csrcomm_optimization.cpp @@ -0,0 +1,315 @@ +#include "CSRComm.h" +#include "CSRMatrix.h" +#include "CommPattern.h" +#include "Distribution.h" +#include +#include + +using namespace csr4mpi; + +// Test with double as the default scalar type +using Scalar = double; + +// Test that the binary search optimization correctly finds elements +TEST(CSRCommOptimizationTest, BinarySearchFindsFirstColumn) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with multiple columns per row + std::vector rowPtr = { 0, 5 }; // One row with 5 elements + std::vector colInd = { 0, 10, 20, 30, 40 }; // Sorted columns + std::vector values = { 1.0, 2.0, 3.0, 4.0, 5.0 }; + cCSRMatrix local(0, 1, 50, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); + + // Assemble to the first column + std::vector> entries; + cRemoteEntry e1; + e1.m_iGlobalRow = 0; + e1.m_iGlobalCol = 0; + e1.m_vValue = 10.0; + entries.push_back(e1); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 10.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[0], 11.0); // 1.0 + 10.0 + EXPECT_DOUBLE_EQ(vals[1], 2.0); // unchanged +} + +TEST(CSRCommOptimizationTest, BinarySearchFindsLastColumn) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with multiple columns per row + std::vector rowPtr = { 0, 5 }; + std::vector colInd = { 0, 10, 20, 30, 40 }; + std::vector values = { 1.0, 2.0, 3.0, 4.0, 5.0 }; + cCSRMatrix local(0, 1, 50, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); + + // Assemble to the last column + std::vector> entries; + cRemoteEntry e1; + e1.m_iGlobalRow = 0; + e1.m_iGlobalCol = 40; + e1.m_vValue = 20.0; + entries.push_back(e1); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 20.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[4], 25.0); // 5.0 + 20.0 + EXPECT_DOUBLE_EQ(vals[0], 1.0); // unchanged +} + +TEST(CSRCommOptimizationTest, BinarySearchFindsMiddleColumn) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with multiple columns per row + std::vector rowPtr = { 0, 5 }; + std::vector colInd = { 0, 10, 20, 30, 40 }; + std::vector values = { 1.0, 2.0, 3.0, 4.0, 5.0 }; + cCSRMatrix local(0, 1, 50, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); + + // Assemble to the middle column + std::vector> entries; + cRemoteEntry e1; + e1.m_iGlobalRow = 0; + e1.m_iGlobalCol = 20; + e1.m_vValue = 30.0; + entries.push_back(e1); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 30.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[2], 33.0); // 3.0 + 30.0 + EXPECT_DOUBLE_EQ(vals[0], 1.0); // unchanged + EXPECT_DOUBLE_EQ(vals[4], 5.0); // unchanged +} + +TEST(CSRCommOptimizationTest, BinarySearchMultipleRows) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with 3 rows, each with different number of columns + std::vector rowPtr = { 0, 2, 5, 7 }; + std::vector colInd = { 0, 5, 1, 3, 7, 2, 9 }; + std::vector values = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 }; + cCSRMatrix local(0, 3, 10, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(3, 1, 0); + + // Assemble to different rows + std::vector> entries; + + cRemoteEntry e1; + e1.m_iGlobalRow = 0; + e1.m_iGlobalCol = 5; + e1.m_vValue = 10.0; + + cRemoteEntry e2; + e2.m_iGlobalRow = 1; + e2.m_iGlobalCol = 3; + e2.m_vValue = 20.0; + + cRemoteEntry e3; + e3.m_iGlobalRow = 2; + e3.m_iGlobalCol = 9; + e3.m_vValue = 30.0; + + entries.push_back(e1); + entries.push_back(e2); + entries.push_back(e3); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 10.0, 20.0, 30.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[1], 12.0); // Row 0, col 5: 2.0 + 10.0 + EXPECT_DOUBLE_EQ(vals[3], 24.0); // Row 1, col 3: 4.0 + 20.0 + EXPECT_DOUBLE_EQ(vals[6], 37.0); // Row 2, col 9: 7.0 + 30.0 +} + +TEST(CSRCommOptimizationTest, LargeRowBinarySearchEfficiency) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Test with a large row to ensure binary search is used (O(log n) vs O(n)) + const int numCols = 1000; + std::vector rowPtr = { 0, numCols }; + std::vector colInd; + std::vector values; + + for (int i = 0; i < numCols; ++i) { + colInd.push_back(i * 10); // Sparse columns + values.push_back(static_cast(i + 1)); + } + + cCSRMatrix local(0, 1, numCols * 10, rowPtr, colInd, values); + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); + + // Add contribution to a column near the end + std::vector> entries; + cRemoteEntry e1; + e1.m_iGlobalRow = 0; + e1.m_iGlobalCol = 9990; // Last column + e1.m_vValue = 999.0; + entries.push_back(e1); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 999.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[999], 1999.0); // 1000.0 + 999.0 +} + +TEST(CSRCommOptimizationTest, MultipleDuplicateAccumulations) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with sorted columns + std::vector rowPtr = { 0, 3, 6 }; + std::vector colInd = { 0, 5, 10, 1, 6, 11 }; + std::vector values = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + cCSRMatrix local(0, 2, 15, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(2, 1, 0); + + // Multiple contributions to the same element + std::vector> entries; + for (int i = 0; i < 5; ++i) { + cRemoteEntry e; + e.m_iGlobalRow = 0; + e.m_iGlobalCol = 5; + e.m_vValue = 10.0; + entries.push_back(e); + } + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib(5, 10.0); + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[1], 52.0); // 2.0 + 5*10.0 +} + +TEST(CSRCommOptimizationTest, EmptyRow) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with an empty row + std::vector rowPtr = { 0, 0, 3 }; // Row 0 is empty + std::vector colInd = { 1, 5, 8 }; + std::vector values = { 1.0, 2.0, 3.0 }; + cCSRMatrix local(0, 2, 10, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(2, 1, 0); + + // Try to assemble to row 1 + std::vector> entries; + cRemoteEntry e1; + e1.m_iGlobalRow = 1; + e1.m_iGlobalCol = 5; + e1.m_vValue = 100.0; + entries.push_back(e1); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 100.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + EXPECT_DOUBLE_EQ(vals[1], 102.0); // 2.0 + 100.0 +} + +TEST(CSRCommOptimizationTest, NonExistentColumnIgnored) +{ + int initFlag = 0; + MPI_Initialized(&initFlag); + if (!initFlag) + MPI_Init(nullptr, nullptr); + + // Matrix with specific columns + std::vector rowPtr = { 0, 3 }; + std::vector colInd = { 0, 5, 10 }; + std::vector values = { 1.0, 2.0, 3.0 }; + cCSRMatrix local(0, 1, 15, rowPtr, colInd, values); + + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); + + // Try to assemble to a column that doesn't exist in the matrix + std::vector> entries; + cRemoteEntry e1; + e1.m_iGlobalRow = 0; + e1.m_iGlobalCol = 7; // Not in colInd + e1.m_vValue = 100.0; + entries.push_back(e1); + + cCommPattern pattern; + pattern.Build(entries, dist, 0, 1); + + std::vector contrib = { 100.0 }; + + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + + const auto& vals = local.vValues(); + // Values should remain unchanged + EXPECT_DOUBLE_EQ(vals[0], 1.0); + EXPECT_DOUBLE_EQ(vals[1], 2.0); + EXPECT_DOUBLE_EQ(vals[2], 3.0); +}