diff --git a/CMakeLists.txt b/CMakeLists.txt index 962946d..af71708 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,9 @@ cmake_minimum_required(VERSION 3.16) project(CSR4MPI LANGUAGES CXX) -set(CSR4MPI_VALUE_TYPE "1" CACHE STRING "Scalar type: 0=float,1=double,2=complex,3=complex") +# CSR4MPI_VALUE_TYPE is now deprecated - all types are compiled into the library. +# This option is kept for backward compatibility but has no effect. +set(CSR4MPI_VALUE_TYPE "1" CACHE STRING "DEPRECATED: All scalar types are now always compiled. This option has no effect.") option(CSR4MPI_ENABLE_BLAS "Enable optional BLAS placeholder path" OFF) option(CSR4MPI_ENABLE_OPENMP "Enable OpenMP parallel kernels" ON) diff --git a/bench/bench_large_spmv.cpp b/bench/bench_large_spmv.cpp index 5652abc..26f6649 100644 --- a/bench/bench_large_spmv.cpp +++ b/bench/bench_large_spmv.cpp @@ -9,8 +9,9 @@ #include using namespace csr4mpi; +using Scalar = double; -static cCSRMatrix ExtractLocal(const cCSRMatrix& full, const cRowDistribution& dist) +static cCSRMatrix ExtractLocal(const cCSRMatrix& full, const cRowDistribution& dist) { iIndex rBeg = dist.iGlobalRowBegin(); iIndex rEnd = dist.iGlobalRowEnd(); @@ -20,7 +21,7 @@ static cCSRMatrix ExtractLocal(const cCSRMatrix& full, const cRowDistribution& d const auto& FV = full.vValues(); std::vector lRP(localRows + 1, 0); std::vector lC; - std::vector lV; + std::vector lV; for (iSize lr = 0; lr < localRows; ++lr) { iIndex g = rBeg + lr; for (iIndex k = FRP[(size_t)g]; k < FRP[(size_t)g + 1]; ++k) { @@ -29,7 +30,7 @@ static cCSRMatrix ExtractLocal(const cCSRMatrix& full, const cRowDistribution& d } lRP[(size_t)lr + 1] = (iIndex)lC.size(); } - cCSRMatrix local(rBeg, rEnd, full.iGlobalColCount(), lRP, lC, lV, full.eSymmetry()); + cCSRMatrix local(rBeg, rEnd, full.iGlobalColCount(), lRP, lC, lV, full.eSymmetry()); local.AttachDistribution(std::make_shared(dist)); return local; } @@ -43,25 +44,25 @@ int main(int argc, char** argv) std::string matrixFile = (argc > 1 ? argv[1] : "bcsstk13.mtx"); std::string path = std::string(CSR4MPI_SOURCE_DIR) + "/tests/data/" + matrixFile; std::vector gRP, gCI; - std::vector gV; + std::vector gV; iIndex gRows = 0, gCols = 0; - if (!LoadMatrixMarket(path, gRP, gCI, gV, gRows, gCols)) { + if (!LoadMatrixMarket(path, gRP, gCI, gV, gRows, gCols)) { if (rank == 0) std::cout << "Matrix not found: " << path << "\n"; MPI_Finalize(); return 0; } auto dist = cRowDistribution::CreateBlockDistribution(gRows, worldSize, rank); - cCSRMatrix local = ExtractLocal(cCSRMatrix(0, gRows, gCols, gRP, gCI, gV), dist); + cCSRMatrix local = ExtractLocal(cCSRMatrix(0, gRows, gCols, gRP, gCI, gV), dist); // Build local slice of x (random deterministic) const auto& offs = dist.vRowOffsets(); iIndex cBeg = offs[(size_t)rank]; iIndex cEnd = offs[(size_t)rank + 1]; - std::vector xLocal; + std::vector xLocal; xLocal.reserve(static_cast(cEnd - cBeg)); for (iIndex g = cBeg; g < cEnd; ++g) - xLocal.push_back(static_cast(1)); - std::vector yLocal; + xLocal.push_back(static_cast(1)); + std::vector yLocal; int iters = (argc > 2 ? std::stoi(argv[2]) : 10); MPI_Barrier(MPI_COMM_WORLD); double t0 = MPI_Wtime(); diff --git a/bench/bench_operations.cpp b/bench/bench_operations.cpp index 6ec29f9..cdf1d1c 100644 --- a/bench/bench_operations.cpp +++ b/bench/bench_operations.cpp @@ -6,8 +6,9 @@ #include using namespace csr4mpi; +using Scalar = double; -static cCSRMatrix BuildRandomUniform(iSize rows, iSize cols, iSize nnzPerRow, unsigned seed) +static cCSRMatrix BuildRandomUniform(iSize rows, iSize cols, iSize nnzPerRow, unsigned seed) { std::mt19937_64 gen(seed); std::uniform_int_distribution colDist(0, cols - 1); @@ -15,7 +16,7 @@ static cCSRMatrix BuildRandomUniform(iSize rows, iSize cols, iSize nnzPerRow, un std::vector vRowPtr(rows + 1, 0); std::vector vColInd; vColInd.reserve(rows * nnzPerRow); - std::vector vValues; + std::vector vValues; vValues.reserve(rows * nnzPerRow); for (iSize r = 0; r < rows; ++r) { std::vector colsChosen; @@ -36,11 +37,11 @@ static cCSRMatrix BuildRandomUniform(iSize rows, iSize cols, iSize nnzPerRow, un for (auto c : colsChosen) { vColInd.push_back(static_cast(c)); double rv = valDist(gen); - vValues.push_back(static_cast(rv)); + vValues.push_back(static_cast(rv)); } vRowPtr[static_cast(r + 1)] = static_cast(vColInd.size()); } - return cCSRMatrix(0, rows, cols, vRowPtr, vColInd, vValues); + return cCSRMatrix(0, rows, cols, vRowPtr, vColInd, vValues); } int main(int argc, char** argv) @@ -62,12 +63,12 @@ int main(int argc, char** argv) repeats = std::stoi(argv[5]); auto A = BuildRandomUniform(rows, cols, nnzPerRow, 42); - std::vector x(static_cast(cols)); + std::vector x(static_cast(cols)); for (size_t i = 0; i < x.size(); ++i) - x[i] = static_cast(1); + x[i] = static_cast(1); // SpMV benchmark - std::vector y; + std::vector y; auto t0 = std::chrono::high_resolution_clock::now(); for (int r = 0; r < repeats; ++r) { SpMV(A, x, y); @@ -76,10 +77,10 @@ int main(int argc, char** argv) std::chrono::duration dtSpMV = t1 - t0; // SpMM benchmark - std::vector X(static_cast(cols * spmmCols)); + std::vector X(static_cast(cols * spmmCols)); for (size_t i = 0; i < X.size(); ++i) - X[i] = static_cast(1); - std::vector Y; + X[i] = static_cast(1); + std::vector Y; auto t2 = std::chrono::high_resolution_clock::now(); for (int r = 0; r < repeats; ++r) { SpMM(A, X, spmmCols, Y); diff --git a/include/BlasAdapter.h b/include/BlasAdapter.h index 2d84c5a..38a829d 100644 --- a/include/BlasAdapter.h +++ b/include/BlasAdapter.h @@ -1,9 +1,31 @@ #pragma once #include "CSRMatrix.h" +#include "Operations.h" #include "Global.h" #include +#include namespace csr4mpi { -bool bBlasEnabled(); -void SpMMBlas(const cCSRMatrix& A, const std::vector& X, iSize nCols, std::vector& Y); + +inline bool bBlasEnabled() +{ +#ifdef CSR4MPI_USE_BLAS + return true; +#else + return false; +#endif +} + +template +void SpMMBlas(const cCSRMatrix& A, const std::vector& X, iSize nCols, std::vector& Y) +{ + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); +#ifdef CSR4MPI_USE_BLAS + // Placeholder: in future, detect dense conversion threshold and call cblas_Xgemm. + // For now, just log once per call and fallback. + std::cerr << "[BLAS placeholder] Falling back to internal SpMM path.\n"; +#endif + SpMM(A, X, nCols, Y); +} + } diff --git a/include/CSRComm.h b/include/CSRComm.h index 0fbcd0c..60e5982 100644 --- a/include/CSRComm.h +++ b/include/CSRComm.h @@ -1,18 +1,163 @@ #pragma once #include "Global.h" +#include "CSRMatrix.h" +#include "CommPattern.h" #include #include namespace csr4mpi { -class cCSRMatrix; -class cCommPattern; +template class cCSRComm { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + public: - static void Assemble(cCSRMatrix& cLocal, - const cCommPattern& cPattern, - const std::vector& vValues, - MPI_Comm cComm); + static void Assemble(cCSRMatrix& cLocal, + const cCommPattern& cPattern, + const std::vector& vValues, + MPI_Comm cComm) + { + const std::vector& vSendRows = cPattern.vSendRows(); + const std::vector& vSendCols = cPattern.vSendCols(); + const std::vector& vSendOffsets = cPattern.vSendOffsets(); + const std::vector& vTargets = cPattern.vTargetRanks(); + + const std::vector& vRowPtr = cLocal.vRowPtr(); + const std::vector& vColInd = cLocal.vColInd(); + std::vector& vLocalVal = cLocal.vValues(); + + (void)vSendOffsets; + (void)vTargets; + + int iRank = 0; + int iWorldSize = 1; + MPI_Comm_rank(cComm, &iRank); + MPI_Comm_size(cComm, &iWorldSize); + + // Prepare send counts per target. + std::vector vSendCounts(iWorldSize, 0); + std::vector vRecvCounts(iWorldSize, 0); + + for (std::size_t i = 0; i < vTargets.size(); ++i) { + int iTarget = vTargets[i]; + iSize iBegin = vSendOffsets[i]; + iSize iEnd = vSendOffsets[i + 1]; + vSendCounts[static_cast(iTarget)] += static_cast(iEnd - iBegin); + } + + // Exchange counts. + MPI_Alltoall(vSendCounts.data(), 1, MPI_INT, + vRecvCounts.data(), 1, MPI_INT, + cComm); + + int iTotalSend = 0; + int iTotalRecv = 0; + for (int i = 0; i < iWorldSize; ++i) { + iTotalSend += vSendCounts[static_cast(i)]; + iTotalRecv += vRecvCounts[static_cast(i)]; + } + + struct cTriplet { + iIndex m_iRow; + iIndex m_iCol; + Scalar m_vVal; + }; + + std::vector vSendBuf(static_cast(iTotalSend)); + std::vector vRecvBuf(static_cast(iTotalRecv)); + + // Fill send buffer in rank-major order. + std::vector vOffsets(iWorldSize, 0); + for (int i = 1; i < iWorldSize; ++i) { + vOffsets[static_cast(i)] = vOffsets[static_cast(i - 1)] + vSendCounts[static_cast(i - 1)]; + } + + for (std::size_t iBucket = 0; iBucket < vTargets.size(); ++iBucket) { + int iTarget = vTargets[iBucket]; + iSize iBegin = vSendOffsets[iBucket]; + iSize iEnd = vSendOffsets[iBucket + 1]; + + for (iSize i = iBegin; i < iEnd; ++i) { + int iPos = vOffsets[static_cast(iTarget)]++; + vSendBuf[static_cast(iPos)].m_iRow = vSendRows[static_cast(i)]; + vSendBuf[static_cast(iPos)].m_iCol = vSendCols[static_cast(i)]; + vSendBuf[static_cast(iPos)].m_vVal = vValues[static_cast(i)]; + } + } + + // Build displacement arrays for Alltoallv. + std::vector vSendDispls(iWorldSize, 0); + std::vector vRecvDispls(iWorldSize, 0); + + for (int i = 1; i < iWorldSize; ++i) { + vSendDispls[static_cast(i)] = vSendDispls[static_cast(i - 1)] + vSendCounts[static_cast(i - 1)]; + vRecvDispls[static_cast(i)] = vRecvDispls[static_cast(i - 1)] + vRecvCounts[static_cast(i - 1)]; + } + + // Define MPI datatype for cTriplet with correct scalar handling. + MPI_Datatype tTripletType; + MPI_Datatype tCreatedScalarType = MPI_DATATYPE_NULL; + MPI_Datatype tScalarType = mpi_helper::GetMPIDatatype(&tCreatedScalarType); + + { + cTriplet cDummy; + int iBlockLengths[3] = { 1, 1, 1 }; + MPI_Aint vDispls[3]; + MPI_Aint iBase; + + MPI_Get_address(&cDummy, &iBase); + MPI_Get_address(&cDummy.m_iRow, &vDispls[0]); + MPI_Get_address(&cDummy.m_iCol, &vDispls[1]); + MPI_Get_address(&cDummy.m_vVal, &vDispls[2]); + + vDispls[0] -= iBase; + vDispls[1] -= iBase; + vDispls[2] -= iBase; + + MPI_Datatype vTypes[3] = { MPI_LONG_LONG, MPI_LONG_LONG, tScalarType }; + MPI_Type_create_struct(3, iBlockLengths, vDispls, vTypes, &tTripletType); + MPI_Type_commit(&tTripletType); + } + + MPI_Alltoallv(vSendBuf.data(), + vSendCounts.data(), + vSendDispls.data(), + tTripletType, + vRecvBuf.data(), + vRecvCounts.data(), + vRecvDispls.data(), + tTripletType, + cComm); + + MPI_Type_free(&tTripletType); + if (tCreatedScalarType != MPI_DATATYPE_NULL) { + MPI_Type_free(&tCreatedScalarType); + } + + // Accumulate received contributions into local CSR matrix. + for (const cTriplet& cT : vRecvBuf) { + iIndex iGlobalRow = cT.m_iRow; + iIndex iGlobalCol = cT.m_iCol; + + iIndex iLocalRow = static_cast(iGlobalRow - cLocal.iGlobalRowBegin()); + 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; + } + } + } + } }; + +// Type aliases for common scalar types +using cCSRCommF = cCSRComm; +using cCSRCommD = cCSRComm; +using cCSRCommCF = cCSRComm>; +using cCSRCommCD = cCSRComm>; + } diff --git a/include/CSRMatrix.h b/include/CSRMatrix.h index 3eabc3a..e9fcc56 100644 --- a/include/CSRMatrix.h +++ b/include/CSRMatrix.h @@ -6,28 +6,48 @@ namespace csr4mpi { class cRowDistribution; -} -namespace csr4mpi { +template class cCSRMatrix { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + public: - cCSRMatrix(); + using scalar_type = Scalar; + + cCSRMatrix() + : m_iGlobalRowBegin(0) + , m_iGlobalRowEnd(0) + , m_iGlobalColCount(0) + , m_pDistribution(nullptr) + { + } cCSRMatrix(iSize iGlobalRowBegin, iSize iGlobalRowEnd, iSize iGlobalColCount, std::vector vRowPtr, std::vector vColInd, - std::vector vValues, - eSymmStorage eSymm = eNone); + std::vector vValues, + eSymmStorage eSymm = eNone) + : m_iGlobalRowBegin(iGlobalRowBegin) + , m_iGlobalRowEnd(iGlobalRowEnd) + , m_iGlobalColCount(iGlobalColCount) + , m_vRowPtr(std::move(vRowPtr)) + , m_vColInd(std::move(vColInd)) + , m_vValues(std::move(vValues)) + , m_pDistribution(nullptr) + , m_eSymm(eSymm) + { + } - iSize iGlobalRowBegin() const; - iSize iGlobalRowEnd() const; + iSize iGlobalRowBegin() const { return m_iGlobalRowBegin; } + iSize iGlobalRowEnd() const { return m_iGlobalRowEnd; } - const std::vector& vRowPtr() const; - const std::vector& vColInd() const; - const std::vector& vValues() const; - iSize iGlobalColCount() const; + const std::vector& vRowPtr() const { return m_vRowPtr; } + const std::vector& vColInd() const { return m_vColInd; } + const std::vector& vValues() const { return m_vValues; } + std::vector& vValues() { return m_vValues; } + iSize iGlobalColCount() const { return m_iGlobalColCount; } eSymmStorage eSymmetry() const { return m_eSymm; } bool bIsSymmetric() const { return m_eSymm != eNone; } @@ -41,8 +61,15 @@ class cCSRMatrix { std::vector m_vRowPtr; std::vector m_vColInd; - std::vector m_vValues; + std::vector m_vValues; std::shared_ptr m_pDistribution; eSymmStorage m_eSymm { eNone }; }; + +// Type aliases for common scalar types +using cCSRMatrixF = cCSRMatrix; +using cCSRMatrixD = cCSRMatrix; +using cCSRMatrixCF = cCSRMatrix>; +using cCSRMatrixCD = cCSRMatrix>; + } diff --git a/include/CommPattern.h b/include/CommPattern.h index 7780f51..3febc45 100644 --- a/include/CommPattern.h +++ b/include/CommPattern.h @@ -1,32 +1,85 @@ #pragma once #include "Global.h" +#include "Distribution.h" #include +#include namespace csr4mpi { + +template struct cRemoteEntry { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); iIndex m_iGlobalRow; iIndex m_iGlobalCol; - vScalar m_vValue; + Scalar m_vValue; }; +// Type aliases for common scalar types +using cRemoteEntryF = cRemoteEntry; +using cRemoteEntryD = cRemoteEntry; +using cRemoteEntryCF = cRemoteEntry>; +using cRemoteEntryCD = cRemoteEntry>; + +template class cCommPattern { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + public: - cCommPattern(); + using scalar_type = Scalar; + + cCommPattern() = default; - void Build(const std::vector& vEntries, - const class cRowDistribution& cDistribution, + void Build(const std::vector>& vEntries, + const cRowDistribution& cDistribution, int iRank, - int iWorldSize); + int iWorldSize) + { + std::unordered_map>> mBuckets; + mBuckets.reserve(static_cast(iWorldSize)); + + for (const cRemoteEntry& cEntry : vEntries) { + int iOwner = cDistribution.iOwnerRank(cEntry.m_iGlobalRow); + mBuckets[iOwner].push_back(cEntry); + } + + m_vTargetRanks.clear(); + m_vSendRows.clear(); + m_vSendCols.clear(); + m_vSendOffsets.clear(); + + m_vSourceRanks.clear(); + m_vRecvRows.clear(); + m_vRecvCols.clear(); + m_vRecvOffsets.clear(); + + m_vSendOffsets.push_back(0); - const std::vector& vTargetRanks() const; - const std::vector& vSendRows() const; - const std::vector& vSendCols() const; - const std::vector& vSendOffsets() const; + for (const auto& cPair : mBuckets) { + int iTarget = cPair.first; + const std::vector>& vBucket = cPair.second; - const std::vector& vRecvRows() const; - const std::vector& vRecvCols() const; - const std::vector& vRecvOffsets() const; + m_vTargetRanks.push_back(iTarget); + + for (const cRemoteEntry& cEntry : vBucket) { + m_vSendRows.push_back(cEntry.m_iGlobalRow); + m_vSendCols.push_back(cEntry.m_iGlobalCol); + } + + m_vSendOffsets.push_back(static_cast(m_vSendRows.size())); + } + + m_vRecvOffsets.push_back(0); + } + + const std::vector& vTargetRanks() const { return m_vTargetRanks; } + const std::vector& vSendRows() const { return m_vSendRows; } + const std::vector& vSendCols() const { return m_vSendCols; } + const std::vector& vSendOffsets() const { return m_vSendOffsets; } + + const std::vector& vRecvRows() const { return m_vRecvRows; } + const std::vector& vRecvCols() const { return m_vRecvCols; } + const std::vector& vRecvOffsets() const { return m_vRecvOffsets; } private: std::vector m_vTargetRanks; @@ -39,4 +92,11 @@ class cCommPattern { std::vector m_vRecvCols; std::vector m_vRecvOffsets; }; + +// Type aliases for common scalar types +using cCommPatternF = cCommPattern; +using cCommPatternD = cCommPattern; +using cCommPatternCF = cCommPattern>; +using cCommPatternCD = cCommPattern>; + } diff --git a/include/DistributedOps.h b/include/DistributedOps.h index efede37..346941c 100644 --- a/include/DistributedOps.h +++ b/include/DistributedOps.h @@ -1,22 +1,63 @@ -/* - * @Author: error: error: git config user.name & please set dead value or install git && error: git config user.email & please set dead value or install git & please set dead value or install git - * @Date: 2025-11-30 10:52:25 - * @LastEditors: error: error: git config user.name & please set dead value or install git && error: git config user.email & please set dead value or install git & please set dead value or install git - * @LastEditTime: 2025-11-30 15:42:14 - * @FilePath: \CSR4MPI\include\DistributedOps.h - * @Description: - * - * Copyright (c) 2025 by error: error: git config user.name & please set dead value or install git && error: git config user.email & please set dead value or install git & please set dead value or install git, All Rights Reserved. - */ #pragma once #include "CSRMatrix.h" #include "Distribution.h" +#include "Operations.h" #include #include +#include +#include namespace csr4mpi { -void DistributedSpMV(const cCSRMatrix& A, const cRowDistribution& dist, - const std::vector& xLocalOrGlobal, - std::vector& y, - bool bReplicatedX, MPI_Comm comm); + +template +void DistributedSpMV(const cCSRMatrix& A, const cRowDistribution& dist, + const std::vector& xLocalOrGlobal, + std::vector& y, + bool bReplicatedX, MPI_Comm comm) +{ + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + + int worldSize = 1, worldRank = 0; + MPI_Comm_size(comm, &worldSize); + MPI_Comm_rank(comm, &worldRank); + + iSize globalCols = A.iGlobalColCount(); + iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); + + std::vector xGlobal; + if (bReplicatedX) { + if (xLocalOrGlobal.size() != static_cast(globalCols)) { + throw std::runtime_error("DistributedSpMV: replicated x size mismatch"); + } + xGlobal = xLocalOrGlobal; + } else { + iSize ownedBegin = dist.iGlobalRowBegin(); + iSize ownedEnd = dist.iGlobalRowEnd(); + iSize ownedCount = ownedEnd - ownedBegin; + if (xLocalOrGlobal.size() != static_cast(ownedCount)) { + throw std::runtime_error("DistributedSpMV: local x slice size mismatch"); + } + std::vector recvCounts(worldSize); + std::vector displs(worldSize); + if (worldSize == 1) { + recvCounts[0] = static_cast(ownedCount); + } else { + throw std::runtime_error("DistributedSpMV: recvCounts inference for multi-rank not implemented (need distribution metadata)"); + } + displs[0] = 0; + for (int r = 1; r < worldSize; ++r) + displs[r] = displs[r - 1] + recvCounts[r - 1]; + xGlobal.resize(static_cast(globalCols)); + MPI_Datatype tCreatedType = MPI_DATATYPE_NULL; + MPI_Datatype dt = mpi_helper::GetMPIDatatype(&tCreatedType); + MPI_Allgatherv(xLocalOrGlobal.data(), recvCounts[worldRank], dt, + xGlobal.data(), recvCounts.data(), displs.data(), dt, comm); + if (tCreatedType != MPI_DATATYPE_NULL) { + MPI_Type_free(&tCreatedType); + } + } + + SpMV(A, xGlobal, y); +} + } diff --git a/include/Global.h b/include/Global.h index c3a6cee..211cc42 100644 --- a/include/Global.h +++ b/include/Global.h @@ -1,6 +1,9 @@ #pragma once #include +#include +#include +#include namespace csr4mpi { enum eValueType { @@ -15,27 +18,76 @@ enum eSymmStorage { eSymmUpper, eSymmFull }; -} -#ifndef CSR4MPI_VALUE_TYPE -#define CSR4MPI_VALUE_TYPE 1 -#endif +using iIndex = std::int64_t; +using iSize = std::int64_t; -#include +// Type traits to detect complex types +template +struct is_complex : std::false_type {}; -namespace csr4mpi { -#if CSR4MPI_VALUE_TYPE == 0 -using vScalar = float; -#elif CSR4MPI_VALUE_TYPE == 1 -using vScalar = double; -#elif CSR4MPI_VALUE_TYPE == 2 -using vScalar = std::complex; -#elif CSR4MPI_VALUE_TYPE == 3 -using vScalar = std::complex; +template +struct is_complex> : std::true_type {}; + +template +inline constexpr bool is_complex_v = is_complex::value; + +// Get the real type underlying a scalar (identity for real, component type for complex) +template +struct real_type { + using type = T; +}; + +template +struct real_type> { + using type = T; +}; + +template +using real_type_t = typename real_type::type; + +// Check if a type is a supported scalar type +template +struct is_supported_scalar : std::bool_constant< + std::is_same_v || + std::is_same_v || + std::is_same_v> || + std::is_same_v>> {}; + +template +inline constexpr bool is_supported_scalar_v = is_supported_scalar::value; + +namespace mpi_helper { + // Helper to get MPI datatype for a scalar type + // pCreatedType is set if a custom MPI type was created (caller must free it) + template + inline MPI_Datatype GetMPIDatatype(MPI_Datatype* pCreatedType = nullptr) { + if constexpr (std::is_same_v) { + return MPI_FLOAT; + } else if constexpr (std::is_same_v) { + return MPI_DOUBLE; + } else if constexpr (std::is_same_v>) { +#ifdef MPI_C_FLOAT_COMPLEX + return MPI_C_FLOAT_COMPLEX; #else -#error "Unsupported CSR4MPI_VALUE_TYPE" + MPI_Datatype dt; + MPI_Type_contiguous(2, MPI_FLOAT, &dt); + MPI_Type_commit(&dt); + if (pCreatedType) *pCreatedType = dt; + return dt; #endif + } else if constexpr (std::is_same_v>) { +#ifdef MPI_C_DOUBLE_COMPLEX + return MPI_C_DOUBLE_COMPLEX; +#else + MPI_Datatype dt; + MPI_Type_contiguous(2, MPI_DOUBLE, &dt); + MPI_Type_commit(&dt); + if (pCreatedType) *pCreatedType = dt; + return dt; +#endif + } + } +} -using iIndex = std::int64_t; -using iSize = std::int64_t; } diff --git a/include/MatrixMarketLoader.h b/include/MatrixMarketLoader.h index c650171..ea52536 100644 --- a/include/MatrixMarketLoader.h +++ b/include/MatrixMarketLoader.h @@ -2,12 +2,136 @@ #include "Global.h" #include #include +#include +#include +#include +#include namespace csr4mpi { + +template bool LoadMatrixMarket(const std::string& sPath, std::vector& vRowPtr, std::vector& vColInd, - std::vector& vValues, + std::vector& vValues, iIndex& iRows, - iIndex& iCols); + iIndex& iCols) +{ + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + + 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 + { + return std::hash()((k.r << 32) ^ k.c); + } + }; + + std::ifstream fin(sPath); + if (!fin.is_open()) { + return false; + } + std::string line; + // Header + if (!std::getline(fin, line)) + return false; + if (line.rfind("%%MatrixMarket", 0) != 0) + return false; + // Expect tokens + bool bSymmetric = false; + bool bGeneral = false; + { + std::istringstream iss(line); + std::string mm, matrix, format, datatype, storage; + iss >> mm >> matrix >> format >> datatype >> storage; + if (matrix != "matrix" || format != "coordinate" || datatype != "real") { + return false; + } + if (storage == "general") + bGeneral = true; + else if (storage == "symmetric") + bSymmetric = true; + else + return false; + } + // Skip comments + while (std::getline(fin, line)) { + if (line.empty()) + continue; + if (line[0] == '%') + continue; + // First non-comment is size line + std::istringstream sz(line); + iIndex nnz; + sz >> iRows >> iCols >> nnz; + if (sz.fail()) + return false; + std::unordered_map mAcc; + mAcc.reserve(static_cast(nnz)); + for (iIndex k = 0; k < nnz; ++k) { + if (!std::getline(fin, line)) + return false; + if (line.empty()) { + --k; + continue; + } + if (line[0] == '%') { + --k; + continue; + } + std::istringstream es(line); + iIndex ir, ic; + double val; + es >> ir >> ic >> val; + if (es.fail()) + return false; + // 1-based to 0-based + cKey key { ir - 1, ic - 1 }; + auto it = mAcc.find(key); + Scalar v = static_cast(val); + if (it == mAcc.end()) + mAcc.emplace(key, v); + else + it->second += v; + if (bSymmetric && ir != ic) { + // 如果输入给的是仅单边三角（常见），不自动插入另一侧，在乘法阶段展开。 + } + } + // Move to vector and sort + struct cEntry { + iIndex r; + iIndex c; + Scalar v; + }; + std::vector entries; + entries.reserve(mAcc.size()); + for (const auto& p : mAcc) { + entries.push_back({ p.first.r, p.first.c, p.second }); + } + std::sort(entries.begin(), entries.end(), [](const cEntry& a, const cEntry& b) { + if (a.r != b.r) + return a.r < b.r; + return a.c < b.c; + }); + vRowPtr.assign(static_cast(iRows + 1), 0); + vColInd.clear(); + vValues.clear(); + for (const auto& e : entries) { + vRowPtr[static_cast(e.r + 1)]++; + vColInd.push_back(e.c); + vValues.push_back(e.v); + } + // prefix sum + for (iIndex r = 0; r < iRows; ++r) { + vRowPtr[static_cast(r + 1)] += vRowPtr[static_cast(r)]; + } + return true; + } + return false; +} + } diff --git a/include/MumpsAdapter.h b/include/MumpsAdapter.h index c8af903..9a3f4a7 100644 --- a/include/MumpsAdapter.h +++ b/include/MumpsAdapter.h @@ -1,28 +1,91 @@ #pragma once #include "Global.h" +#include "CSRMatrix.h" +#include "Distribution.h" #include namespace csr4mpi { -class cCSRMatrix; -class cRowDistribution; +template class cMumpsAdapter { + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + public: // Export local owned CSR block to MUMPS COO arrays (1-based indices). // Minimizes reallocations by reserving exact capacity before appending. - static void ExportLocalBlock(const cCSRMatrix& cLocal, + static void ExportLocalBlock(const cCSRMatrix& cLocal, const cRowDistribution& cDistribution, std::vector& vIRN, std::vector& vJCN, - std::vector& vA); + std::vector& vA) + { + const std::vector& vRowPtr = cLocal.vRowPtr(); + const std::vector& vColInd = cLocal.vColInd(); + const std::vector& vValues = cLocal.vValues(); + + const iSize iRowBegin = cDistribution.iGlobalRowBegin(); + const iSize iRowEnd = cDistribution.iGlobalRowEnd(); + const iSize iLocalRows = iRowEnd - iRowBegin; + const iSize iLocalNNZ = static_cast(vRowPtr[static_cast(iLocalRows)]); + + // Minimize reallocations: reserve exact additional capacity for this append. + vIRN.reserve(vIRN.size() + static_cast(iLocalNNZ)); + vJCN.reserve(vJCN.size() + static_cast(iLocalNNZ)); + vA.reserve(vA.size() + static_cast(iLocalNNZ)); + + for (iSize iGlobalRow = iRowBegin; iGlobalRow < iRowEnd; ++iGlobalRow) { + const iIndex iLocalRow = static_cast(iGlobalRow - iRowBegin); + const iIndex iStart = vRowPtr[static_cast(iLocalRow)]; + const iIndex iEnd = vRowPtr[static_cast(iLocalRow + 1)]; + + for (iIndex i = iStart; i < iEnd; ++i) { + const iIndex iCol = vColInd[static_cast(i)]; + vIRN.push_back(static_cast(iGlobalRow + 1)); + vJCN.push_back(static_cast(iCol + 1)); + vA.push_back(vValues[static_cast(i)]); + } + } + } // Export into pre-allocated buffers provided by the caller to avoid extra allocations. // Returns the number of nonzeros written. Buffers must have capacity >= local nnz. - static iSize ExportLocalBlockInto(const cCSRMatrix& cLocal, + static iSize ExportLocalBlockInto(const cCSRMatrix& cLocal, const cRowDistribution& cDistribution, int* pIRN, int* pJCN, - vScalar* pA); + Scalar* pA) + { + const std::vector& vRowPtr = cLocal.vRowPtr(); + const std::vector& vColInd = cLocal.vColInd(); + const std::vector& vValues = cLocal.vValues(); + + const iSize iRowBegin = cDistribution.iGlobalRowBegin(); + const iSize iRowEnd = cDistribution.iGlobalRowEnd(); + const iSize iLocalRows = iRowEnd - iRowBegin; + const iSize iLocalNNZ = static_cast(vRowPtr[static_cast(iLocalRows)]); + + iSize k = 0; + for (iSize iGlobalRow = iRowBegin; iGlobalRow < iRowEnd; ++iGlobalRow) { + const iIndex iLocalRow = static_cast(iGlobalRow - iRowBegin); + const iIndex iStart = vRowPtr[static_cast(iLocalRow)]; + const iIndex iEnd = vRowPtr[static_cast(iLocalRow + 1)]; + for (iIndex i = iStart; i < iEnd; ++i, ++k) { + const iIndex iCol = vColInd[static_cast(i)]; + pIRN[k] = static_cast(iGlobalRow + 1); + pJCN[k] = static_cast(iCol + 1); + pA[k] = vValues[static_cast(i)]; + } + } + + return iLocalNNZ; + } }; + +// Type aliases for common scalar types +using cMumpsAdapterF = cMumpsAdapter; +using cMumpsAdapterD = cMumpsAdapter; +using cMumpsAdapterCF = cMumpsAdapter>; +using cMumpsAdapterCD = cMumpsAdapter>; + } diff --git a/include/Operations.h b/include/Operations.h index f05cde5..6df1569 100644 --- a/include/Operations.h +++ b/include/Operations.h @@ -1,11 +1,446 @@ #pragma once #include "CSRMatrix.h" +#include "Distribution.h" +#include #include +#include +#include +#ifdef CSR4MPI_USE_OPENMP +#include +#endif namespace csr4mpi { -void SpMV(const cCSRMatrix& A, const std::vector& x, std::vector& y); -void SpMV(const cCSRMatrix& A, std::vector& x); -void SpMM(const cCSRMatrix& A, const std::vector& X, iSize nCols, std::vector& Y); -void SpMMInPlace(const cCSRMatrix& A, std::vector& X, iSize nCols); + +template +void SpMV(const cCSRMatrix& A, const std::vector& x, std::vector& y) +{ + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + + iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); + const auto& rowPtr = A.vRowPtr(); + const auto& colInd = A.vColInd(); + const auto& values = A.vValues(); + eSymmStorage eSymm = A.eSymmetry(); + auto pDist = A.pDistribution(); + int rank = 0, worldSize = 1; + if (pDist) { + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &worldSize); + } + // 构造全局 x (必要时聚合) + std::vector xGlobalBuf; + const std::vector* pXFull = &x; + if (x.size() != static_cast(A.iGlobalColCount())) { + if (!pDist) + throw std::runtime_error("SpMV: input vector length mismatch and no distribution"); + if (worldSize < 2) + throw std::runtime_error("SpMV: distribution mismatch in single-rank build"); + const auto& offsets = pDist->vRowOffsets(); + if (offsets.size() != static_cast(worldSize + 1)) + throw std::runtime_error("SpMV: offsets size mismatch"); + if (A.iGlobalColCount() != pDist->iGlobalRowCount()) + throw std::runtime_error("SpMV: non-square distributed matrix unsupported"); + iSize expected = offsets[static_cast(rank + 1)] - offsets[static_cast(rank)]; + if (x.size() != static_cast(expected)) + throw std::runtime_error("SpMV: local slice size mismatch"); + xGlobalBuf.resize(static_cast(A.iGlobalColCount())); + std::vector counts(worldSize), displs(worldSize); + for (int r = 0; r < worldSize; ++r) + counts[r] = static_cast(offsets[static_cast(r + 1)] - offsets[static_cast(r)]); + displs[0] = 0; + for (int r = 1; r < worldSize; ++r) + displs[r] = displs[r - 1] + counts[r - 1]; + MPI_Datatype tCreatedType = MPI_DATATYPE_NULL; + MPI_Datatype dt = mpi_helper::GetMPIDatatype(&tCreatedType); + MPI_Allgatherv(x.data(), counts[rank], dt, xGlobalBuf.data(), counts.data(), displs.data(), dt, MPI_COMM_WORLD); + if (tCreatedType != MPI_DATATYPE_NULL) { + MPI_Type_free(&tCreatedType); + } + pXFull = &xGlobalBuf; + } + const std::vector& xFull = *pXFull; + y.assign(static_cast(localRows), static_cast(0)); + iSize rowBase = A.iGlobalRowBegin(); + bool symTriangular = (eSymm == eSymmLower || eSymm == eSymmUpper); + if (!symTriangular || eSymm == eSymmFull) { +#ifdef CSR4MPI_USE_OPENMP +#pragma omp parallel for schedule(static) +#endif + for (iSize r = 0; r < localRows; ++r) { + iIndex start = rowPtr[static_cast(r)]; + iIndex end = rowPtr[static_cast(r + 1)]; + Scalar sum = static_cast(0); + for (iIndex k = start; k < end; ++k) { + iIndex c = colInd[static_cast(k)]; + sum += values[static_cast(k)] * xFull[static_cast(c)]; + } + y[static_cast(r)] = sum; + } + return; + } + // Triangular symmetric expansion + bool hasRemote = (pDist && worldSize > 1); + std::vector sendRows; + std::vector sendVals; + std::vector sendCounts(hasRemote ? worldSize : 1, 0); + // 主循环：并行仅在无远程通信时启用 + if (!hasRemote) { +#ifdef CSR4MPI_USE_OPENMP +#pragma omp parallel for schedule(static) +#endif + for (iSize r = 0; r < localRows; ++r) { + iIndex start = rowPtr[static_cast(r)]; + iIndex end = rowPtr[static_cast(r + 1)]; + Scalar diagAcc = static_cast(0); + iIndex gRow = r + rowBase; + struct cUpd { + iIndex idx; + Scalar val; + }; + std::vector vAdds; + vAdds.reserve(static_cast(end - start)); + for (iIndex k = start; k < end; ++k) { + iIndex c = colInd[static_cast(k)]; + Scalar a = values[static_cast(k)]; + if (c == gRow) { + diagAcc += a * xFull[static_cast(c)]; + continue; + } + bool lowerCond = (gRow >= c); + if ((eSymm == eSymmLower && lowerCond) || (eSymm == eSymmUpper && !lowerCond)) { + diagAcc += a * xFull[static_cast(c)]; + iIndex maybeLocal = c - rowBase; + if (maybeLocal >= 0 && maybeLocal < localRows) + vAdds.push_back({ maybeLocal, a * xFull[static_cast(gRow)] }); + } + } + if (!vAdds.empty()) { +#ifdef CSR4MPI_USE_OPENMP +#pragma omp critical +#endif + for (const auto& u : vAdds) + y[static_cast(u.idx)] += u.val; + } + y[static_cast(r)] += diagAcc; + } + return; + } + // 有远程通信时，顺序执行收集远程对称展开贡献 + for (iSize r = 0; r < localRows; ++r) { + iIndex start = rowPtr[static_cast(r)]; + iIndex end = rowPtr[static_cast(r + 1)]; + Scalar diagAcc = static_cast(0); + iIndex gRow = r + rowBase; + for (iIndex k = start; k < end; ++k) { + iIndex c = colInd[static_cast(k)]; + Scalar a = values[static_cast(k)]; + if (c == gRow) { + diagAcc += a * xFull[static_cast(c)]; + continue; + } + bool lowerCond = (gRow >= c); + if ((eSymm == eSymmLower && lowerCond) || (eSymm == eSymmUpper && !lowerCond)) { + diagAcc += a * xFull[static_cast(c)]; + int owner = pDist->iOwnerRank(c); + if (owner == rank) { + iIndex maybeLocal = c - rowBase; + if (maybeLocal >= 0 && maybeLocal < localRows) + y[static_cast(maybeLocal)] += a * xFull[static_cast(gRow)]; + } else { + sendRows.push_back(c); + sendVals.push_back(a * xFull[static_cast(gRow)]); + } + } + } + y[static_cast(r)] += diagAcc; + } + if (worldSize > 1) { + // 统计发送计数（允许 0） + for (size_t i = 0; i < sendRows.size(); ++i) { + int tgt = pDist->iOwnerRank(sendRows[i]); + sendCounts[tgt]++; + } + std::vector recvCounts(worldSize, 0); + MPI_Alltoall(sendCounts.data(), 1, MPI_INT, recvCounts.data(), 1, MPI_INT, MPI_COMM_WORLD); + int totalSend = 0, totalRecv = 0; + for (int i = 0; i < worldSize; ++i) { + totalSend += sendCounts[i]; + totalRecv += recvCounts[i]; + } + struct cPair { + iIndex row; + Scalar val; + }; + std::vector sendBuf(static_cast(totalSend)); + std::vector recvBuf(static_cast(totalRecv)); + std::vector sendDispls(worldSize, 0), recvDispls(worldSize, 0), offsets(worldSize, 0); + for (int i = 1; i < worldSize; ++i) { + sendDispls[i] = sendDispls[i - 1] + sendCounts[i - 1]; + recvDispls[i] = recvDispls[i - 1] + recvCounts[i - 1]; + } + for (size_t i = 0; i < sendRows.size(); ++i) { + int tgt = pDist->iOwnerRank(sendRows[i]); + int pos = sendDispls[tgt] + offsets[tgt]++; + sendBuf[static_cast(pos)].row = sendRows[i]; + sendBuf[static_cast(pos)].val = sendVals[i]; + } + MPI_Datatype tPairType; + MPI_Datatype tCreatedScalarType = MPI_DATATYPE_NULL; + { + cPair dummy; + int bl[2] = { 1, 1 }; + MPI_Aint disp[2]; + MPI_Aint base; + MPI_Get_address(&dummy, &base); + MPI_Get_address(&dummy.row, &disp[0]); + MPI_Get_address(&dummy.val, &disp[1]); + disp[0] -= base; + disp[1] -= base; + MPI_Datatype types[2]; + types[0] = MPI_LONG_LONG; + types[1] = mpi_helper::GetMPIDatatype(&tCreatedScalarType); + MPI_Type_create_struct(2, bl, disp, types, &tPairType); + MPI_Type_commit(&tPairType); + } + MPI_Alltoallv(sendBuf.data(), sendCounts.data(), sendDispls.data(), tPairType, + recvBuf.data(), recvCounts.data(), recvDispls.data(), tPairType, + MPI_COMM_WORLD); + MPI_Type_free(&tPairType); + if (tCreatedScalarType != MPI_DATATYPE_NULL) { + MPI_Type_free(&tCreatedScalarType); + } + iSize rowBase2 = A.iGlobalRowBegin(); + for (const cPair& pr : recvBuf) { + iIndex local = pr.row - rowBase2; + if (local >= 0 && local < localRows) + y[static_cast(local)] += pr.val; + } + } +} + +template +void SpMV(const cCSRMatrix& A, std::vector& x) +{ + std::vector y; + SpMV(A, x, y); + x.swap(y); +} + +template +void SpMM(const cCSRMatrix& A, const std::vector& X, iSize nCols, std::vector& Y) +{ + static_assert(is_supported_scalar_v, "Scalar must be float, double, std::complex, or std::complex"); + + iSize globalCols = A.iGlobalColCount(); + iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); + if (X.size() != static_cast(globalCols * nCols)) { + throw std::runtime_error("SpMM: input dense matrix size mismatch"); + } + Y.assign(static_cast(localRows * nCols), static_cast(0)); + const auto& rowPtr = A.vRowPtr(); + const auto& colInd = A.vColInd(); + const auto& values = A.vValues(); + eSymmStorage eSymm = A.eSymmetry(); + iSize rowBase = A.iGlobalRowBegin(); + auto pDist = A.pDistribution(); + int rank = 0, worldSize = 1; + if (pDist) { + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &worldSize); + } + if (eSymm == eNone || eSymm == eSymmFull) { +#ifdef CSR4MPI_USE_OPENMP +#pragma omp parallel for schedule(static) +#endif + 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 c = colInd[static_cast(k)]; + Scalar a = values[static_cast(k)]; + const Scalar* xColBase = &X[static_cast(c)]; + Scalar* yRowBase = &Y[static_cast(r)]; + for (iSize j = 0; j < nCols; ++j) { + yRowBase[static_cast(j * localRows)] += a * xColBase[static_cast(j * globalCols)]; + } + } + } + } else { + bool isLowerStored = (eSymm == eSymmLower); + bool isUpperStored = (eSymm == eSymmUpper); + bool hasRemote = (pDist && worldSize > 1); + // 如果无远程通信：批量本地更新，按行收集再一次性写入 + if (!hasRemote) { +#ifdef CSR4MPI_USE_OPENMP +#pragma omp parallel for schedule(static) +#endif + for (iSize r = 0; r < localRows; ++r) { + iIndex start = rowPtr[static_cast(r)]; + iIndex end = rowPtr[static_cast(r + 1)]; + iIndex gRow = r + rowBase; + struct cUpd { + iIndex idx; + std::vector vals; + }; + std::vector vAdds; + for (iIndex k = start; k < end; ++k) { + iIndex c = colInd[static_cast(k)]; + Scalar a = values[static_cast(k)]; + bool lowerCond = (gRow >= c); + if (c == gRow) { + const Scalar* xColBase = &X[static_cast(c)]; + Scalar* yRowBase = &Y[static_cast(r)]; + for (iSize j = 0; j < nCols; ++j) + yRowBase[static_cast(j * localRows)] += a * xColBase[static_cast(j * globalCols)]; + continue; + } + if ((isLowerStored && lowerCond) || (isUpperStored && !lowerCond)) { + const Scalar* xColBase_c = &X[static_cast(c)]; + Scalar* yRowBase_r = &Y[static_cast(r)]; + for (iSize j = 0; j < nCols; ++j) + yRowBase_r[static_cast(j * localRows)] += a * xColBase_c[static_cast(j * globalCols)]; + iIndex maybeLocal = c - rowBase; + if (maybeLocal >= 0 && maybeLocal < localRows) { + const Scalar* xColBase_r = &X[static_cast(gRow)]; + cUpd upd; + upd.idx = maybeLocal; + upd.vals.resize(static_cast(nCols)); + for (iSize j = 0; j < nCols; ++j) + upd.vals[static_cast(j)] = a * xColBase_r[static_cast(j * globalCols)]; + vAdds.push_back(std::move(upd)); + } + } + } + if (!vAdds.empty()) { +#ifdef CSR4MPI_USE_OPENMP +#pragma omp critical +#endif + for (const auto& u : vAdds) { + Scalar* yRowBase_c = &Y[static_cast(u.idx)]; + for (iSize j = 0; j < nCols; ++j) + yRowBase_c[static_cast(j * localRows)] += u.vals[static_cast(j)]; + } + } + } + return; + } + // 远程路径：行块聚合，针对每个远程镜像行收集长度为 nCols 的向量，减少元数据 + std::unordered_map> mRemoteBlocks; + for (iSize r = 0; r < localRows; ++r) { + iIndex start = rowPtr[static_cast(r)]; + iIndex end = rowPtr[static_cast(r + 1)]; + iIndex gRow = r + rowBase; + for (iIndex k = start; k < end; ++k) { + iIndex c = colInd[static_cast(k)]; + Scalar a = values[static_cast(k)]; + bool lowerCond = (gRow >= c); + if (c == gRow) { + const Scalar* xColBase = &X[static_cast(c)]; + Scalar* yRowBase = &Y[static_cast(r)]; + for (iSize j = 0; j < nCols; ++j) + yRowBase[static_cast(j * localRows)] += a * xColBase[static_cast(j * globalCols)]; + continue; + } + if ((isLowerStored && lowerCond) || (isUpperStored && !lowerCond)) { + const Scalar* xColBase_c = &X[static_cast(c)]; + Scalar* yRowBase_r = &Y[static_cast(r)]; + for (iSize j = 0; j < nCols; ++j) + yRowBase_r[static_cast(j * localRows)] += a * xColBase_c[static_cast(j * globalCols)]; + int owner = pDist->iOwnerRank(c); + if (owner == rank) { + iIndex maybeLocal = c - rowBase; + if (maybeLocal >= 0 && maybeLocal < localRows) { + const Scalar* xColBase_r = &X[static_cast(gRow)]; + Scalar* yRowBase_c = &Y[static_cast(maybeLocal)]; + for (iSize j = 0; j < nCols; ++j) + yRowBase_c[static_cast(j * localRows)] += a * xColBase_r[static_cast(j * globalCols)]; + } + } else { + const Scalar* xColBase_r = &X[static_cast(gRow)]; + auto& vecRef = mRemoteBlocks[c]; + if (vecRef.empty()) + vecRef.assign(static_cast(nCols), static_cast(0)); + for (iSize j = 0; j < nCols; ++j) + vecRef[static_cast(j)] += a * xColBase_r[static_cast(j * globalCols)]; + } + } + } + } + // 统计远程行块数量 + std::vector sendRowCounts(worldSize, 0); + for (const auto& kv : mRemoteBlocks) { + int tgt = pDist->iOwnerRank(kv.first); + sendRowCounts[tgt]++; + } + std::vector recvCounts(worldSize, 0); + MPI_Alltoall(sendRowCounts.data(), 1, MPI_INT, recvCounts.data(), 1, MPI_INT, MPI_COMM_WORLD); + int totalSend = 0, totalRecv = 0; + for (int i = 0; i < worldSize; ++i) { + totalSend += sendRowCounts[i]; + totalRecv += recvCounts[i]; + } + std::vector sendRowsList; + sendRowsList.reserve(static_cast(totalSend)); + std::vector sendValsList; + sendValsList.reserve(static_cast(totalSend) * static_cast(nCols)); + std::vector> perRankRows(worldSize); + std::vector> perRankVals(worldSize); + for (auto& kv : mRemoteBlocks) { + int tgt = pDist->iOwnerRank(kv.first); + perRankRows[tgt].push_back(kv.first); + auto& vec = kv.second; + perRankVals[tgt].insert(perRankVals[tgt].end(), vec.begin(), vec.end()); + } + for (int rnk = 0; rnk < worldSize; ++rnk) { + sendRowsList.insert(sendRowsList.end(), perRankRows[rnk].begin(), perRankRows[rnk].end()); + sendValsList.insert(sendValsList.end(), perRankVals[rnk].begin(), perRankVals[rnk].end()); + } + std::vector sendRowDispls(worldSize, 0), recvRowDispls(worldSize, 0); + for (int i = 1; i < worldSize; ++i) { + sendRowDispls[i] = sendRowDispls[i - 1] + sendRowCounts[i - 1]; + recvRowDispls[i] = recvRowDispls[i - 1] + recvCounts[i - 1]; + } + std::vector recvRowsList(static_cast(totalRecv)); + MPI_Alltoallv(sendRowsList.data(), sendRowCounts.data(), sendRowDispls.data(), MPI_LONG_LONG, + recvRowsList.data(), recvCounts.data(), recvRowDispls.data(), MPI_LONG_LONG, MPI_COMM_WORLD); + std::vector sendValCounts(worldSize, 0), recvValCounts(worldSize, 0), sendValDispls(worldSize, 0), recvValDispls(worldSize, 0); + for (int i = 0; i < worldSize; ++i) { + sendValCounts[i] = sendRowCounts[i] * static_cast(nCols); + recvValCounts[i] = recvCounts[i] * static_cast(nCols); + } + for (int i = 1; i < worldSize; ++i) { + sendValDispls[i] = sendValDispls[i - 1] + sendValCounts[i - 1]; + recvValDispls[i] = recvValDispls[i - 1] + recvValCounts[i - 1]; + } + std::vector recvValsList(static_cast(totalRecv) * static_cast(nCols)); + MPI_Datatype tCreatedScalarType = MPI_DATATYPE_NULL; + MPI_Datatype scalarType = mpi_helper::GetMPIDatatype(&tCreatedScalarType); + MPI_Alltoallv(sendValsList.data(), sendValCounts.data(), sendValDispls.data(), scalarType, + recvValsList.data(), recvValCounts.data(), recvValDispls.data(), scalarType, + MPI_COMM_WORLD); + for (size_t blk = 0; blk < recvRowsList.size(); ++blk) { + iIndex local = recvRowsList[blk] - rowBase; + if (local >= 0 && local < localRows) { + Scalar* yRowBase = &Y[static_cast(local)]; + const Scalar* vals = &recvValsList[blk * static_cast(nCols)]; + for (iSize j = 0; j < nCols; ++j) { + yRowBase[static_cast(j * localRows)] += vals[static_cast(j)]; + } + } + } + if (tCreatedScalarType != MPI_DATATYPE_NULL) { + MPI_Type_free(&tCreatedScalarType); + } + } +} + +template +void SpMMInPlace(const cCSRMatrix& A, std::vector& X, iSize nCols) +{ + std::vector Y; + SpMM(A, X, nCols, Y); + X.swap(Y); +} + } diff --git a/src/BlasAdapter.cpp b/src/BlasAdapter.cpp index c98ca81..a42d97b 100644 --- a/src/BlasAdapter.cpp +++ b/src/BlasAdapter.cpp @@ -1,29 +1,7 @@ -#include "BlasAdapter.h" -#include "Operations.h" -#include - -namespace csr4mpi { +// BlasAdapter is now a header-only implementation. +// This file is kept for build system compatibility but contains no implementation. +// See BlasAdapter.h for the complete implementation. -bool bBlasEnabled() -{ -#ifdef CSR4MPI_USE_BLAS - return true; -#else - return false; -#endif -} - -void SpMMBlas(const cCSRMatrix& A, const std::vector& X, iSize nCols, std::vector& Y) -{ -#ifdef CSR4MPI_USE_BLAS - // Placeholder: in future, detect dense conversion threshold and call cblas_Xgemm. - // For now, just log once per call and fallback. - // Avoid performance penalty of i/o in real scenarios; kept minimal here. - // (Could be guarded by an environment variable.) - // std::cerr used intentionally to separate from normal output. - std::cerr << "[BLAS placeholder] Falling back to internal SpMM path.\n"; -#endif - SpMM(A, X, nCols, Y); -} +#include "BlasAdapter.h" -} +// Explicit template instantiations are not needed since the functions are header-only. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 686e402..af8ce1b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,12 +7,13 @@ set(SRC_LIST Operations.cpp BlasAdapter.cpp MatrixMarketLoader.cpp + DistributedOps.cpp ) add_library(csr4mpi STATIC ${SRC_LIST}) target_include_directories(csr4mpi PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../include) -target_compile_definitions(csr4mpi PUBLIC CSR4MPI_VALUE_TYPE=${CSR4MPI_VALUE_TYPE}) +# CSR4MPI_VALUE_TYPE is no longer needed - all types are templated and compiled into headers if(CSR4MPI_ENABLE_BLAS) target_compile_definitions(csr4mpi PUBLIC CSR4MPI_USE_BLAS) endif() diff --git a/src/CSRComm.cpp b/src/CSRComm.cpp index 764811f..b7a6a43 100644 --- a/src/CSRComm.cpp +++ b/src/CSRComm.cpp @@ -1,178 +1,7 @@ -#include "CSRComm.h" -#include "CSRMatrix.h" -#include "CommPattern.h" - -#include - -namespace csr4mpi { -void cCSRComm::Assemble(cCSRMatrix& cLocal, - const cCommPattern& cPattern, - const std::vector& vValues, - MPI_Comm cComm) -{ - const std::vector& vSendRows = cPattern.vSendRows(); - const std::vector& vSendCols = cPattern.vSendCols(); - const std::vector& vSendOffsets = cPattern.vSendOffsets(); - const std::vector& vTargets = cPattern.vTargetRanks(); - - const std::vector& vRowPtr = cLocal.vRowPtr(); - const std::vector& vColInd = cLocal.vColInd(); - std::vector& vLocalVal = const_cast&>(cLocal.vValues()); - - (void)vSendOffsets; - (void)vTargets; - - int iRank = 0; - int iWorldSize = 1; - MPI_Comm_rank(cComm, &iRank); - MPI_Comm_size(cComm, &iWorldSize); - - // Prepare send counts per target. - std::vector vSendCounts(iWorldSize, 0); - std::vector vRecvCounts(iWorldSize, 0); - - for (std::size_t i = 0; i < vTargets.size(); ++i) { - int iTarget = vTargets[i]; - iSize iBegin = vSendOffsets[i]; - iSize iEnd = vSendOffsets[i + 1]; - vSendCounts[static_cast(iTarget)] += static_cast(iEnd - iBegin); - } - - // Exchange counts. - MPI_Alltoall(vSendCounts.data(), 1, MPI_INT, - vRecvCounts.data(), 1, MPI_INT, - cComm); - - int iTotalSend = 0; - int iTotalRecv = 0; - for (int i = 0; i < iWorldSize; ++i) { - iTotalSend += vSendCounts[static_cast(i)]; - iTotalRecv += vRecvCounts[static_cast(i)]; - } - - struct cTriplet { - iIndex m_iRow; - iIndex m_iCol; - vScalar m_vVal; - }; - - std::vector vSendBuf(static_cast(iTotalSend)); - std::vector vRecvBuf(static_cast(iTotalRecv)); - - // Fill send buffer in rank-major order. - std::vector vOffsets(iWorldSize, 0); - for (int i = 1; i < iWorldSize; ++i) { - vOffsets[static_cast(i)] = vOffsets[static_cast(i - 1)] + vSendCounts[static_cast(i - 1)]; - } - - for (std::size_t iBucket = 0; iBucket < vTargets.size(); ++iBucket) { - int iTarget = vTargets[iBucket]; - iSize iBegin = vSendOffsets[iBucket]; - iSize iEnd = vSendOffsets[iBucket + 1]; +// CSRComm is now a header-only template class. +// This file is kept for build system compatibility but contains no implementation. +// See CSRComm.h for the complete implementation. - for (iSize i = iBegin; i < iEnd; ++i) { - int iPos = vOffsets[static_cast(iTarget)]++; - vSendBuf[static_cast(iPos)].m_iRow = vSendRows[static_cast(i)]; - vSendBuf[static_cast(iPos)].m_iCol = vSendCols[static_cast(i)]; - vSendBuf[static_cast(iPos)].m_vVal = vValues[static_cast(i)]; - } - } - - // Build displacement arrays for Alltoallv. - std::vector vSendDispls(iWorldSize, 0); - std::vector vRecvDispls(iWorldSize, 0); - - for (int i = 1; i < iWorldSize; ++i) { - vSendDispls[static_cast(i)] = vSendDispls[static_cast(i - 1)] + vSendCounts[static_cast(i - 1)]; - vRecvDispls[static_cast(i)] = vRecvDispls[static_cast(i - 1)] + vRecvCounts[static_cast(i - 1)]; - } - - // Define MPI datatype for cTriplet with correct scalar handling. - MPI_Datatype tTripletType; - MPI_Datatype tScalarType; - bool bCreatedScalarType = false; - { - // Map vScalar -> MPI_Datatype -#if CSR4MPI_VALUE_TYPE == 0 - tScalarType = MPI_FLOAT; -#elif CSR4MPI_VALUE_TYPE == 1 - tScalarType = MPI_DOUBLE; -#elif CSR4MPI_VALUE_TYPE == 2 -#ifdef MPI_C_FLOAT_COMPLEX - tScalarType = MPI_C_FLOAT_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_FLOAT, &tScalarType); - MPI_Type_commit(&tScalarType); - bCreatedScalarType = true; -#endif -#elif CSR4MPI_VALUE_TYPE == 3 -#ifdef MPI_C_DOUBLE_COMPLEX - tScalarType = MPI_C_DOUBLE_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_DOUBLE, &tScalarType); - MPI_Type_commit(&tScalarType); - bCreatedScalarType = true; -#endif -#else -#error "Unsupported CSR4MPI_VALUE_TYPE in CSRComm.cpp" -#endif - - cTriplet cDummy; - int iBlockLengths[3] = { 1, 1, 1 }; - MPI_Aint vDispls[3]; - MPI_Aint iBase; - - MPI_Get_address(&cDummy, &iBase); - MPI_Get_address(&cDummy.m_iRow, &vDispls[0]); - MPI_Get_address(&cDummy.m_iCol, &vDispls[1]); - MPI_Get_address(&cDummy.m_vVal, &vDispls[2]); - - vDispls[0] -= iBase; - vDispls[1] -= iBase; - vDispls[2] -= iBase; - - MPI_Datatype vTypes[3] = { MPI_LONG_LONG, MPI_LONG_LONG, tScalarType }; - MPI_Type_create_struct(3, iBlockLengths, vDispls, vTypes, &tTripletType); - MPI_Type_commit(&tTripletType); - } - - MPI_Alltoallv(vSendBuf.data(), - vSendCounts.data(), - vSendDispls.data(), - tTripletType, - vRecvBuf.data(), - vRecvCounts.data(), - vRecvDispls.data(), - tTripletType, - cComm); - - MPI_Type_free(&tTripletType); - if (bCreatedScalarType) { - MPI_Type_free(&tScalarType); - } - - // Accumulate received contributions into local CSR matrix. - for (const cTriplet& cT : vRecvBuf) { - iIndex iGlobalRow = cT.m_iRow; - iIndex iGlobalCol = cT.m_iCol; - - // Map global row to local row using cRowDistribution on the caller side - // when building vValues and pattern. - // Here we assume that global row belongs to this rank and local row is - // simply offset from cLocal.iGlobalRowBegin(). - - iIndex iLocalRow = static_cast(iGlobalRow - cLocal.iGlobalRowBegin()); - 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; - } - } - } +#include "CSRComm.h" - -} -} +// Explicit template instantiations are not needed since the class is header-only. diff --git a/src/CSRMatrix.cpp b/src/CSRMatrix.cpp index 2f38544..69f1d4a 100644 --- a/src/CSRMatrix.cpp +++ b/src/CSRMatrix.cpp @@ -1,59 +1,7 @@ -#include "CSRMatrix.h" - -namespace csr4mpi { -cCSRMatrix::cCSRMatrix() - : m_iGlobalRowBegin(0) - , m_iGlobalRowEnd(0) - , m_iGlobalColCount(0) - , m_pDistribution(nullptr) -{ -} - -cCSRMatrix::cCSRMatrix(iSize iGlobalRowBegin, - iSize iGlobalRowEnd, - iSize iGlobalColCount, - std::vector vRowPtr, - std::vector vColInd, - std::vector vValues, - eSymmStorage eSymm) - : m_iGlobalRowBegin(iGlobalRowBegin) - , m_iGlobalRowEnd(iGlobalRowEnd) - , m_iGlobalColCount(iGlobalColCount) - , m_vRowPtr(std::move(vRowPtr)) - , m_vColInd(std::move(vColInd)) - , m_vValues(std::move(vValues)) - , m_pDistribution(nullptr) - , m_eSymm(eSymm) -{ -} - -iSize cCSRMatrix::iGlobalRowBegin() const -{ - return m_iGlobalRowBegin; -} +// CSRMatrix is now a header-only template class. +// This file is kept for build system compatibility but contains no implementation. +// See CSRMatrix.h for the complete implementation. -iSize cCSRMatrix::iGlobalRowEnd() const -{ - return m_iGlobalRowEnd; -} - -const std::vector& cCSRMatrix::vRowPtr() const -{ - return m_vRowPtr; -} - -const std::vector& cCSRMatrix::vColInd() const -{ - return m_vColInd; -} - -const std::vector& cCSRMatrix::vValues() const -{ - return m_vValues; -} +#include "CSRMatrix.h" -iSize cCSRMatrix::iGlobalColCount() const -{ - return m_iGlobalColCount; -} -} +// Explicit template instantiations are not needed since the class is header-only. diff --git a/src/CommPattern.cpp b/src/CommPattern.cpp index e2aabab..870ba52 100644 --- a/src/CommPattern.cpp +++ b/src/CommPattern.cpp @@ -1,87 +1,7 @@ -#include "CommPattern.h" -#include "Distribution.h" -#include - -namespace csr4mpi { -cCommPattern::cCommPattern() = default; - -void cCommPattern::Build(const std::vector& vEntries, - const cRowDistribution& cDistribution, - int iRank, - int iWorldSize) -{ - std::unordered_map> mBuckets; - mBuckets.reserve(static_cast(iWorldSize)); - - for (const cRemoteEntry& cEntry : vEntries) { - int iOwner = cDistribution.iOwnerRank(cEntry.m_iGlobalRow); - // Include local and remote entries; local entries will be handled without MPI. - mBuckets[iOwner].push_back(cEntry); - } - - m_vTargetRanks.clear(); - m_vSendRows.clear(); - m_vSendCols.clear(); - m_vSendOffsets.clear(); - - m_vSourceRanks.clear(); - m_vRecvRows.clear(); - m_vRecvCols.clear(); - m_vRecvOffsets.clear(); - - m_vSendOffsets.push_back(0); - - for (const auto& cPair : mBuckets) { - int iTarget = cPair.first; - const std::vector& vBucket = cPair.second; - - m_vTargetRanks.push_back(iTarget); +// CommPattern is now a header-only template class. +// This file is kept for build system compatibility but contains no implementation. +// See CommPattern.h for the complete implementation. - for (const cRemoteEntry& cEntry : vBucket) { - m_vSendRows.push_back(cEntry.m_iGlobalRow); - m_vSendCols.push_back(cEntry.m_iGlobalCol); - } - - m_vSendOffsets.push_back(static_cast(m_vSendRows.size())); - } - - // Receive side will be filled during MPI setup inside communication layer. - // Here we only ensure empty but valid prefix array. - m_vRecvOffsets.push_back(0); -} - -const std::vector& cCommPattern::vTargetRanks() const -{ - return m_vTargetRanks; -} - -const std::vector& cCommPattern::vSendRows() const -{ - return m_vSendRows; -} - -const std::vector& cCommPattern::vSendCols() const -{ - return m_vSendCols; -} - -const std::vector& cCommPattern::vSendOffsets() const -{ - return m_vSendOffsets; -} - -const std::vector& cCommPattern::vRecvRows() const -{ - return m_vRecvRows; -} - -const std::vector& cCommPattern::vRecvCols() const -{ - return m_vRecvCols; -} +#include "CommPattern.h" -const std::vector& cCommPattern::vRecvOffsets() const -{ - return m_vRecvOffsets; -} -} +// Explicit template instantiations are not needed since the class is header-only. diff --git a/src/DistributedOps.cpp b/src/DistributedOps.cpp index 707ffdc..0014c4a 100644 --- a/src/DistributedOps.cpp +++ b/src/DistributedOps.cpp @@ -1,64 +1,7 @@ -#include "DistributedOps.h" -#include "Operations.h" -#include -#include -#include - -namespace csr4mpi { -void DistributedSpMV(const cCSRMatrix& A, const cRowDistribution& dist, - const std::vector& xLocalOrGlobal, - std::vector& y, - bool bReplicatedX, MPI_Comm comm) -{ - int worldSize = 1, worldRank = 0; - MPI_Comm_size(comm, &worldSize); - MPI_Comm_rank(comm, &worldRank); +// DistributedOps is now a header-only template implementation. +// This file is kept for build system compatibility but contains no implementation. +// See DistributedOps.h for the complete implementation. - iSize globalCols = A.iGlobalColCount(); - iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); - - std::vector xGlobal; - if (bReplicatedX) { - if (xLocalOrGlobal.size() != static_cast(globalCols)) { - throw std::runtime_error("DistributedSpMV: replicated x size mismatch"); - } - // Use directly - xGlobal = xLocalOrGlobal; // copy to allow unified code path - } else { - // Assume column distribution mirrors row distribution (square or compatible). - // Each rank provides its owned segment; gather all. - iSize ownedBegin = dist.iGlobalRowBegin(); // local process row begin - iSize ownedEnd = dist.iGlobalRowEnd(); - iSize ownedCount = ownedEnd - ownedBegin; - if (xLocalOrGlobal.size() != static_cast(ownedCount)) { - throw std::runtime_error("DistributedSpMV: local x slice size mismatch"); - } - std::vector recvCounts(worldSize); - std::vector displs(worldSize); - // Reconstruct counts by assuming uniform block distribution from constructor parameters not stored. - // Fallback: infer per-rank sizes from ownedCount when size==1; when size>1 require external provision (not implemented yet). - if (worldSize == 1) { - recvCounts[0] = static_cast(ownedCount); - } else { - throw std::runtime_error("DistributedSpMV: recvCounts inference for multi-rank not implemented (need distribution metadata)"); - } - displs[0] = 0; - for (int r = 1; r < worldSize; ++r) - displs[r] = displs[r - 1] + recvCounts[r - 1]; - xGlobal.resize(static_cast(globalCols)); - // 选择 MPI_Datatype（仅支持实数 float/double；复杂与其它需扩展） - MPI_Datatype dt; - if constexpr (std::is_same_v) - dt = MPI_FLOAT; - else if constexpr (std::is_same_v) - dt = MPI_DOUBLE; - else { - throw std::runtime_error("DistributedSpMV: MPI gather not yet implemented for complex types"); - } - MPI_Allgatherv(xLocalOrGlobal.data(), recvCounts[worldRank], dt, - xGlobal.data(), recvCounts.data(), displs.data(), dt, comm); - } +#include "DistributedOps.h" - SpMV(A, xGlobal, y); -} -} +// Explicit template instantiations are not needed since the function is header-only. diff --git a/src/MatrixMarketLoader.cpp b/src/MatrixMarketLoader.cpp index c961ed1..0dd2195 100644 --- a/src/MatrixMarketLoader.cpp +++ b/src/MatrixMarketLoader.cpp @@ -1,132 +1,7 @@ -#include "MatrixMarketLoader.h" -#include -#include -#include -#include -#include +// MatrixMarketLoader is now a header-only template implementation. +// This file is kept for build system compatibility but contains no implementation. +// See MatrixMarketLoader.h for the complete implementation. -namespace csr4mpi { -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 - { - return std::hash()((k.r << 32) ^ k.c); - } -}; +#include "MatrixMarketLoader.h" -bool LoadMatrixMarket(const std::string& sPath, - std::vector& vRowPtr, - std::vector& vColInd, - std::vector& vValues, - iIndex& iRows, - iIndex& iCols) -{ - std::ifstream fin(sPath); - if (!fin.is_open()) { - return false; - } - std::string line; - // Header - if (!std::getline(fin, line)) - return false; - if (line.rfind("%%MatrixMarket", 0) != 0) - return false; - // Expect tokens - bool bSymmetric = false; - bool bGeneral = false; - { - std::istringstream iss(line); - std::string mm, matrix, format, datatype, storage; - iss >> mm >> matrix >> format >> datatype >> storage; - if (matrix != "matrix" || format != "coordinate" || datatype != "real") { - return false; // 限制到最简单情况 - } - if (storage == "general") - bGeneral = true; - else if (storage == "symmetric") - bSymmetric = true; - else - return false; - } - // Skip comments - while (std::getline(fin, line)) { - if (line.empty()) - continue; - if (line[0] == '%') - continue; - // First non-comment is size line - std::istringstream sz(line); - iIndex nnz; - sz >> iRows >> iCols >> nnz; - if (sz.fail()) - return false; - std::unordered_map mAcc; - mAcc.reserve(static_cast(nnz)); - for (iIndex k = 0; k < nnz; ++k) { - if (!std::getline(fin, line)) - return false; - if (line.empty()) { - --k; - continue; - } - if (line[0] == '%') { - --k; - continue; - } - std::istringstream es(line); - iIndex ir, ic; - double val; - es >> ir >> ic >> val; // real value - if (es.fail()) - return false; - // 1-based to 0-based - cKey key { ir - 1, ic - 1 }; - auto it = mAcc.find(key); - vScalar v = static_cast(val); - if (it == mAcc.end()) - mAcc.emplace(key, v); - else - it->second += v; - if (bSymmetric && ir != ic) { - // 如果输入给的是仅单边三角（常见），不自动插入另一侧，在乘法阶段展开。 - // 若输入已经包含另一侧则此处将自然形成独立 key。 - } - } - // Move to vector and sort - struct cEntry { - iIndex r; - iIndex c; - vScalar v; - }; - std::vector entries; - entries.reserve(mAcc.size()); - for (const auto& p : mAcc) { - entries.push_back({ p.first.r, p.first.c, p.second }); - } - std::sort(entries.begin(), entries.end(), [](const cEntry& a, const cEntry& b) { - if (a.r != b.r) - return a.r < b.r; - return a.c < b.c; - }); - vRowPtr.assign(static_cast(iRows + 1), 0); - vColInd.clear(); - vValues.clear(); - for (const auto& e : entries) { - vRowPtr[static_cast(e.r + 1)]++; - vColInd.push_back(e.c); - vValues.push_back(e.v); - } - // prefix sum - for (iIndex r = 0; r < iRows; ++r) { - vRowPtr[static_cast(r + 1)] += vRowPtr[static_cast(r)]; - } - // 暂不返回对称类型枚举；调用侧可二次判定。此处只装载原始条目。 - return true; - } - return false; -} -} +// Explicit template instantiations are not needed since the function is header-only. diff --git a/src/MumpsAdapter.cpp b/src/MumpsAdapter.cpp index ceecf60..2875c77 100644 --- a/src/MumpsAdapter.cpp +++ b/src/MumpsAdapter.cpp @@ -1,70 +1,7 @@ -#include "MumpsAdapter.h" -#include "CSRMatrix.h" -#include "Distribution.h" - -namespace csr4mpi { -void cMumpsAdapter::ExportLocalBlock(const cCSRMatrix& cLocal, - const cRowDistribution& cDistribution, - std::vector& vIRN, - std::vector& vJCN, - std::vector& vA) -{ - const std::vector& vRowPtr = cLocal.vRowPtr(); - const std::vector& vColInd = cLocal.vColInd(); - const std::vector& vValues = cLocal.vValues(); - - const iSize iRowBegin = cDistribution.iGlobalRowBegin(); - const iSize iRowEnd = cDistribution.iGlobalRowEnd(); - const iSize iLocalRows = iRowEnd - iRowBegin; - const iSize iLocalNNZ = static_cast(vRowPtr[static_cast(iLocalRows)]); - - // Minimize reallocations: reserve exact additional capacity for this append. - vIRN.reserve(vIRN.size() + static_cast(iLocalNNZ)); - vJCN.reserve(vJCN.size() + static_cast(iLocalNNZ)); - vA.reserve(vA.size() + static_cast(iLocalNNZ)); - - for (iSize iGlobalRow = iRowBegin; iGlobalRow < iRowEnd; ++iGlobalRow) { - const iIndex iLocalRow = static_cast(iGlobalRow - iRowBegin); - const iIndex iStart = vRowPtr[static_cast(iLocalRow)]; - const iIndex iEnd = vRowPtr[static_cast(iLocalRow + 1)]; +// MumpsAdapter is now a header-only template class. +// This file is kept for build system compatibility but contains no implementation. +// See MumpsAdapter.h for the complete implementation. - for (iIndex i = iStart; i < iEnd; ++i) { - const iIndex iCol = vColInd[static_cast(i)]; - vIRN.push_back(static_cast(iGlobalRow + 1)); - vJCN.push_back(static_cast(iCol + 1)); - vA.push_back(vValues[static_cast(i)]); - } - } -} - -iSize cMumpsAdapter::ExportLocalBlockInto(const cCSRMatrix& cLocal, - const cRowDistribution& cDistribution, - int* pIRN, - int* pJCN, - vScalar* pA) -{ - const std::vector& vRowPtr = cLocal.vRowPtr(); - const std::vector& vColInd = cLocal.vColInd(); - const std::vector& vValues = cLocal.vValues(); - - const iSize iRowBegin = cDistribution.iGlobalRowBegin(); - const iSize iRowEnd = cDistribution.iGlobalRowEnd(); - const iSize iLocalRows = iRowEnd - iRowBegin; - const iSize iLocalNNZ = static_cast(vRowPtr[static_cast(iLocalRows)]); - - iSize k = 0; - for (iSize iGlobalRow = iRowBegin; iGlobalRow < iRowEnd; ++iGlobalRow) { - const iIndex iLocalRow = static_cast(iGlobalRow - iRowBegin); - const iIndex iStart = vRowPtr[static_cast(iLocalRow)]; - const iIndex iEnd = vRowPtr[static_cast(iLocalRow + 1)]; - for (iIndex i = iStart; i < iEnd; ++i, ++k) { - const iIndex iCol = vColInd[static_cast(i)]; - pIRN[k] = static_cast(iGlobalRow + 1); - pJCN[k] = static_cast(iCol + 1); - pA[k] = vValues[static_cast(i)]; - } - } +#include "MumpsAdapter.h" - return iLocalNNZ; -} -} +// Explicit template instantiations are not needed since the class is header-only. diff --git a/src/Operations.cpp b/src/Operations.cpp index b01d5d6..5695c6d 100644 --- a/src/Operations.cpp +++ b/src/Operations.cpp @@ -1,524 +1,7 @@ -/* - * @Author: error: error: git config user.name & please set dead value or install git && error: git config user.email & please set dead value or install git & please set dead value or install git - * @Date: 2025-11-30 04:30:36 - * @LastEditors: error: error: git config user.name & please set dead value or install git && error: git config user.email & please set dead value or install git & please set dead value or install git - * @LastEditTime: 2025-11-30 14:13:39 - * @FilePath: \CSR4MPI\src\Operations.cpp - * @Description: - * - * Copyright (c) 2025 by error: error: git config user.name & please set dead value or install git && error: git config user.email & please set dead value or install git & please set dead value or install git, All Rights Reserved. - */ -#include "Operations.h" -#include "Distribution.h" -#include -#include -#ifdef CSR4MPI_USE_OPENMP -#include -#endif -#include - -namespace csr4mpi { -void SpMV(const cCSRMatrix& A, const std::vector& x, std::vector& y) -{ - iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); - const auto& rowPtr = A.vRowPtr(); - const auto& colInd = A.vColInd(); - const auto& values = A.vValues(); - eSymmStorage eSymm = A.eSymmetry(); - auto pDist = A.pDistribution(); - int rank = 0, worldSize = 1; - if (pDist) { - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &worldSize); - } - // 构造全局 x (必要时聚合) - std::vector xGlobalBuf; - const std::vector* pXFull = &x; - if (x.size() != static_cast(A.iGlobalColCount())) { - if (!pDist) - throw std::runtime_error("SpMV: input vector length mismatch and no distribution"); - if (worldSize < 2) - throw std::runtime_error("SpMV: distribution mismatch in single-rank build"); - const auto& offsets = pDist->vRowOffsets(); - if (offsets.size() != static_cast(worldSize + 1)) - throw std::runtime_error("SpMV: offsets size mismatch"); - if (A.iGlobalColCount() != pDist->iGlobalRowCount()) - throw std::runtime_error("SpMV: non-square distributed matrix unsupported"); - iSize expected = offsets[static_cast(rank + 1)] - offsets[static_cast(rank)]; - if (x.size() != static_cast(expected)) - throw std::runtime_error("SpMV: local slice size mismatch"); - xGlobalBuf.resize(static_cast(A.iGlobalColCount())); - std::vector counts(worldSize), displs(worldSize); - for (int r = 0; r < worldSize; ++r) - counts[r] = static_cast(offsets[static_cast(r + 1)] - offsets[static_cast(r)]); - displs[0] = 0; - for (int r = 1; r < worldSize; ++r) - displs[r] = displs[r - 1] + counts[r - 1]; - MPI_Datatype dt; - if constexpr (std::is_same_v) - dt = MPI_FLOAT; - else if constexpr (std::is_same_v) - dt = MPI_DOUBLE; - else if constexpr (std::is_same_v>) { -#ifdef MPI_C_FLOAT_COMPLEX - dt = MPI_C_FLOAT_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_FLOAT, &dt); - MPI_Type_commit(&dt); -#endif - } else if constexpr (std::is_same_v>) { -#ifdef MPI_C_DOUBLE_COMPLEX - dt = MPI_C_DOUBLE_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_DOUBLE, &dt); - MPI_Type_commit(&dt); -#endif - } else { - throw std::runtime_error("SpMV: unsupported scalar type"); - } - MPI_Allgatherv(x.data(), counts[rank], dt, xGlobalBuf.data(), counts.data(), displs.data(), dt, MPI_COMM_WORLD); -#if !defined(MPI_C_FLOAT_COMPLEX) - if constexpr (std::is_same_v>) { - if (dt != MPI_C_FLOAT_COMPLEX) - MPI_Type_free(&dt); - } -#endif -#if !defined(MPI_C_DOUBLE_COMPLEX) - if constexpr (std::is_same_v>) { - if (dt != MPI_C_DOUBLE_COMPLEX) - MPI_Type_free(&dt); - } -#endif - pXFull = &xGlobalBuf; - } - const std::vector& xFull = *pXFull; - y.assign(static_cast(localRows), static_cast(0)); - iSize rowBase = A.iGlobalRowBegin(); - bool symTriangular = (eSymm == eSymmLower || eSymm == eSymmUpper); - if (!symTriangular || eSymm == eSymmFull) { -#ifdef CSR4MPI_USE_OPENMP -#pragma omp parallel for schedule(static) -#endif - for (iSize r = 0; r < localRows; ++r) { - iIndex start = rowPtr[static_cast(r)]; - iIndex end = rowPtr[static_cast(r + 1)]; - vScalar sum = static_cast(0); - for (iIndex k = start; k < end; ++k) { - iIndex c = colInd[static_cast(k)]; - sum += values[static_cast(k)] * xFull[static_cast(c)]; - } - y[static_cast(r)] = sum; - } - return; - } - // Triangular symmetric expansion - bool hasRemote = (pDist && worldSize > 1); - std::vector sendRows; - std::vector sendVals; - std::vector sendCounts(hasRemote ? worldSize : 1, 0); - // 主循环：并行仅在无远程通信时启用 - if (!hasRemote) { -#ifdef CSR4MPI_USE_OPENMP -#pragma omp parallel for schedule(static) -#endif - for (iSize r = 0; r < localRows; ++r) { - iIndex start = rowPtr[static_cast(r)]; - iIndex end = rowPtr[static_cast(r + 1)]; - vScalar diagAcc = static_cast(0); - iIndex gRow = r + rowBase; - struct cUpd { - iIndex idx; - vScalar val; - }; - std::vector vAdds; - vAdds.reserve(static_cast(end - start)); - for (iIndex k = start; k < end; ++k) { - iIndex c = colInd[static_cast(k)]; - vScalar a = values[static_cast(k)]; - if (c == gRow) { - diagAcc += a * xFull[static_cast(c)]; - continue; - } - bool lowerCond = (gRow >= c); - if ((eSymm == eSymmLower && lowerCond) || (eSymm == eSymmUpper && !lowerCond)) { - diagAcc += a * xFull[static_cast(c)]; - iIndex maybeLocal = c - rowBase; - if (maybeLocal >= 0 && maybeLocal < localRows) - vAdds.push_back({ maybeLocal, a * xFull[static_cast(gRow)] }); - } - } - if (!vAdds.empty()) { -#ifdef CSR4MPI_USE_OPENMP -#pragma omp critical -#endif - for (const auto& u : vAdds) - y[static_cast(u.idx)] += u.val; - } - y[static_cast(r)] += diagAcc; - } - return; - } - // 有远程通信时，顺序执行收集远程对称展开贡献 - for (iSize r = 0; r < localRows; ++r) { - iIndex start = rowPtr[static_cast(r)]; - iIndex end = rowPtr[static_cast(r + 1)]; - vScalar diagAcc = static_cast(0); - iIndex gRow = r + rowBase; - for (iIndex k = start; k < end; ++k) { - iIndex c = colInd[static_cast(k)]; - vScalar a = values[static_cast(k)]; - if (c == gRow) { - diagAcc += a * xFull[static_cast(c)]; - continue; - } - bool lowerCond = (gRow >= c); - if ((eSymm == eSymmLower && lowerCond) || (eSymm == eSymmUpper && !lowerCond)) { - diagAcc += a * xFull[static_cast(c)]; - int owner = pDist->iOwnerRank(c); - if (owner == rank) { - iIndex maybeLocal = c - rowBase; - if (maybeLocal >= 0 && maybeLocal < localRows) - y[static_cast(maybeLocal)] += a * xFull[static_cast(gRow)]; - } else { - sendRows.push_back(c); - sendVals.push_back(a * xFull[static_cast(gRow)]); - } - } - } - y[static_cast(r)] += diagAcc; - } - if (worldSize > 1) { - // 统计发送计数（允许 0） - for (size_t i = 0; i < sendRows.size(); ++i) { - int tgt = pDist->iOwnerRank(sendRows[i]); - sendCounts[tgt]++; - } - std::vector recvCounts(worldSize, 0); - MPI_Alltoall(sendCounts.data(), 1, MPI_INT, recvCounts.data(), 1, MPI_INT, MPI_COMM_WORLD); - int totalSend = 0, totalRecv = 0; - for (int i = 0; i < worldSize; ++i) { - totalSend += sendCounts[i]; - totalRecv += recvCounts[i]; - } - struct cPair { - iIndex row; - vScalar val; - }; - std::vector sendBuf(static_cast(totalSend)); - std::vector recvBuf(static_cast(totalRecv)); - std::vector sendDispls(worldSize, 0), recvDispls(worldSize, 0), offsets(worldSize, 0); - for (int i = 1; i < worldSize; ++i) { - sendDispls[i] = sendDispls[i - 1] + sendCounts[i - 1]; - recvDispls[i] = recvDispls[i - 1] + recvCounts[i - 1]; - } - for (size_t i = 0; i < sendRows.size(); ++i) { - int tgt = pDist->iOwnerRank(sendRows[i]); - int pos = sendDispls[tgt] + offsets[tgt]++; - sendBuf[static_cast(pos)].row = sendRows[i]; - sendBuf[static_cast(pos)].val = sendVals[i]; - } - MPI_Datatype tPairType; - { - cPair dummy; - int bl[2] = { 1, 1 }; - MPI_Aint disp[2]; - MPI_Aint base; - MPI_Get_address(&dummy, &base); - MPI_Get_address(&dummy.row, &disp[0]); - MPI_Get_address(&dummy.val, &disp[1]); - disp[0] -= base; - disp[1] -= base; - MPI_Datatype types[2]; - types[0] = MPI_LONG_LONG; - if constexpr (std::is_same_v) - types[1] = MPI_FLOAT; - else if constexpr (std::is_same_v) - types[1] = MPI_DOUBLE; - else if constexpr (std::is_same_v>) { -#ifdef MPI_C_FLOAT_COMPLEX - types[1] = MPI_C_FLOAT_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_FLOAT, &types[1]); - MPI_Type_commit(&types[1]); -#endif - } else if constexpr (std::is_same_v>) { -#ifdef MPI_C_DOUBLE_COMPLEX - types[1] = MPI_C_DOUBLE_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_DOUBLE, &types[1]); - MPI_Type_commit(&types[1]); -#endif - } else { - throw std::runtime_error("SpMV symmetric remote: unsupported scalar type"); - } - MPI_Type_create_struct(2, bl, disp, types, &tPairType); - MPI_Type_commit(&tPairType); - } - MPI_Alltoallv(sendBuf.data(), sendCounts.data(), sendDispls.data(), tPairType, - recvBuf.data(), recvCounts.data(), recvDispls.data(), tPairType, - MPI_COMM_WORLD); - MPI_Type_free(&tPairType); - iSize rowBase2 = A.iGlobalRowBegin(); - for (const cPair& pr : recvBuf) { - iIndex local = pr.row - rowBase2; - if (local >= 0 && local < localRows) - y[static_cast(local)] += pr.val; - } - } -} +// Operations is now a header-only template implementation. +// This file is kept for build system compatibility but contains no implementation. +// See Operations.h for the complete implementation. -void SpMV(const cCSRMatrix& A, std::vector& x) -{ - std::vector y; - SpMV(A, x, y); - x.swap(y); // x now holds local result -} - -void SpMM(const cCSRMatrix& A, const std::vector& X, iSize nCols, std::vector& Y) -{ - iSize globalCols = A.iGlobalColCount(); - iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); - if (X.size() != static_cast(globalCols * nCols)) { - throw std::runtime_error("SpMM: input dense matrix size mismatch"); - } - Y.assign(static_cast(localRows * nCols), static_cast(0)); - const auto& rowPtr = A.vRowPtr(); - const auto& colInd = A.vColInd(); - const auto& values = A.vValues(); - eSymmStorage eSymm = A.eSymmetry(); - iSize rowBase = A.iGlobalRowBegin(); - auto pDist = A.pDistribution(); - int rank = 0, worldSize = 1; - if (pDist) { - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &worldSize); - } - if (eSymm == eNone || eSymm == eSymmFull) { -#ifdef CSR4MPI_USE_OPENMP -#pragma omp parallel for schedule(static) -#endif - 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 c = colInd[static_cast(k)]; - vScalar a = values[static_cast(k)]; - const vScalar* xColBase = &X[static_cast(c)]; - vScalar* yRowBase = &Y[static_cast(r)]; - for (iSize j = 0; j < nCols; ++j) { - yRowBase[static_cast(j * localRows)] += a * xColBase[static_cast(j * globalCols)]; - } - } - } - } else { - bool isLowerStored = (eSymm == eSymmLower); - bool isUpperStored = (eSymm == eSymmUpper); - bool hasRemote = (pDist && worldSize > 1); - // 如果无远程通信：批量本地更新，按行收集再一次性写入 - if (!hasRemote) { -#ifdef CSR4MPI_USE_OPENMP -#pragma omp parallel for schedule(static) -#endif - for (iSize r = 0; r < localRows; ++r) { - iIndex start = rowPtr[static_cast(r)]; - iIndex end = rowPtr[static_cast(r + 1)]; - iIndex gRow = r + rowBase; - struct cUpd { - iIndex idx; - std::vector vals; - }; - std::vector vAdds; - for (iIndex k = start; k < end; ++k) { - iIndex c = colInd[static_cast(k)]; - vScalar a = values[static_cast(k)]; - bool lowerCond = (gRow >= c); - if (c == gRow) { - const vScalar* xColBase = &X[static_cast(c)]; - vScalar* yRowBase = &Y[static_cast(r)]; - for (iSize j = 0; j < nCols; ++j) - yRowBase[static_cast(j * localRows)] += a * xColBase[static_cast(j * globalCols)]; - continue; - } - if ((isLowerStored && lowerCond) || (isUpperStored && !lowerCond)) { - const vScalar* xColBase_c = &X[static_cast(c)]; - vScalar* yRowBase_r = &Y[static_cast(r)]; - for (iSize j = 0; j < nCols; ++j) - yRowBase_r[static_cast(j * localRows)] += a * xColBase_c[static_cast(j * globalCols)]; - iIndex maybeLocal = c - rowBase; - if (maybeLocal >= 0 && maybeLocal < localRows) { - const vScalar* xColBase_r = &X[static_cast(gRow)]; - cUpd upd; - upd.idx = maybeLocal; - upd.vals.resize(static_cast(nCols)); - for (iSize j = 0; j < nCols; ++j) - upd.vals[static_cast(j)] = a * xColBase_r[static_cast(j * globalCols)]; - vAdds.push_back(std::move(upd)); - } - } - } - if (!vAdds.empty()) { -#ifdef CSR4MPI_USE_OPENMP -#pragma omp critical -#endif - for (const auto& u : vAdds) { - vScalar* yRowBase_c = &Y[static_cast(u.idx)]; - for (iSize j = 0; j < nCols; ++j) - yRowBase_c[static_cast(j * localRows)] += u.vals[static_cast(j)]; - } - } - } - return; - } - // 远程路径：行块聚合，针对每个远程镜像行收集长度为 nCols 的向量，减少元数据 - std::unordered_map> mRemoteBlocks; // row -> contributions[nCols] - for (iSize r = 0; r < localRows; ++r) { - iIndex start = rowPtr[static_cast(r)]; - iIndex end = rowPtr[static_cast(r + 1)]; - iIndex gRow = r + rowBase; - for (iIndex k = start; k < end; ++k) { - iIndex c = colInd[static_cast(k)]; - vScalar a = values[static_cast(k)]; - bool lowerCond = (gRow >= c); - if (c == gRow) { - const vScalar* xColBase = &X[static_cast(c)]; - vScalar* yRowBase = &Y[static_cast(r)]; - for (iSize j = 0; j < nCols; ++j) - yRowBase[static_cast(j * localRows)] += a * xColBase[static_cast(j * globalCols)]; - continue; - } - if ((isLowerStored && lowerCond) || (isUpperStored && !lowerCond)) { - // y_r += a * x_c - const vScalar* xColBase_c = &X[static_cast(c)]; - vScalar* yRowBase_r = &Y[static_cast(r)]; - for (iSize j = 0; j < nCols; ++j) - yRowBase_r[static_cast(j * localRows)] += a * xColBase_c[static_cast(j * globalCols)]; - int owner = pDist->iOwnerRank(c); - if (owner == rank) { - iIndex maybeLocal = c - rowBase; - if (maybeLocal >= 0 && maybeLocal < localRows) { - const vScalar* xColBase_r = &X[static_cast(gRow)]; - vScalar* yRowBase_c = &Y[static_cast(maybeLocal)]; - for (iSize j = 0; j < nCols; ++j) - yRowBase_c[static_cast(j * localRows)] += a * xColBase_r[static_cast(j * globalCols)]; - } - } else { - const vScalar* xColBase_r = &X[static_cast(gRow)]; - auto& vecRef = mRemoteBlocks[c]; - if (vecRef.empty()) - vecRef.assign(static_cast(nCols), static_cast(0)); - for (iSize j = 0; j < nCols; ++j) - vecRef[static_cast(j)] += a * xColBase_r[static_cast(j * globalCols)]; - } - } - } - } - // 统计远程行块数量 - std::vector sendRowCounts(worldSize, 0); - for (const auto& kv : mRemoteBlocks) { - int tgt = pDist->iOwnerRank(kv.first); - sendRowCounts[tgt]++; - } - std::vector recvCounts(worldSize, 0); - MPI_Alltoall(sendRowCounts.data(), 1, MPI_INT, recvCounts.data(), 1, MPI_INT, MPI_COMM_WORLD); - int totalSend = 0, totalRecv = 0; - for (int i = 0; i < worldSize; ++i) { - totalSend += sendRowCounts[i]; - totalRecv += recvCounts[i]; - } - // 打包行索引与块值两个向量；行数量与块数量相同 - std::vector sendRowsList; - sendRowsList.reserve(static_cast(totalSend)); - std::vector sendValsList; - sendValsList.reserve(static_cast(totalSend) * static_cast(nCols)); - // rank 分桶（保持顺序，不需稳定，仅需与计数一致） - std::vector bucketInserted(worldSize, 0); - // 临时：按目标 rank 聚合，再按 rank 复制到线性数组 - std::vector> perRankRows(worldSize); - std::vector> perRankVals(worldSize); - for (auto& kv : mRemoteBlocks) { - int tgt = pDist->iOwnerRank(kv.first); - perRankRows[tgt].push_back(kv.first); - auto& vec = kv.second; - perRankVals[tgt].insert(perRankVals[tgt].end(), vec.begin(), vec.end()); - } - for (int rnk = 0; rnk < worldSize; ++rnk) { - sendRowsList.insert(sendRowsList.end(), perRankRows[rnk].begin(), perRankRows[rnk].end()); - sendValsList.insert(sendValsList.end(), perRankVals[rnk].begin(), perRankVals[rnk].end()); - } - // 行索引 Alltoallv - std::vector sendRowDispls(worldSize, 0), recvRowDispls(worldSize, 0); - for (int i = 1; i < worldSize; ++i) { - sendRowDispls[i] = sendRowDispls[i - 1] + sendRowCounts[i - 1]; - recvRowDispls[i] = recvRowDispls[i - 1] + recvCounts[i - 1]; - } - std::vector recvRowsList(static_cast(totalRecv)); - MPI_Alltoallv(sendRowsList.data(), sendRowCounts.data(), sendRowDispls.data(), MPI_LONG_LONG, - recvRowsList.data(), recvCounts.data(), recvRowDispls.data(), MPI_LONG_LONG, MPI_COMM_WORLD); - // 块值 Alltoallv：每个块有 nCols 标量 - std::vector sendValCounts(worldSize, 0), recvValCounts(worldSize, 0), sendValDispls(worldSize, 0), recvValDispls(worldSize, 0); - for (int i = 0; i < worldSize; ++i) { - sendValCounts[i] = sendRowCounts[i] * static_cast(nCols); - recvValCounts[i] = recvCounts[i] * static_cast(nCols); - } - for (int i = 1; i < worldSize; ++i) { - sendValDispls[i] = sendValDispls[i - 1] + sendValCounts[i - 1]; - recvValDispls[i] = recvValDispls[i - 1] + recvValCounts[i - 1]; - } - std::vector recvValsList(static_cast(totalRecv) * static_cast(nCols)); - MPI_Datatype scalarType; - if constexpr (std::is_same_v) - scalarType = MPI_FLOAT; - else if constexpr (std::is_same_v) - scalarType = MPI_DOUBLE; - else if constexpr (std::is_same_v>) { -#ifdef MPI_C_FLOAT_COMPLEX - scalarType = MPI_C_FLOAT_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_FLOAT, &scalarType); - MPI_Type_commit(&scalarType); -#endif - } else if constexpr (std::is_same_v>) { -#ifdef MPI_C_DOUBLE_COMPLEX - scalarType = MPI_C_DOUBLE_COMPLEX; -#else - MPI_Type_contiguous(2, MPI_DOUBLE, &scalarType); - MPI_Type_commit(&scalarType); -#endif - } else { - throw std::runtime_error("SpMM symmetric remote block: unsupported scalar type"); - } - MPI_Alltoallv(sendValsList.data(), sendValCounts.data(), sendValDispls.data(), scalarType, - recvValsList.data(), recvValCounts.data(), recvValDispls.data(), scalarType, - MPI_COMM_WORLD); - // 应用接收块 - for (size_t blk = 0; blk < recvRowsList.size(); ++blk) { - iIndex local = recvRowsList[blk] - rowBase; - if (local >= 0 && local < localRows) { - vScalar* yRowBase = &Y[static_cast(local)]; - const vScalar* vals = &recvValsList[blk * static_cast(nCols)]; - for (iSize j = 0; j < nCols; ++j) { - yRowBase[static_cast(j * localRows)] += vals[static_cast(j)]; - } - } - } - // 释放临时派生类型（若创建了复数 contiguous 类型） -#if !defined(MPI_C_FLOAT_COMPLEX) - if constexpr (std::is_same_v>) { - MPI_Type_free(&scalarType); - } -#endif -#if !defined(MPI_C_DOUBLE_COMPLEX) - if constexpr (std::is_same_v>) { - MPI_Type_free(&scalarType); - } -#endif - } -} +#include "Operations.h" -void SpMMInPlace(const cCSRMatrix& A, std::vector& X, iSize nCols) -{ - std::vector Y; - SpMM(A, X, nCols, Y); - X.swap(Y); // X now holds localRows x nCols result -} -} +// Explicit template instantiations are not needed since the functions are header-only. diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6f5a541..ca0c31f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -39,10 +39,11 @@ add_executable(csr4mpi_tests test_matrix_market.cpp test_symmetric.cpp test_mpi_symmetric_spmm.cpp + test_all_scalar_types.cpp ) target_link_libraries(csr4mpi_tests PRIVATE csr4mpi GTest::gtest_main) -target_compile_definitions(csr4mpi_tests PRIVATE CSR4MPI_SOURCE_DIR="${CMAKE_SOURCE_DIR}" CSR4MPI_VALUE_TYPE=${CSR4MPI_VALUE_TYPE}) +target_compile_definitions(csr4mpi_tests PRIVATE CSR4MPI_SOURCE_DIR="${CMAKE_SOURCE_DIR}") if(CSR4MPI_ENABLE_BLAS) target_compile_definitions(csr4mpi_tests PRIVATE CSR4MPI_USE_BLAS) endif() @@ -229,7 +230,7 @@ if(CSR4MPI_ENABLE_MUMPS) if (DEFINED CSR4MPI_MPI_LINK_TARGET) target_link_libraries(csr4mpi_mumps_tests PRIVATE ${CSR4MPI_MPI_LINK_TARGET}) endif() - target_compile_definitions(csr4mpi_mumps_tests PRIVATE CSR4MPI_SOURCE_DIR="${CMAKE_SOURCE_DIR}" CSR4MPI_VALUE_TYPE=${CSR4MPI_VALUE_TYPE} CSR4MPI_USE_MUMPS) + target_compile_definitions(csr4mpi_mumps_tests PRIVATE CSR4MPI_SOURCE_DIR="${CMAKE_SOURCE_DIR}" CSR4MPI_USE_MUMPS) if (CSR4MPI_MUMPS_DEBUG) target_compile_definitions(csr4mpi_mumps_tests PRIVATE CSR4MPI_MUMPS_DEBUG) endif() diff --git a/tests/test_all_scalar_types.cpp b/tests/test_all_scalar_types.cpp new file mode 100644 index 0000000..467f20a --- /dev/null +++ b/tests/test_all_scalar_types.cpp @@ -0,0 +1,214 @@ +// Test file to verify all 4 scalar types work correctly +// This ensures all scalar types (float, double, complex, complex) compile and run correctly + +#include "CSRMatrix.h" +#include "Operations.h" +#include "MumpsAdapter.h" +#include "Distribution.h" +#include +#include +#include + +using namespace csr4mpi; + +// Template test fixture for testing all scalar types +template +class AllScalarTypesTest : public ::testing::Test { +protected: + using scalar_type = Scalar; + + // Helper to build a simple 3x3 local CSR matrix with rows [0,3) + static cCSRMatrix BuildSimple3x3() + { + iSize iGlobalRowBegin = 0; + iSize iGlobalRowEnd = 3; + iSize iGlobalColCount = 3; + // Row 0: cols 0,2 -> values 1,2 + // Row 1: col 1 -> value 3 + // Row 2: cols 0,1,2 -> values 4,5,6 + std::vector vRowPtr = { 0, 2, 3, 6 }; + std::vector vColInd = { 0, 2, 1, 0, 1, 2 }; + std::vector vValues = { + static_cast(1), static_cast(2), + static_cast(3), + static_cast(4), static_cast(5), static_cast(6) + }; + return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); + } + + // Helper function to check approximate equality for floating-point and complex types + static bool ApproxEqual(Scalar a, Scalar b, double tolerance = 1e-10) { + // std::abs works for both real and complex types + return std::abs(a - b) < tolerance; + } +}; + +// Define type list for all supported scalar types +using ScalarTypes = ::testing::Types, std::complex>; +TYPED_TEST_SUITE(AllScalarTypesTest, ScalarTypes); + +// Test basic matrix construction for all scalar types +TYPED_TEST(AllScalarTypesTest, BasicConstruction) +{ + using Scalar = TypeParam; + + std::vector vRowPtr = { 0, 2, 3 }; + std::vector vColInd = { 0, 2, 1 }; + std::vector vValues = { + static_cast(1.0), + static_cast(2.0), + static_cast(3.0) + }; + + cCSRMatrix cLocal(0, 2, 3, vRowPtr, vColInd, vValues); + + EXPECT_EQ(cLocal.iGlobalRowBegin(), 0); + EXPECT_EQ(cLocal.iGlobalRowEnd(), 2); + EXPECT_EQ(cLocal.vRowPtr().size(), 3u); + EXPECT_EQ(cLocal.vColInd().size(), 3u); + EXPECT_EQ(cLocal.vValues().size(), 3u); +} + +// Test SpMV for all scalar types +TYPED_TEST(AllScalarTypesTest, SpMVBasic) +{ + using Scalar = TypeParam; + + auto A = this->BuildSimple3x3(); + std::vector x = { + static_cast(10), + static_cast(20), + static_cast(30) + }; + std::vector y; + SpMV(A, x, y); + + ASSERT_EQ(y.size(), 3u); + // Row0: 1*10 + 2*30 = 70 + // Row1: 3*20 = 60 + // Row2: 4*10 + 5*20 + 6*30 = 320 + EXPECT_TRUE(this->ApproxEqual(y[0], static_cast(70))); + EXPECT_TRUE(this->ApproxEqual(y[1], static_cast(60))); + EXPECT_TRUE(this->ApproxEqual(y[2], static_cast(320))); +} + +// Test SpMM for all scalar types +TYPED_TEST(AllScalarTypesTest, SpMMBasic) +{ + using Scalar = TypeParam; + + auto A = this->BuildSimple3x3(); + // X: 3x2 column-major -> columns x0={10,20,30}, x1={1,2,3} + std::vector X = { + static_cast(10), static_cast(20), static_cast(30), + static_cast(1), static_cast(2), static_cast(3) + }; + std::vector Y; + SpMM(A, X, 2, Y); + + iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); + ASSERT_EQ(Y.size(), static_cast(localRows * 2)); + + // Column 0 result equals SpMV with x0 + EXPECT_TRUE(this->ApproxEqual(Y[0], static_cast(70))); + EXPECT_TRUE(this->ApproxEqual(Y[1], static_cast(60))); + EXPECT_TRUE(this->ApproxEqual(Y[2], static_cast(320))); + + // Column 1: multiply with x1={1,2,3} + // Row0: 1*1 + 2*3 = 7 + // Row1: 3*2 = 6 + // Row2: 4*1 + 5*2 + 6*3 = 32 + EXPECT_TRUE(this->ApproxEqual(Y[3], static_cast(7))); + EXPECT_TRUE(this->ApproxEqual(Y[4], static_cast(6))); + EXPECT_TRUE(this->ApproxEqual(Y[5], static_cast(32))); +} + +// Test MumpsAdapter export for all scalar types +TYPED_TEST(AllScalarTypesTest, MumpsAdapterExport) +{ + using Scalar = TypeParam; + + std::vector rowPtr = { 0, 2, 3 }; + std::vector colInd = { 0, 2, 1 }; + std::vector values = { + static_cast(10.0), + static_cast(20.0), + static_cast(30.0) + }; + + cCSRMatrix local(0, 2, 3, rowPtr, colInd, values); + cRowDistribution dist = cRowDistribution::CreateBlockDistribution(2, 1, 0); + + std::vector IRN; + std::vector JCN; + std::vector A; + + cMumpsAdapter::ExportLocalBlock(local, dist, IRN, JCN, A); + + ASSERT_EQ(IRN.size(), 3u); + ASSERT_EQ(JCN.size(), 3u); + ASSERT_EQ(A.size(), 3u); + + // Verify 1-based indices + EXPECT_EQ(IRN[0], 1); + EXPECT_EQ(JCN[0], 1); + EXPECT_TRUE(this->ApproxEqual(A[0], static_cast(10.0))); +} + +// Test empty matrix for all scalar types +TYPED_TEST(AllScalarTypesTest, EmptyMatrix) +{ + using Scalar = TypeParam; + + std::vector vRowPtr = { 0 }; + std::vector vColInd; + std::vector vValues; + + cCSRMatrix A(0, 0, 3, vRowPtr, vColInd, vValues); + + EXPECT_EQ(A.iGlobalRowBegin(), 0); + EXPECT_EQ(A.iGlobalRowEnd(), 0); + EXPECT_EQ(A.vValues().size(), 0u); +} + +// Test type traits +TEST(TypeTraitsTest, IsComplex) +{ + EXPECT_FALSE(is_complex_v); + EXPECT_FALSE(is_complex_v); + EXPECT_TRUE(is_complex_v>); + EXPECT_TRUE(is_complex_v>); +} + +TEST(TypeTraitsTest, IsSupportedScalar) +{ + EXPECT_TRUE(is_supported_scalar_v); + EXPECT_TRUE(is_supported_scalar_v); + EXPECT_TRUE(is_supported_scalar_v>); + EXPECT_TRUE(is_supported_scalar_v>); + EXPECT_FALSE(is_supported_scalar_v); + EXPECT_FALSE(is_supported_scalar_v); +} + +TEST(TypeTraitsTest, RealType) +{ + static_assert(std::is_same_v, float>, "real_type_t should be float"); + static_assert(std::is_same_v, double>, "real_type_t should be double"); + static_assert(std::is_same_v>, float>, "real_type_t> should be float"); + static_assert(std::is_same_v>, double>, "real_type_t> should be double"); +} + +// Test type aliases +TEST(TypeAliasTest, MatrixAliases) +{ + // Verify type aliases compile correctly + cCSRMatrixF matF; + cCSRMatrixD matD; + cCSRMatrixCF matCF; + cCSRMatrixCD matCD; + + EXPECT_EQ(matF.iGlobalRowBegin(), 0); + EXPECT_EQ(matD.iGlobalRowBegin(), 0); + EXPECT_EQ(matCF.iGlobalRowBegin(), 0); + EXPECT_EQ(matCD.iGlobalRowBegin(), 0); +} diff --git a/tests/test_blas_placeholder.cpp b/tests/test_blas_placeholder.cpp index 67f7c49..66a8263 100644 --- a/tests/test_blas_placeholder.cpp +++ b/tests/test_blas_placeholder.cpp @@ -4,6 +4,7 @@ #include using namespace csr4mpi; +using Scalar = double; TEST(TestBlasPlaceholder, FlagReflectsOption) { @@ -20,12 +21,12 @@ TEST(TestBlasPlaceholder, SpMMBlasMatchesSpMM) iSize rBegin = 0, rEnd = 3, cCount = 4; std::vector rowPtr = { 0, 2, 3, 5 }; std::vector colInd = { 0, 2, 1, 1, 3 }; - std::vector values = { 1, 2, 3, 4, 5 }; - cCSRMatrix A(rBegin, rEnd, cCount, rowPtr, colInd, values); + std::vector values = { 1, 2, 3, 4, 5 }; + cCSRMatrix A(rBegin, rEnd, cCount, rowPtr, colInd, values); iSize nCols = 2; // Dense RHS (4x2): columns {1,2,3,4} and {5,6,7,8} - std::vector X = { 1, 2, 3, 4, 5, 6, 7, 8 }; - std::vector Y1, Y2; + std::vector X = { 1, 2, 3, 4, 5, 6, 7, 8 }; + std::vector Y1, Y2; SpMM(A, X, nCols, Y1); SpMMBlas(A, X, nCols, Y2); ASSERT_EQ(Y1.size(), Y2.size()); diff --git a/tests/test_csr_matrix.cpp b/tests/test_csr_matrix.cpp index 3a8c257..2a7af3a 100644 --- a/tests/test_csr_matrix.cpp +++ b/tests/test_csr_matrix.cpp @@ -6,6 +6,9 @@ using namespace csr4mpi; +// Test with double as the default scalar type +using Scalar = double; + TEST(CSRMatrixTest, BasicConstruction) { std::vector vRowPtr; @@ -18,12 +21,12 @@ TEST(CSRMatrixTest, BasicConstruction) vColInd.push_back(2); vColInd.push_back(1); - std::vector vValues; - vValues.push_back(static_cast(1.0)); - vValues.push_back(static_cast(2.0)); - vValues.push_back(static_cast(3.0)); + std::vector vValues; + vValues.push_back(static_cast(1.0)); + vValues.push_back(static_cast(2.0)); + vValues.push_back(static_cast(3.0)); - cCSRMatrix cLocal(0, 2, 3, vRowPtr, vColInd, vValues); + cCSRMatrix cLocal(0, 2, 3, vRowPtr, vColInd, vValues); EXPECT_EQ(cLocal.iGlobalRowBegin(), 0); EXPECT_EQ(cLocal.iGlobalRowEnd(), 2); @@ -43,40 +46,40 @@ TEST(CSRCommTest, LocalAssembleNoMPI) vColInd.push_back(0); vColInd.push_back(1); - std::vector vValues; - vValues.push_back(static_cast(1.0)); - vValues.push_back(static_cast(2.0)); + std::vector vValues; + vValues.push_back(static_cast(1.0)); + vValues.push_back(static_cast(2.0)); - cCSRMatrix cLocal(0, 2, 2, vRowPtr, vColInd, vValues); + cCSRMatrix cLocal(0, 2, 2, vRowPtr, vColInd, vValues); cRowDistribution cDist = cRowDistribution::CreateBlockDistribution(2, 1, 0); - std::vector vEntries; - cRemoteEntry cEntry0; + std::vector> vEntries; + cRemoteEntry cEntry0; cEntry0.m_iGlobalRow = 0; cEntry0.m_iGlobalCol = 0; - cEntry0.m_vValue = static_cast(3.0); + cEntry0.m_vValue = static_cast(3.0); vEntries.push_back(cEntry0); - cRemoteEntry cEntry1; + cRemoteEntry cEntry1; cEntry1.m_iGlobalRow = 1; cEntry1.m_iGlobalCol = 1; - cEntry1.m_vValue = static_cast(4.0); + cEntry1.m_vValue = static_cast(4.0); vEntries.push_back(cEntry1); - cCommPattern cPattern; + cCommPattern cPattern; cPattern.Build(vEntries, cDist, 0, 1); - std::vector vContrib; - vContrib.push_back(static_cast(3.0)); - vContrib.push_back(static_cast(4.0)); + std::vector vContrib; + vContrib.push_back(static_cast(3.0)); + vContrib.push_back(static_cast(4.0)); MPI_Init(nullptr, nullptr); - cCSRComm::Assemble(cLocal, cPattern, vContrib, MPI_COMM_WORLD); + cCSRComm::Assemble(cLocal, cPattern, vContrib, MPI_COMM_WORLD); MPI_Finalize(); - const std::vector& vNewValues = cLocal.vValues(); + const std::vector& vNewValues = cLocal.vValues(); - EXPECT_EQ(vNewValues[0], static_cast(4.0)); - EXPECT_EQ(vNewValues[1], static_cast(6.0)); + EXPECT_EQ(vNewValues[0], static_cast(4.0)); + EXPECT_EQ(vNewValues[1], static_cast(6.0)); } diff --git a/tests/test_distribution_pattern.cpp b/tests/test_distribution_pattern.cpp index 3e67b84..480200a 100644 --- a/tests/test_distribution_pattern.cpp +++ b/tests/test_distribution_pattern.cpp @@ -6,6 +6,9 @@ using namespace csr4mpi; +// Test with double as the default scalar type +using Scalar = double; + TEST(RowDistributionTest, BlockEvenSplit) { cRowDistribution dist = cRowDistribution::CreateBlockDistribution(8, 4, 2); @@ -43,63 +46,63 @@ TEST(CommPatternTest, DuplicateEntriesAccumulate) // Local 2x2 matrix with both rows owned. std::vector rowPtr = { 0, 2, 4 }; std::vector colInd = { 0, 1, 0, 1 }; - std::vector values = { 1.0, 2.0, 3.0, 4.0 }; - cCSRMatrix local(0, 2, 2, rowPtr, colInd, values); + std::vector values = { 1.0, 2.0, 3.0, 4.0 }; + cCSRMatrix local(0, 2, 2, rowPtr, colInd, values); cRowDistribution dist = cRowDistribution::CreateBlockDistribution(2, 1, 0); - std::vector entries; - cRemoteEntry e1; + std::vector> entries; + cRemoteEntry e1; e1.m_iGlobalRow = 0; e1.m_iGlobalCol = 0; - e1.m_vValue = static_cast(5.0); - cRemoteEntry e2; + e1.m_vValue = static_cast(5.0); + cRemoteEntry e2; e2.m_iGlobalRow = 0; e2.m_iGlobalCol = 0; - e2.m_vValue = static_cast(6.0); - cRemoteEntry e3; + e2.m_vValue = static_cast(6.0); + cRemoteEntry e3; e3.m_iGlobalRow = 1; e3.m_iGlobalCol = 1; - e3.m_vValue = static_cast(7.0); + e3.m_vValue = static_cast(7.0); entries.push_back(e1); entries.push_back(e2); entries.push_back(e3); - cCommPattern pattern; + cCommPattern pattern; pattern.Build(entries, dist, 0, 1); - std::vector contrib; - contrib.push_back(static_cast(5.0)); - contrib.push_back(static_cast(6.0)); - contrib.push_back(static_cast(7.0)); + std::vector contrib; + contrib.push_back(static_cast(5.0)); + contrib.push_back(static_cast(6.0)); + contrib.push_back(static_cast(7.0)); MPI_Init(nullptr, nullptr); - cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); MPI_Finalize(); - const std::vector& newVals = local.vValues(); + const std::vector& newVals = local.vValues(); // Original (row0,col0)=1 plus 5 + 6 => 12; (row1,col1)=4 plus 7 => 11 - EXPECT_EQ(newVals[0], static_cast(12.0)); - EXPECT_EQ(newVals[3], static_cast(11.0)); + EXPECT_EQ(newVals[0], static_cast(12.0)); + EXPECT_EQ(newVals[3], static_cast(11.0)); } TEST(CommPatternTest, EmptyPatternNoChange) { std::vector rowPtr = { 0, 1 }; std::vector colInd = { 0 }; - std::vector values = { 2.5 }; - cCSRMatrix local(0, 1, 1, rowPtr, colInd, values); + std::vector values = { 2.5 }; + cCSRMatrix local(0, 1, 1, rowPtr, colInd, values); cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); - std::vector entries; // empty + std::vector> entries; // empty - cCommPattern pattern; + cCommPattern pattern; pattern.Build(entries, dist, 0, 1); - std::vector contrib; // empty + std::vector contrib; // empty MPI_Init(nullptr, nullptr); - cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); + cCSRComm::Assemble(local, pattern, contrib, MPI_COMM_WORLD); MPI_Finalize(); - EXPECT_EQ(local.vValues()[0], static_cast(2.5)); + EXPECT_EQ(local.vValues()[0], static_cast(2.5)); } diff --git a/tests/test_matrix_market.cpp b/tests/test_matrix_market.cpp index 58defd4..2009a7e 100644 --- a/tests/test_matrix_market.cpp +++ b/tests/test_matrix_market.cpp @@ -5,31 +5,32 @@ #include using namespace csr4mpi; +using Scalar = double; -static cCSRMatrix BuildFromMM(const std::string& path) +static cCSRMatrix BuildFromMM(const std::string& path) { std::vector rowPtr, colInd; - std::vector values; + std::vector values; iIndex rows = 0, cols = 0; bool ok = LoadMatrixMarket(path, rowPtr, colInd, values, rows, cols); if (!ok) throw std::runtime_error("Failed to load MatrixMarket file: " + path); - return cCSRMatrix(0, rows, cols, rowPtr, colInd, values); + return cCSRMatrix(0, rows, cols, rowPtr, colInd, values); } TEST(MatrixMarketTest, Lap5SpMVOnes) { std::string path = std::string(CSR4MPI_SOURCE_DIR) + "/tests/data/lap5.mtx"; auto A = BuildFromMM(path); - std::vector x(static_cast(A.iGlobalColCount()), static_cast(1)); - std::vector y; + std::vector x(static_cast(A.iGlobalColCount()), static_cast(1)); + std::vector y; SpMV(A, x, y); ASSERT_EQ(y.size(), static_cast(5)); - EXPECT_EQ(y[0], static_cast(1)); - EXPECT_EQ(y[1], static_cast(0)); - EXPECT_EQ(y[2], static_cast(0)); - EXPECT_EQ(y[3], static_cast(0)); - EXPECT_EQ(y[4], static_cast(1)); + EXPECT_EQ(y[0], static_cast(1)); + EXPECT_EQ(y[1], static_cast(0)); + EXPECT_EQ(y[2], static_cast(0)); + EXPECT_EQ(y[3], static_cast(0)); + EXPECT_EQ(y[4], static_cast(1)); } TEST(MatrixMarketTest, Small4DuplicateAccumulation) @@ -45,7 +46,7 @@ TEST(MatrixMarketTest, Small4DuplicateAccumulation) iIndex start = rowPtr[static_cast(r)]; iIndex end = rowPtr[static_cast(r + 1)]; bool found = false; - vScalar val14 = static_cast(0); + Scalar val14 = static_cast(0); for (iIndex k = start; k < end; ++k) { if (colInd[static_cast(k)] == 3) { found = true; @@ -54,16 +55,16 @@ TEST(MatrixMarketTest, Small4DuplicateAccumulation) } } ASSERT_TRUE(found); - EXPECT_EQ(val14, static_cast(14)); + EXPECT_EQ(val14, static_cast(14)); - std::vector x { static_cast(1), static_cast(2), static_cast(3), static_cast(4) }; - std::vector y; + std::vector x { static_cast(1), static_cast(2), static_cast(3), static_cast(4) }; + std::vector y; SpMV(A, x, y); ASSERT_EQ(y.size(), static_cast(4)); - EXPECT_EQ(y[0], static_cast(19)); // 10*1 + 3*3 - EXPECT_EQ(y[1], static_cast(10)); // 5*2 - EXPECT_EQ(y[2], static_cast(2 + 14 * 4)); // 2*1 +14*4=58 - EXPECT_EQ(y[3], static_cast(4)); // 1*4 + EXPECT_EQ(y[0], static_cast(19)); // 10*1 + 3*3 + EXPECT_EQ(y[1], static_cast(10)); // 5*2 + EXPECT_EQ(y[2], static_cast(2 + 14 * 4)); // 2*1 +14*4=58 + EXPECT_EQ(y[3], static_cast(4)); // 1*4 } TEST(MatrixMarketTest, Dup3AccumulationAndSpMV) @@ -77,7 +78,7 @@ TEST(MatrixMarketTest, Dup3AccumulationAndSpMV) // locate (0,0) iIndex start = rowPtr[0]; iIndex end = rowPtr[1]; - vScalar v00 = static_cast(0); + Scalar v00 = static_cast(0); bool f = false; for (iIndex k = start; k < end; ++k) { if (colInd[static_cast(k)] == 0) { @@ -87,12 +88,12 @@ TEST(MatrixMarketTest, Dup3AccumulationAndSpMV) } } ASSERT_TRUE(f); - EXPECT_EQ(v00, static_cast(3)); - std::vector x { static_cast(1), static_cast(1), static_cast(1) }; - std::vector y; + EXPECT_EQ(v00, static_cast(3)); + std::vector x { static_cast(1), static_cast(1), static_cast(1) }; + std::vector y; SpMV(A, x, y); ASSERT_EQ(y.size(), static_cast(3)); - EXPECT_EQ(y[0], static_cast(7)); - EXPECT_EQ(y[1], static_cast(5)); - EXPECT_EQ(y[2], static_cast(13)); + EXPECT_EQ(y[0], static_cast(7)); + EXPECT_EQ(y[1], static_cast(5)); + EXPECT_EQ(y[2], static_cast(13)); } diff --git a/tests/test_mpi_symmetric_spmm.cpp b/tests/test_mpi_symmetric_spmm.cpp index 8722c07..319a1a3 100644 --- a/tests/test_mpi_symmetric_spmm.cpp +++ b/tests/test_mpi_symmetric_spmm.cpp @@ -15,6 +15,7 @@ #include using namespace csr4mpi; +using Scalar = double; // 分布式对称（下三角存储）矩阵 4x4： // Lower stored entries per row (owned rows only): @@ -57,49 +58,49 @@ TEST(MpiSymmetricSpMMTest, DistributedLowerExpansionGatherMM) // Build CSR for owned rows only (lower triangle entries on those rows) std::vector vRowPtr; std::vector vColInd; - std::vector vValues; + std::vector vValues; vRowPtr.push_back(0); for (iIndex gRow = iGlobalRowBegin; gRow < iGlobalRowEnd; ++gRow) { if (gRow == 0) { vColInd.push_back(0); - vValues.push_back(static_cast(2)); + vValues.push_back(static_cast(2)); } else if (gRow == 1) { vColInd.push_back(0); - vValues.push_back(static_cast(1)); + vValues.push_back(static_cast(1)); vColInd.push_back(1); - vValues.push_back(static_cast(3)); + vValues.push_back(static_cast(3)); } else if (gRow == 2) { vColInd.push_back(1); - vValues.push_back(static_cast(4)); + vValues.push_back(static_cast(4)); vColInd.push_back(2); - vValues.push_back(static_cast(5)); + vValues.push_back(static_cast(5)); } else if (gRow == 3) { vColInd.push_back(2); - vValues.push_back(static_cast(6)); + vValues.push_back(static_cast(6)); vColInd.push_back(3); - vValues.push_back(static_cast(7)); + vValues.push_back(static_cast(7)); } vRowPtr.push_back(static_cast(vColInd.size())); } - cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, 4, vRowPtr, vColInd, vValues, eSymmLower); + cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, 4, vRowPtr, vColInd, vValues, eSymmLower); A.AttachDistribution(dist); // Dense X 4x2 col-major - std::vector X = { static_cast(1), static_cast(1), static_cast(1), static_cast(1), - static_cast(2), static_cast(3), static_cast(4), static_cast(5) }; - std::vector Y; + std::vector X = { static_cast(1), static_cast(1), static_cast(1), static_cast(1), + static_cast(2), static_cast(3), static_cast(4), static_cast(5) }; + std::vector Y; SpMM(A, X, 2, Y); ASSERT_EQ(Y.size(), static_cast(localRows * 2)); // Expected slices if (rank == 0) { - EXPECT_EQ(Y[0], static_cast(3)); - EXPECT_EQ(Y[1], static_cast(8)); - EXPECT_EQ(Y[2], static_cast(7)); - EXPECT_EQ(Y[3], static_cast(27)); + EXPECT_EQ(Y[0], static_cast(3)); + EXPECT_EQ(Y[1], static_cast(8)); + EXPECT_EQ(Y[2], static_cast(7)); + EXPECT_EQ(Y[3], static_cast(27)); } else { - EXPECT_EQ(Y[0], static_cast(15)); - EXPECT_EQ(Y[1], static_cast(13)); - EXPECT_EQ(Y[2], static_cast(62)); - EXPECT_EQ(Y[3], static_cast(59)); + EXPECT_EQ(Y[0], static_cast(15)); + EXPECT_EQ(Y[1], static_cast(13)); + EXPECT_EQ(Y[2], static_cast(62)); + EXPECT_EQ(Y[3], static_cast(59)); } int finFlag = 0; MPI_Finalized(&finFlag); @@ -130,55 +131,55 @@ TEST(MpiSymmetricSpMMTest, DistributedLowerExpansionGatherMM4Proc) iSize localRows = iGlobalRowEnd - iGlobalRowBegin; std::vector vRowPtr; std::vector vColInd; - std::vector vValues; + std::vector vValues; vRowPtr.push_back(0); for (iIndex gRow = iGlobalRowBegin; gRow < iGlobalRowEnd; ++gRow) { if (gRow == 0) { vColInd.push_back(0); - vValues.push_back(static_cast(2)); + vValues.push_back(static_cast(2)); } else if (gRow == 1) { vColInd.push_back(0); - vValues.push_back(static_cast(1)); + vValues.push_back(static_cast(1)); vColInd.push_back(1); - vValues.push_back(static_cast(3)); + vValues.push_back(static_cast(3)); } else if (gRow == 2) { vColInd.push_back(1); - vValues.push_back(static_cast(4)); + vValues.push_back(static_cast(4)); vColInd.push_back(2); - vValues.push_back(static_cast(5)); + vValues.push_back(static_cast(5)); } else if (gRow == 3) { vColInd.push_back(2); - vValues.push_back(static_cast(6)); + vValues.push_back(static_cast(6)); vColInd.push_back(3); - vValues.push_back(static_cast(7)); + vValues.push_back(static_cast(7)); } vRowPtr.push_back(static_cast(vColInd.size())); } - cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, 4, vRowPtr, vColInd, vValues, eSymmLower); + cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, 4, vRowPtr, vColInd, vValues, eSymmLower); A.AttachDistribution(dist); // Dense X 4x2 col-major - std::vector X = { static_cast(1), static_cast(1), static_cast(1), static_cast(1), - static_cast(2), static_cast(3), static_cast(4), static_cast(5) }; - std::vector Y; + std::vector X = { static_cast(1), static_cast(1), static_cast(1), static_cast(1), + static_cast(2), static_cast(3), static_cast(4), static_cast(5) }; + std::vector Y; SpMM(A, X, 2, Y); ASSERT_EQ(Y.size(), static_cast(localRows * 2)); // Expected per-rank slice switch (rank) { case 0: - EXPECT_EQ(Y[0], static_cast(3)); - EXPECT_EQ(Y[1], static_cast(7)); + EXPECT_EQ(Y[0], static_cast(3)); + EXPECT_EQ(Y[1], static_cast(7)); break; case 1: - EXPECT_EQ(Y[0], static_cast(8)); - EXPECT_EQ(Y[1], static_cast(27)); + EXPECT_EQ(Y[0], static_cast(8)); + EXPECT_EQ(Y[1], static_cast(27)); break; case 2: - EXPECT_EQ(Y[0], static_cast(15)); - EXPECT_EQ(Y[1], static_cast(62)); + EXPECT_EQ(Y[0], static_cast(15)); + EXPECT_EQ(Y[1], static_cast(62)); break; case 3: - EXPECT_EQ(Y[0], static_cast(13)); - EXPECT_EQ(Y[1], static_cast(59)); + EXPECT_EQ(Y[0], static_cast(13)); + EXPECT_EQ(Y[1], static_cast(59)); break; } int finFlag = 0; diff --git a/tests/test_mumps_adapter.cpp b/tests/test_mumps_adapter.cpp index 3585ad3..4e66d5b 100644 --- a/tests/test_mumps_adapter.cpp +++ b/tests/test_mumps_adapter.cpp @@ -4,25 +4,26 @@ #include using namespace csr4mpi; +using Scalar = double; TEST(MumpsAdapterTest, ExportLocalBlockSimple) { // 2x3 local block (rows 0..2 global with 0-based in matrix definition). std::vector rowPtr = { 0, 2, 3 }; std::vector colInd = { 0, 2, 1 }; - std::vector values; - values.push_back(static_cast(10.0)); - values.push_back(static_cast(20.0)); - values.push_back(static_cast(30.0)); + std::vector values; + values.push_back(static_cast(10.0)); + values.push_back(static_cast(20.0)); + values.push_back(static_cast(30.0)); - cCSRMatrix local(0, 2, 3, rowPtr, colInd, values); + cCSRMatrix local(0, 2, 3, rowPtr, colInd, values); cRowDistribution dist = cRowDistribution::CreateBlockDistribution(2, 1, 0); std::vector IRN; std::vector JCN; - std::vector A; + std::vector A; - cMumpsAdapter::ExportLocalBlock(local, dist, IRN, JCN, A); + cMumpsAdapter::ExportLocalBlock(local, dist, IRN, JCN, A); ASSERT_EQ(IRN.size(), A.size()); ASSERT_EQ(JCN.size(), A.size()); @@ -31,13 +32,13 @@ TEST(MumpsAdapterTest, ExportLocalBlockSimple) // Row 1: (1,1)=10, (1,3)=20; Row 2: (2,2)=30 (ordering follows CSR traversal) EXPECT_EQ(IRN[0], 1); EXPECT_EQ(JCN[0], 1); - EXPECT_EQ(A[0], static_cast(10.0)); + EXPECT_EQ(A[0], static_cast(10.0)); EXPECT_EQ(IRN[1], 1); EXPECT_EQ(JCN[1], 3); - EXPECT_EQ(A[1], static_cast(20.0)); + EXPECT_EQ(A[1], static_cast(20.0)); EXPECT_EQ(IRN[2], 2); EXPECT_EQ(JCN[2], 2); - EXPECT_EQ(A[2], static_cast(30.0)); + EXPECT_EQ(A[2], static_cast(30.0)); } TEST(MumpsAdapterTest, ExportLocalBlockNonZeroOffset) @@ -47,17 +48,17 @@ TEST(MumpsAdapterTest, ExportLocalBlockNonZeroOffset) // Row 2: cols 1,3 (values 5,7) Row 3: col 2 (value 9) std::vector vRowPtr { 0, 2, 3 }; std::vector vColInd { 1, 3, 2 }; // 0-based column indices - std::vector vValues; - vValues.push_back(static_cast(5.0)); - vValues.push_back(static_cast(7.0)); - vValues.push_back(static_cast(9.0)); - cCSRMatrix cLocal(2, 4, 6, vRowPtr, vColInd, vValues); + std::vector vValues; + vValues.push_back(static_cast(5.0)); + vValues.push_back(static_cast(7.0)); + vValues.push_back(static_cast(9.0)); + cCSRMatrix cLocal(2, 4, 6, vRowPtr, vColInd, vValues); cRowDistribution dist = cRowDistribution::CreateBlockDistribution(6, 3, 1); // rank=1 owns [2,4) std::vector vIRN; std::vector vJCN; - std::vector vA; - cMumpsAdapter::ExportLocalBlock(cLocal, dist, vIRN, vJCN, vA); + std::vector vA; + cMumpsAdapter::ExportLocalBlock(cLocal, dist, vIRN, vJCN, vA); ASSERT_EQ(vIRN.size(), 3u); ASSERT_EQ(vJCN.size(), 3u); ASSERT_EQ(vA.size(), 3u); @@ -69,24 +70,24 @@ TEST(MumpsAdapterTest, ExportLocalBlockNonZeroOffset) EXPECT_EQ(vJCN[0], 2); EXPECT_EQ(vJCN[1], 4); EXPECT_EQ(vJCN[2], 3); - EXPECT_EQ(vA[0], static_cast(5.0)); - EXPECT_EQ(vA[1], static_cast(7.0)); - EXPECT_EQ(vA[2], static_cast(9.0)); + EXPECT_EQ(vA[0], static_cast(5.0)); + EXPECT_EQ(vA[1], static_cast(7.0)); + EXPECT_EQ(vA[2], static_cast(9.0)); } TEST(MumpsAdapterTest, ExportLocalBlockIntoPreallocated) { std::vector vRowPtr { 0, 1, 3 }; // two rows: first has 1 nnz, second has 2 nnz std::vector vColInd { 0, 1, 2 }; - std::vector vValues { static_cast(2.0), static_cast(4.0), static_cast(8.0) }; - cCSRMatrix cLocal(4, 6, 10, vRowPtr, vColInd, vValues); // global rows [4,6) + std::vector vValues { static_cast(2.0), static_cast(4.0), static_cast(8.0) }; + cCSRMatrix cLocal(4, 6, 10, vRowPtr, vColInd, vValues); // global rows [4,6) cRowDistribution dist = cRowDistribution::CreateBlockDistribution(10, 5, 2); // Need rank=2 block starting at row 4? For 10 rows & 5 ranks each gets 2 rows => rank2 owns [4,6) // Preallocate int iLocalNNZ = static_cast(vRowPtr.back()); std::vector vIRN(iLocalNNZ); std::vector vJCN(iLocalNNZ); - std::vector vA(iLocalNNZ); - iSize iWritten = cMumpsAdapter::ExportLocalBlockInto(cLocal, dist, vIRN.data(), vJCN.data(), vA.data()); + std::vector vA(iLocalNNZ); + iSize iWritten = cMumpsAdapter::ExportLocalBlockInto(cLocal, dist, vIRN.data(), vJCN.data(), vA.data()); ASSERT_EQ(iWritten, vRowPtr.back()); // Expected rows: global 4,5 -> 1-based 5,5,6 EXPECT_EQ(vIRN[0], 5); @@ -99,47 +100,47 @@ TEST(MumpsAdapterTest, ExportLocalBlockIntoPreallocated) EXPECT_EQ(vJCN[1], 2); EXPECT_EQ(vIRN[2], 6); EXPECT_EQ(vJCN[2], 3); - EXPECT_EQ(vA[0], static_cast(2.0)); - EXPECT_EQ(vA[1], static_cast(4.0)); - EXPECT_EQ(vA[2], static_cast(8.0)); + EXPECT_EQ(vA[0], static_cast(2.0)); + EXPECT_EQ(vA[1], static_cast(4.0)); + EXPECT_EQ(vA[2], static_cast(8.0)); } TEST(MumpsAdapterTest, ExportLocalBlockAppendPreservesExisting) { std::vector vRowPtr { 0, 2 }; std::vector vColInd { 0, 2 }; - std::vector vValues { static_cast(1.0), static_cast(3.0) }; - cCSRMatrix cLocal(0, 1, 5, vRowPtr, vColInd, vValues); + std::vector vValues { static_cast(1.0), static_cast(3.0) }; + cCSRMatrix cLocal(0, 1, 5, vRowPtr, vColInd, vValues); cRowDistribution dist = cRowDistribution::CreateBlockDistribution(1, 1, 0); std::vector vIRN { -999 }; std::vector vJCN { -888 }; - std::vector vA { static_cast(-1.0) }; - cMumpsAdapter::ExportLocalBlock(cLocal, dist, vIRN, vJCN, vA); + std::vector vA { static_cast(-1.0) }; + cMumpsAdapter::ExportLocalBlock(cLocal, dist, vIRN, vJCN, vA); ASSERT_EQ(vIRN.size(), 1u + 2u); // Existing sentinel preserved EXPECT_EQ(vIRN[0], -999); EXPECT_EQ(vJCN[0], -888); - EXPECT_EQ(vA[0], static_cast(-1.0)); + EXPECT_EQ(vA[0], static_cast(-1.0)); // Appended entries (row 0 -> 1-based 1) EXPECT_EQ(vIRN[1], 1); EXPECT_EQ(vJCN[1], 1); - EXPECT_EQ(vA[1], static_cast(1.0)); + EXPECT_EQ(vA[1], static_cast(1.0)); EXPECT_EQ(vIRN[2], 1); EXPECT_EQ(vJCN[2], 3); - EXPECT_EQ(vA[2], static_cast(3.0)); + EXPECT_EQ(vA[2], static_cast(3.0)); } TEST(MumpsAdapterTest, ExportLocalBlockEmpty) { std::vector vRowPtr { 0 }; // zero rows std::vector vColInd; - std::vector vValues; - cCSRMatrix cLocal(0, 0, 4, vRowPtr, vColInd, vValues); + std::vector vValues; + cCSRMatrix cLocal(0, 0, 4, vRowPtr, vColInd, vValues); cRowDistribution dist = cRowDistribution::CreateBlockDistribution(0, 1, 0); std::vector vIRN; std::vector vJCN; - std::vector vA; - cMumpsAdapter::ExportLocalBlock(cLocal, dist, vIRN, vJCN, vA); + std::vector vA; + cMumpsAdapter::ExportLocalBlock(cLocal, dist, vIRN, vJCN, vA); EXPECT_TRUE(vIRN.empty()); EXPECT_TRUE(vJCN.empty()); EXPECT_TRUE(vA.empty()); diff --git a/tests/test_operations.cpp b/tests/test_operations.cpp index 3230cfa..6ccbb62 100644 --- a/tests/test_operations.cpp +++ b/tests/test_operations.cpp @@ -4,8 +4,11 @@ using namespace csr4mpi; +// Test with double as the default scalar type +using Scalar = double; + // Helper to build a simple 3x3 local CSR matrix with rows [0,3) -static cCSRMatrix BuildSimple3x3() +static cCSRMatrix BuildSimple3x3() { iSize iGlobalRowBegin = 0; iSize iGlobalRowEnd = 3; @@ -15,34 +18,34 @@ static cCSRMatrix BuildSimple3x3() // Row 2: cols 0,1,2 -> values 4,5,6 std::vector vRowPtr = { 0, 2, 3, 6 }; std::vector vColInd = { 0, 2, 1, 0, 1, 2 }; - std::vector vValues = { 1, 2, 3, 4, 5, 6 }; - return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); + std::vector vValues = { 1, 2, 3, 4, 5, 6 }; + return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); } TEST(TestOperations, SpMVBasic) { auto A = BuildSimple3x3(); - std::vector x = { 10, 20, 30 }; - std::vector y; + std::vector x = { 10, 20, 30 }; + std::vector y; SpMV(A, x, y); ASSERT_EQ(y.size(), 3u); // Row0: 1*10 + 2*30 = 70 // Row1: 3*20 = 60 // Row2: 4*10 +5*20 +6*30 = 4*10 +5*20 +180 = 40+100+180 = 320 - EXPECT_EQ(y[0], static_cast(70)); - EXPECT_EQ(y[1], static_cast(60)); - EXPECT_EQ(y[2], static_cast(320)); + EXPECT_EQ(y[0], static_cast(70)); + EXPECT_EQ(y[1], static_cast(60)); + EXPECT_EQ(y[2], static_cast(320)); } TEST(TestOperations, SpMVInPlace) { auto A = BuildSimple3x3(); - std::vector x = { 10, 20, 30 }; + std::vector x = { 10, 20, 30 }; SpMV(A, x); ASSERT_EQ(x.size(), 3u); - EXPECT_EQ(x[0], static_cast(70)); - EXPECT_EQ(x[1], static_cast(60)); - EXPECT_EQ(x[2], static_cast(320)); + EXPECT_EQ(x[0], static_cast(70)); + EXPECT_EQ(x[1], static_cast(60)); + EXPECT_EQ(x[2], static_cast(320)); } TEST(TestOperations, SpMMBasic) @@ -50,35 +53,35 @@ TEST(TestOperations, SpMMBasic) auto A = BuildSimple3x3(); // X: 3x2 column-major -> columns x0={10,20,30}, x1={1,2,3} // Stored as [10,20,30, 1,2,3] - std::vector X = { 10, 20, 30, 1, 2, 3 }; - std::vector Y; + std::vector X = { 10, 20, 30, 1, 2, 3 }; + std::vector Y; SpMM(A, X, 2, Y); iSize localRows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); ASSERT_EQ(Y.size(), static_cast(localRows * 2)); // Column 0 result equals SpMV with x0 - EXPECT_EQ(Y[0], static_cast(70)); - EXPECT_EQ(Y[1], static_cast(60)); - EXPECT_EQ(Y[2], static_cast(320)); + EXPECT_EQ(Y[0], static_cast(70)); + EXPECT_EQ(Y[1], static_cast(60)); + EXPECT_EQ(Y[2], static_cast(320)); // Column 1: multiply with x1={1,2,3} // Row0: 1*1 + 2*3 = 7 // Row1: 3*2 = 6 // Row2: 4*1 +5*2 +6*3 = 4 +10 +18 = 32 - EXPECT_EQ(Y[3], static_cast(7)); - EXPECT_EQ(Y[4], static_cast(6)); - EXPECT_EQ(Y[5], static_cast(32)); + EXPECT_EQ(Y[3], static_cast(7)); + EXPECT_EQ(Y[4], static_cast(6)); + EXPECT_EQ(Y[5], static_cast(32)); } TEST(TestOperations, SpMMInPlace) { auto A = BuildSimple3x3(); - std::vector X = { 10, 20, 30, 1, 2, 3 }; + std::vector X = { 10, 20, 30, 1, 2, 3 }; SpMMInPlace(A, X, 2); - EXPECT_EQ(X[0], static_cast(70)); - EXPECT_EQ(X[1], static_cast(60)); - EXPECT_EQ(X[2], static_cast(320)); - EXPECT_EQ(X[3], static_cast(7)); - EXPECT_EQ(X[4], static_cast(6)); - EXPECT_EQ(X[5], static_cast(32)); + EXPECT_EQ(X[0], static_cast(70)); + EXPECT_EQ(X[1], static_cast(60)); + EXPECT_EQ(X[2], static_cast(320)); + EXPECT_EQ(X[3], static_cast(7)); + EXPECT_EQ(X[4], static_cast(6)); + EXPECT_EQ(X[5], static_cast(32)); } TEST(TestOperations, SpMVEmptyRow) @@ -87,14 +90,14 @@ TEST(TestOperations, SpMVEmptyRow) iSize iGlobalRowBegin = 0, iGlobalRowEnd = 2, iGlobalColCount = 3; std::vector vRowPtr = { 0, 2, 2 }; std::vector vColInd = { 0, 2 }; - std::vector vValues = { 5, 6 }; - cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); - std::vector x = { 1, 2, 3 }; - std::vector y; + std::vector vValues = { 5, 6 }; + cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); + std::vector x = { 1, 2, 3 }; + std::vector y; SpMV(A, x, y); ASSERT_EQ(y.size(), 2u); - EXPECT_EQ(y[0], static_cast(5 * 1 + 6 * 3)); - EXPECT_EQ(y[1], static_cast(0)); + EXPECT_EQ(y[0], static_cast(5 * 1 + 6 * 3)); + EXPECT_EQ(y[1], static_cast(0)); } TEST(TestOperations, SpMMEmptyRow) @@ -102,17 +105,17 @@ TEST(TestOperations, SpMMEmptyRow) iSize iGlobalRowBegin = 0, iGlobalRowEnd = 2, iGlobalColCount = 3; std::vector vRowPtr = { 0, 2, 2 }; std::vector vColInd = { 0, 2 }; - std::vector vValues = { 5, 6 }; - cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); + std::vector vValues = { 5, 6 }; + cCSRMatrix A(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); // 3x2 dense - std::vector X = { 1, 2, 3, 4, 5, 6 }; - std::vector Y; + std::vector X = { 1, 2, 3, 4, 5, 6 }; + std::vector Y; SpMM(A, X, 2, Y); ASSERT_EQ(Y.size(), 4u); // Column0: row0 5*1 +6*3=23; row1 0 - EXPECT_EQ(Y[0], static_cast(23)); - EXPECT_EQ(Y[1], static_cast(0)); + EXPECT_EQ(Y[0], static_cast(23)); + EXPECT_EQ(Y[1], static_cast(0)); // Column1: row0 5*4 +6*6 = 20+36=56; row1 0 - EXPECT_EQ(Y[2], static_cast(56)); - EXPECT_EQ(Y[3], static_cast(0)); + EXPECT_EQ(Y[2], static_cast(56)); + EXPECT_EQ(Y[3], static_cast(0)); } diff --git a/tests/test_operations_correctness.cpp b/tests/test_operations_correctness.cpp index fc82b09..051fa1a 100644 --- a/tests/test_operations_correctness.cpp +++ b/tests/test_operations_correctness.cpp @@ -3,51 +3,52 @@ #include using namespace csr4mpi; +using Scalar = double; // Build a 4x5 sparse matrix with pattern: // Row0: 1 0 2 0 0 // Row1: 0 0 0 3 0 // Row2: 4 5 0 0 6 // Row3: 0 0 7 0 8 -static cCSRMatrix BuildMatrix4x5() +static cCSRMatrix BuildMatrix4x5() { iSize iGlobalRowBegin = 0; iSize iGlobalRowEnd = 4; iSize iGlobalColCount = 5; std::vector vRowPtr = { 0, 2, 3, 6, 8 }; std::vector vColInd = { 0, 2, 3, 0, 1, 4, 2, 4 }; - std::vector vValues = { 1, 2, 3, 4, 5, 6, 7, 8 }; - return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); + std::vector vValues = { 1, 2, 3, 4, 5, 6, 7, 8 }; + return cCSRMatrix(iGlobalRowBegin, iGlobalRowEnd, iGlobalColCount, vRowPtr, vColInd, vValues); } TEST(TestOperationsCorrectness, SpMVMatchesDense) { auto A = BuildMatrix4x5(); - std::vector x = { 10, 11, 12, 13, 14 }; + std::vector x = { 10, 11, 12, 13, 14 }; // Dense reference computation iSize rows = A.iGlobalRowEnd() - A.iGlobalRowBegin(); iSize cols = A.iGlobalColCount(); // Reconstruct dense matrix row-major for reference // Using the same integer values ensures exact arithmetic for real types - std::vector dense(rows * cols, static_cast(0)); + std::vector dense(rows * cols, static_cast(0)); // Fill from pattern - dense[0 * cols + 0] = static_cast(1); - dense[0 * cols + 2] = static_cast(2); - dense[1 * cols + 3] = static_cast(3); - dense[2 * cols + 0] = static_cast(4); - dense[2 * cols + 1] = static_cast(5); - dense[2 * cols + 4] = static_cast(6); - dense[3 * cols + 2] = static_cast(7); - dense[3 * cols + 4] = static_cast(8); - std::vector y_ref(rows, static_cast(0)); + dense[0 * cols + 0] = static_cast(1); + dense[0 * cols + 2] = static_cast(2); + dense[1 * cols + 3] = static_cast(3); + dense[2 * cols + 0] = static_cast(4); + dense[2 * cols + 1] = static_cast(5); + dense[2 * cols + 4] = static_cast(6); + dense[3 * cols + 2] = static_cast(7); + dense[3 * cols + 4] = static_cast(8); + std::vector y_ref(rows, static_cast(0)); for (iSize r = 0; r < rows; ++r) { - vScalar acc = static_cast(0); + Scalar acc = static_cast(0); for (iSize c = 0; c < cols; ++c) { acc += dense[r * cols + c] * x[static_cast(c)]; } y_ref[static_cast(r)] = acc; } - std::vector y_spmv; + std::vector y_spmv; SpMV(A, x, y_spmv); ASSERT_EQ(y_spmv.size(), y_ref.size()); for (size_t i = 0; i < y_ref.size(); ++i) { @@ -65,34 +66,34 @@ TEST(TestOperationsCorrectness, SpMMMatchesDense) // col0: 1 2 3 4 5 // col1: 6 7 8 9 10 // col2: 11 12 13 14 15 - std::vector X(cols * nCols, static_cast(0)); + std::vector X(cols * nCols, static_cast(0)); for (iSize c = 0; c < cols; ++c) { - X[static_cast(c + 0 * cols)] = static_cast(1 + c); - X[static_cast(c + 1 * cols)] = static_cast(6 + c); - X[static_cast(c + 2 * cols)] = static_cast(11 + c); + X[static_cast(c + 0 * cols)] = static_cast(1 + c); + X[static_cast(c + 1 * cols)] = static_cast(6 + c); + X[static_cast(c + 2 * cols)] = static_cast(11 + c); } // Dense matrix reconstruction same as previous test - std::vector dense(rows * cols, static_cast(0)); - dense[0 * cols + 0] = static_cast(1); - dense[0 * cols + 2] = static_cast(2); - dense[1 * cols + 3] = static_cast(3); - dense[2 * cols + 0] = static_cast(4); - dense[2 * cols + 1] = static_cast(5); - dense[2 * cols + 4] = static_cast(6); - dense[3 * cols + 2] = static_cast(7); - dense[3 * cols + 4] = static_cast(8); + std::vector dense(rows * cols, static_cast(0)); + dense[0 * cols + 0] = static_cast(1); + dense[0 * cols + 2] = static_cast(2); + dense[1 * cols + 3] = static_cast(3); + dense[2 * cols + 0] = static_cast(4); + dense[2 * cols + 1] = static_cast(5); + dense[2 * cols + 4] = static_cast(6); + dense[3 * cols + 2] = static_cast(7); + dense[3 * cols + 4] = static_cast(8); // Reference multiplication Y_ref (column-major rows x nCols) - std::vector Y_ref(rows * nCols, static_cast(0)); + std::vector Y_ref(rows * nCols, static_cast(0)); for (iSize j = 0; j < nCols; ++j) { for (iSize r = 0; r < rows; ++r) { - vScalar acc = static_cast(0); + Scalar acc = static_cast(0); for (iSize c = 0; c < cols; ++c) { acc += dense[r * cols + c] * X[static_cast(c + j * cols)]; } Y_ref[static_cast(r + j * rows)] = acc; } } - std::vector Y_spmm; + std::vector Y_spmm; SpMM(A, X, nCols, Y_spmm); ASSERT_EQ(Y_spmm.size(), Y_ref.size()); for (size_t idx = 0; idx < Y_ref.size(); ++idx) { diff --git a/tests/test_symmetric.cpp b/tests/test_symmetric.cpp index befcccc..840122c 100644 --- a/tests/test_symmetric.cpp +++ b/tests/test_symmetric.cpp @@ -4,12 +4,13 @@ #include using namespace csr4mpi; +using Scalar = double; -static cCSRMatrix LoadSymLower(const std::string& rel) +static cCSRMatrix LoadSymLower(const std::string& rel) { std::string path = std::string(CSR4MPI_SOURCE_DIR) + rel; std::vector rp, ci; - std::vector vv; + std::vector vv; iIndex rows = 0, cols = 0; bool ok = LoadMatrixMarket(path, rp, ci, vv, rows, cols); EXPECT_TRUE(ok); @@ -27,22 +28,22 @@ static cCSRMatrix LoadSymLower(const std::string& rel) if (eSymm == eSymmFull) break; } - return cCSRMatrix(0, rows, cols, rp, ci, vv, eSymm); + return cCSRMatrix(0, rows, cols, rp, ci, vv, eSymm); } TEST(SymmetricTest, SpMVLowerTriangleExpansion) { auto A = LoadSymLower("/tests/data/sym3_lower.mtx"); ASSERT_TRUE(A.bIsSymmetric()); - std::vector x { static_cast(1), static_cast(2), static_cast(3) }; - std::vector y; + std::vector x { static_cast(1), static_cast(2), static_cast(3) }; + std::vector y; SpMV(A, x, y); // Full matrix equivalent: // [4 1 0; 1 3 2; 0 2 5] * [1 2 3] = [4*1+1*2=6, 1*1+3*2+2*3=1+6+6=13, 2*2+5*3=4+15=19] ASSERT_EQ(y.size(), (size_t)3); - EXPECT_EQ(y[0], static_cast(6)); - EXPECT_EQ(y[1], static_cast(13)); - EXPECT_EQ(y[2], static_cast(19)); + EXPECT_EQ(y[0], static_cast(6)); + EXPECT_EQ(y[1], static_cast(13)); + EXPECT_EQ(y[2], static_cast(19)); } TEST(SymmetricTest, SpMMLowerTriangleExpansion) @@ -50,17 +51,17 @@ TEST(SymmetricTest, SpMMLowerTriangleExpansion) auto A = LoadSymLower("/tests/data/sym3_lower.mtx"); iSize nCols = 2; // X shape 3 x 2, column-major; columns: [1 2 3]^T and [2 0 1]^T - std::vector X { static_cast(1), static_cast(2), static_cast(3), - static_cast(2), static_cast(0), static_cast(1) }; - std::vector Y; + std::vector X { static_cast(1), static_cast(2), static_cast(3), + static_cast(2), static_cast(0), static_cast(1) }; + std::vector Y; SpMM(A, X, nCols, Y); // First column: same as SpMV with x:[1 2 3] => [6,13,19] // Second column: A * [2 0 1] = [4*2 +1*0=8, 1*2+3*0+2*1=2+0+2=4, 2*0+5*1=5] ASSERT_EQ(Y.size(), (size_t)(3 * 2)); - EXPECT_EQ(Y[0], static_cast(6)); - EXPECT_EQ(Y[1], static_cast(13)); - EXPECT_EQ(Y[2], static_cast(19)); - EXPECT_EQ(Y[3], static_cast(8)); - EXPECT_EQ(Y[4], static_cast(4)); - EXPECT_EQ(Y[5], static_cast(5)); + EXPECT_EQ(Y[0], static_cast(6)); + EXPECT_EQ(Y[1], static_cast(13)); + EXPECT_EQ(Y[2], static_cast(19)); + EXPECT_EQ(Y[3], static_cast(8)); + EXPECT_EQ(Y[4], static_cast(4)); + EXPECT_EQ(Y[5], static_cast(5)); }