From 6c65da5eaa682ca5af79db099bf8b982c2fa157b Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Sun, 9 Aug 2026 10:22:39 -0700 Subject: [PATCH 1/4] [schur-complement] add the examples --- src/CMakeLists.txt | 2 + src/schur-complement-cuda/CMakeLists.txt | 8 + src/schur-complement-cuda/LICENSE | 43 +++++ src/schur-complement-cuda/Makefile | 56 +++++++ src/schur-complement-cuda/main.cu | 204 +++++++++++++++++++++++ src/schur-complement-cuda/reference.h | 123 ++++++++++++++ src/schur-complement-hip/CMakeLists.txt | 11 ++ src/schur-complement-hip/LICENSE | 43 +++++ src/schur-complement-hip/Makefile | 56 +++++++ src/schur-complement-hip/main.cu | 202 ++++++++++++++++++++++ src/schur-complement-omp/CMakeLists.txt | 11 ++ src/schur-complement-omp/LICENSE | 43 +++++ src/schur-complement-omp/Makefile | 62 +++++++ src/schur-complement-omp/Makefile.aomp | 66 ++++++++ src/schur-complement-omp/Makefile.nvc | 63 +++++++ src/schur-complement-omp/main.cpp | 181 ++++++++++++++++++++ src/schur-complement-sycl/CMakeLists.txt | 11 ++ src/schur-complement-sycl/LICENSE | 43 +++++ src/schur-complement-sycl/Makefile | 82 +++++++++ src/schur-complement-sycl/main.cpp | 180 ++++++++++++++++++++ 20 files changed, 1490 insertions(+) create mode 100644 src/schur-complement-cuda/CMakeLists.txt create mode 100644 src/schur-complement-cuda/LICENSE create mode 100644 src/schur-complement-cuda/Makefile create mode 100644 src/schur-complement-cuda/main.cu create mode 100644 src/schur-complement-cuda/reference.h create mode 100644 src/schur-complement-hip/CMakeLists.txt create mode 100644 src/schur-complement-hip/LICENSE create mode 100644 src/schur-complement-hip/Makefile create mode 100644 src/schur-complement-hip/main.cu create mode 100644 src/schur-complement-omp/CMakeLists.txt create mode 100644 src/schur-complement-omp/LICENSE create mode 100644 src/schur-complement-omp/Makefile create mode 100644 src/schur-complement-omp/Makefile.aomp create mode 100644 src/schur-complement-omp/Makefile.nvc create mode 100644 src/schur-complement-omp/main.cpp create mode 100644 src/schur-complement-sycl/CMakeLists.txt create mode 100644 src/schur-complement-sycl/LICENSE create mode 100644 src/schur-complement-sycl/Makefile create mode 100644 src/schur-complement-sycl/main.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a1fcf804b..374c7e71a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -370,6 +370,7 @@ set(HECBENCH_POC_BENCHMARKS scatterAdd scatterThrust scel + schur-complement score sddmm-batch seam-carving @@ -391,6 +392,7 @@ set(HECBENCH_POC_BENCHMARKS sortKV sosfil sparkler + sparseMLA spaxpby spd2s spgeam diff --git a/src/schur-complement-cuda/CMakeLists.txt b/src/schur-complement-cuda/CMakeLists.txt new file mode 100644 index 000000000..e2fd29122 --- /dev/null +++ b/src/schur-complement-cuda/CMakeLists.txt @@ -0,0 +1,8 @@ +# schur-complement-cuda/CMakeLists.txt + +add_hecbench_benchmark( + NAME schur-complement + MODEL cuda + SOURCES main.cu + CATEGORIES math +) diff --git a/src/schur-complement-cuda/LICENSE b/src/schur-complement-cuda/LICENSE new file mode 100644 index 000000000..78a030c8d --- /dev/null +++ b/src/schur-complement-cuda/LICENSE @@ -0,0 +1,43 @@ +BSD 3-Clause License. Please also read “Additional BSD Notice” below. + +Copyright (c) 2017, Lawrence Livermore National Laboratory +LLNL-CODE-742473.All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: +i. Redistributions of source code must retain the above copyright notice, this list +of conditions and the disclaimer below. +ii. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the disclaimer (as noted below) in the documentation and/or +other materials provided with the distribution. +iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +endorse or promote products derived from this software without specific prior written +permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Additional BSD Notice +1. This notice is required to be provided under our contract with the U.S. Department +of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +Contract No. DE-AC52-07NA27344 with the DOE. +2. Neither the United States Government nor Lawrence Livermore National Security, LLC +nor any of their employees, makes any warranty, express or implied, or assumes any +liability or responsibility for the accuracy, completeness, or usefulness of any +information, apparatus, product, or process disclosed, or represents that its use would +not infringe privately-owned rights. +3. Also, reference herein to any specific commercial products, process, or services by +trade name, trademark, manufacturer or otherwise does not necessarily constitute or +imply its endorsement, recommendation, or favoring by the United States Government or +Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +herein do not necessarily state or reflect those of the United States Government or +Lawrence Livermore National Security, LLC, and shall not be used for advertising or +product endorsement purposes. diff --git a/src/schur-complement-cuda/Makefile b/src/schur-complement-cuda/Makefile new file mode 100644 index 000000000..df77c9c1c --- /dev/null +++ b/src/schur-complement-cuda/Makefile @@ -0,0 +1,56 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = nvcc +OPTIMIZE = yes +DEBUG = no +ARCH = sm_60 +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cu + +obj = $(source:.cu=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +# Standard Flags +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall -arch=$(ARCH) + +# Linker Flags +LDFLAGS = + +# Debug Flags +ifeq ($(DEBUG),yes) + CFLAGS += -g -DDEBUG + LDFLAGS += -g +endif + +# Optimization Flags +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) Makefile + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cu reference.h Makefile + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) + +run: $(program) + $(LAUNCHER) ./$(program) 4096 32 100 diff --git a/src/schur-complement-cuda/main.cu b/src/schur-complement-cuda/main.cu new file mode 100644 index 000000000..953f165c8 --- /dev/null +++ b/src/schur-complement-cuda/main.cu @@ -0,0 +1,204 @@ +#include +#include +#include +#include +#include +#include +#include "reference.h" + +#define CHECK(call) \ + do { \ + const cudaError_t err = (call); \ + if (err != cudaSuccess) { \ + fprintf(stderr, "CUDA error %s:%d '%s': %s\n", __FILE__, __LINE__, \ + #call, cudaGetErrorString(err)); \ + exit(EXIT_FAILURE); \ + } \ + } while (0) + +// Original HiOp mapping is one thread per row i, looping over all partner rows +// j -- only O(m) work items, which leaves the GPU underused. Here each thread +// instead owns a single (i,j) row pair and performs one merge, raising the +// parallelism to O(m^2). Threads are laid out 2D with the fast (x) dimension on +// j so the writes to W[i][j] are coalesced; BLOCK_X is warp-aligned. +#define BLOCK_X 32 +#define BLOCK_Y 8 + +// W += alpha * J * D^{-1} * J^T : thread (i,j) computes one upper-triangle entry +// (j >= i) of the symmetric diagonal block. +__global__ void mdinvmtrans_diag(int nrows, + const int* __restrict__ row_start, + const int* __restrict__ jcol, + const double* __restrict__ values, + const double* __restrict__ D, + int row_dest_start, int col_dest_start, + double alpha, double* __restrict__ W, int m_W) +{ + const int j = blockIdx.x * blockDim.x + threadIdx.x; + const int i = blockIdx.y * blockDim.y + threadIdx.y; + if (i >= nrows || j >= nrows || j < i) return; + + double acc = 0.0; + if (i == j) { + for (int k = row_start[i]; k < row_start[i + 1]; k++) + acc += values[k] / D[jcol[k]] * values[k]; + } else { + int ki = row_start[i], kj = row_start[j]; + const int kie = row_start[i + 1], kje = row_start[j + 1]; + while (ki < kie && kj < kje) { + const int ci = jcol[ki], cj = jcol[kj]; + if (ci == cj) { acc += values[ki] / D[ci] * values[kj]; ki++; kj++; } + else if (ci < cj) ki++; + else kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; +} + +// W += alpha * J1 * D^{-1} * J2^T : thread (i,j) computes one entry of the full +// off-diagonal block by merging row i of J1 against row j of J2. +__global__ void mdinvntrans(int m1, int m2, + const int* __restrict__ rs1, const int* __restrict__ jc1, + const double* __restrict__ v1, + const int* __restrict__ rs2, const int* __restrict__ jc2, + const double* __restrict__ v2, + const double* __restrict__ D, + int row_dest_start, int col_dest_start, + double alpha, double* __restrict__ W, int m_W) +{ + const int j = blockIdx.x * blockDim.x + threadIdx.x; + const int i = blockIdx.y * blockDim.y + threadIdx.y; + if (i >= m1 || j >= m2) return; + + double acc = 0.0; + int ki = rs1[i], kj = rs2[j]; + const int kie = rs1[i + 1], kje = rs2[j + 1]; + while (ki < kie && kj < kje) { + const int ci = jc1[ki], cj = jc2[kj]; + if (ci == cj) { acc += v1[ki] / D[ci] * v2[kj]; ki++; kj++; } + else if (ci < cj) ki++; + else kj++; + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; +} + +int main(int argc, char* argv[]) +{ + if (argc != 4) { + printf("Usage: %s \n", argv[0]); + return 1; + } + + const int m = atoi(argv[1]); + const int nnz_row = atoi(argv[2]); + const int repeat = atoi(argv[3]); + + // number of variables (columns of J); a few times wider than nnz_row so the + // per-row column lists overlap only partially, exercising the merge loop + const int nx = 8 * nnz_row + 1024; + const double alpha = -1.0; // HiOp assembles the Schur complement with -1 + + std::vector h_rs1, h_jc1, h_rs2, h_jc2; + std::vector h_v1, h_v2, h_D; + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + + const int nnz = (int)h_v1.size(); + const size_t w_elems = (size_t)m * m; + const size_t w_bytes = w_elems * sizeof(double); + + int *d_rs1, *d_jc1, *d_rs2, *d_jc2; + double *d_v1, *d_v2, *d_D, *d_W; + CHECK(cudaMalloc((void**)&d_rs1, (m + 1) * sizeof(int))); + CHECK(cudaMalloc((void**)&d_jc1, nnz * sizeof(int))); + CHECK(cudaMalloc((void**)&d_v1, nnz * sizeof(double))); + CHECK(cudaMalloc((void**)&d_rs2, (m + 1) * sizeof(int))); + CHECK(cudaMalloc((void**)&d_jc2, nnz * sizeof(int))); + CHECK(cudaMalloc((void**)&d_v2, nnz * sizeof(double))); + CHECK(cudaMalloc((void**)&d_D, nx * sizeof(double))); + CHECK(cudaMalloc((void**)&d_W, w_bytes)); + + CHECK(cudaMemcpy(d_rs1, h_rs1.data(), (m + 1) * sizeof(int), cudaMemcpyHostToDevice)); + CHECK(cudaMemcpy(d_jc1, h_jc1.data(), nnz * sizeof(int), cudaMemcpyHostToDevice)); + CHECK(cudaMemcpy(d_v1, h_v1.data(), nnz * sizeof(double), cudaMemcpyHostToDevice)); + CHECK(cudaMemcpy(d_rs2, h_rs2.data(), (m + 1) * sizeof(int), cudaMemcpyHostToDevice)); + CHECK(cudaMemcpy(d_jc2, h_jc2.data(), nnz * sizeof(int), cudaMemcpyHostToDevice)); + CHECK(cudaMemcpy(d_v2, h_v2.data(), nnz * sizeof(double), cudaMemcpyHostToDevice)); + CHECK(cudaMemcpy(d_D, h_D.data(), nx * sizeof(double), cudaMemcpyHostToDevice)); + + const dim3 block(BLOCK_X, BLOCK_Y); + const dim3 grid((m + BLOCK_X - 1) / BLOCK_X, (m + BLOCK_Y - 1) / BLOCK_Y); + + std::vector h_W(w_elems); + std::vector h_ref(w_elems); + int errors = 0; + + // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- + + // host/device correctness check (run once, verify against reference) before timing + CHECK(cudaMemset(d_W, 0, w_bytes)); + mdinvmtrans_diag<<>>(m, d_rs1, d_jc1, d_v1, d_D, 0, 0, alpha, d_W, m); + CHECK(cudaMemcpy(h_W.data(), d_W, w_bytes, cudaMemcpyDeviceToHost)); + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvmtrans_diag(m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (int i = 0; i < m; i++) + for (int j = i; j < m; j++) + if (!close_enough(h_W[(size_t)i * m + j], h_ref[(size_t)i * m + j], 1e-10)) { + errors++; i = m; break; + } + + // benchmark + CHECK(cudaMemset(d_W, 0, w_bytes)); + CHECK(cudaDeviceSynchronize()); + auto start = std::chrono::steady_clock::now(); + + for (int i = 0; i < repeat; i++) + mdinvmtrans_diag<<>>( + m, d_rs1, d_jc1, d_v1, d_D, 0, 0, alpha, d_W, m); + + CHECK(cudaDeviceSynchronize()); + auto end = std::chrono::steady_clock::now(); + auto time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvMtrans (diag block): %f (us)\n", + time * 1e-3 / repeat); + + // --- off-diagonal block: W += alpha * J1 D^{-1} J2^T --------------------- + + // host/device correctness check before timing + CHECK(cudaMemset(d_W, 0, w_bytes)); + mdinvntrans<<>>(m, m, d_rs1, d_jc1, d_v1, d_rs2, d_jc2, d_v2, d_D, 0, 0, alpha, d_W, m); + CHECK(cudaMemcpy(h_W.data(), d_W, w_bytes, cudaMemcpyDeviceToHost)); + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvntrans(m, m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), + h_rs2.data(), h_jc2.data(), h_v2.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (size_t k = 0; k < w_elems; k++) + if (!close_enough(h_W[k], h_ref[k], 1e-10)) { errors++; break; } + + // benchmark + CHECK(cudaMemset(d_W, 0, w_bytes)); + CHECK(cudaDeviceSynchronize()); + start = std::chrono::steady_clock::now(); + + for (int i = 0; i < repeat; i++) + mdinvntrans<<>>( + m, m, d_rs1, d_jc1, d_v1, d_rs2, d_jc2, d_v2, d_D, 0, 0, alpha, d_W, m); + + CHECK(cudaDeviceSynchronize()); + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvNtrans (off-diag): %f (us)\n", + time * 1e-3 / repeat); + + printf("%s\n", errors == 0 ? "PASS" : "FAIL"); + + CHECK(cudaFree(d_rs1)); CHECK(cudaFree(d_jc1)); CHECK(cudaFree(d_v1)); + CHECK(cudaFree(d_rs2)); CHECK(cudaFree(d_jc2)); CHECK(cudaFree(d_v2)); + CHECK(cudaFree(d_D)); CHECK(cudaFree(d_W)); + + return 0; +} diff --git a/src/schur-complement-cuda/reference.h b/src/schur-complement-cuda/reference.h new file mode 100644 index 000000000..b55fde3fb --- /dev/null +++ b/src/schur-complement-cuda/reference.h @@ -0,0 +1,123 @@ +#ifndef REFERENCE_H +#define REFERENCE_H + +#include +#include +#include +#include + +// https://github.com/ORNL/hiop +// +// addMDinvMtransToDiagBlockOfSymDeMatUTri : W += alpha * J * D^{-1} * J^T +// addMDinvNtransToSymDeMatUTri : W += alpha * J1 * D^{-1} * J2^T +// +// J is a sparse matrix stored CSR-style with per-row sorted column indices, +// D is a diagonal (stored as a vector), and W is a dense matrix into whose +// upper triangle the (symmetric) result is scattered. + +// Build a CSR matrix with `m` rows, `nx` columns and exactly `nnz_row` +// nonzeros per row (column indices sorted ascending within each row). +static void gen_csr(int m, int nx, int nnz_row, unsigned seed, + std::vector& row_start, + std::vector& jcol, + std::vector& values) +{ + if (nnz_row > nx) nnz_row = nx; + srand(seed); + + row_start.resize(m + 1); + jcol.resize((size_t)m * nnz_row); + values.resize((size_t)m * nnz_row); + + std::vector cols(nx); + for (int c = 0; c < nx; c++) cols[c] = c; + + row_start[0] = 0; + for (int i = 0; i < m; i++) { + // partial Fisher-Yates: pick the first nnz_row distinct columns + for (int k = 0; k < nnz_row; k++) { + int r = k + rand() % (nx - k); + std::swap(cols[k], cols[r]); + } + std::sort(cols.begin(), cols.begin() + nnz_row); + const size_t base = (size_t)i * nnz_row; + for (int k = 0; k < nnz_row; k++) { + jcol[base + k] = cols[k]; + values[base + k] = 2.0 * (rand() / (double)RAND_MAX) - 1.0; + } + row_start[i + 1] = row_start[i] + nnz_row; + } +} + +// D is the diagonal; keep entries well away from zero so 1/D is well behaved. +static void gen_diag(int nx, unsigned seed, std::vector& D) +{ + srand(seed); + D.resize(nx); + for (int c = 0; c < nx; c++) D[c] = 0.1 + rand() / (double)RAND_MAX; +} + +// W += alpha * J * D^{-1} * J^T (upper triangle of the m x m diagonal block) +static void reference_mdinvmtrans_diag(int m, int m_W, + const int* row_start, const int* jcol, + const double* values, const double* D, + int row_dest_start, int col_dest_start, + double alpha, double* W) +{ + for (int i = 0; i < m; i++) { + double acc = 0.0; + for (int k = row_start[i]; k < row_start[i + 1]; k++) + acc += values[k] / D[jcol[k]] * values[k]; + W[(size_t)(i + row_dest_start) * m_W + i + col_dest_start] += alpha * acc; + + for (int j = i + 1; j < m; j++) { + acc = 0.0; + int ki = row_start[i], kj = row_start[j]; + while (ki < row_start[i + 1] && kj < row_start[j + 1]) { + if (jcol[ki] == jcol[kj]) { + acc += values[ki] / D[jcol[ki]] * values[kj]; + ki++; kj++; + } else if (jcol[ki] < jcol[kj]) { + ki++; + } else { + kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + } + } +} + +// W += alpha * J1 * D^{-1} * J2^T (full m1 x m2 off-diagonal block) +static void reference_mdinvntrans(int m1, int m2, int m_W, + const int* rs1, const int* jc1, const double* v1, + const int* rs2, const int* jc2, const double* v2, + const double* D, + int row_dest_start, int col_dest_start, + double alpha, double* W) +{ + for (int i = 0; i < m1; i++) { + for (int j = 0; j < m2; j++) { + double acc = 0.0; + int ki = rs1[i], kj = rs2[j]; + while (ki < rs1[i + 1] && kj < rs2[j + 1]) { + if (jc1[ki] == jc2[kj]) { + acc += v1[ki] / D[jc1[ki]] * v2[kj]; + ki++; kj++; + } else if (jc1[ki] < jc2[kj]) { + ki++; + } else { + kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + } + } +} + +static bool close_enough(double a, double b, double tol) +{ + return fabs(a - b) <= tol * (1.0 + fabs(b)); +} + +#endif diff --git a/src/schur-complement-hip/CMakeLists.txt b/src/schur-complement-hip/CMakeLists.txt new file mode 100644 index 000000000..072306a19 --- /dev/null +++ b/src/schur-complement-hip/CMakeLists.txt @@ -0,0 +1,11 @@ +# schur-complement-hip/CMakeLists.txt + +set(INC_DIR "${CMAKE_CURRENT_LIST_DIR}/../schur-complement-cuda") + +add_hecbench_benchmark( + NAME schur-complement + MODEL hip + SOURCES main.cu + CATEGORIES math + INCLUDE_DIRS ${INC_DIR} +) diff --git a/src/schur-complement-hip/LICENSE b/src/schur-complement-hip/LICENSE new file mode 100644 index 000000000..78a030c8d --- /dev/null +++ b/src/schur-complement-hip/LICENSE @@ -0,0 +1,43 @@ +BSD 3-Clause License. Please also read “Additional BSD Notice” below. + +Copyright (c) 2017, Lawrence Livermore National Laboratory +LLNL-CODE-742473.All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: +i. Redistributions of source code must retain the above copyright notice, this list +of conditions and the disclaimer below. +ii. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the disclaimer (as noted below) in the documentation and/or +other materials provided with the distribution. +iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +endorse or promote products derived from this software without specific prior written +permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Additional BSD Notice +1. This notice is required to be provided under our contract with the U.S. Department +of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +Contract No. DE-AC52-07NA27344 with the DOE. +2. Neither the United States Government nor Lawrence Livermore National Security, LLC +nor any of their employees, makes any warranty, express or implied, or assumes any +liability or responsibility for the accuracy, completeness, or usefulness of any +information, apparatus, product, or process disclosed, or represents that its use would +not infringe privately-owned rights. +3. Also, reference herein to any specific commercial products, process, or services by +trade name, trademark, manufacturer or otherwise does not necessarily constitute or +imply its endorsement, recommendation, or favoring by the United States Government or +Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +herein do not necessarily state or reflect those of the United States Government or +Lawrence Livermore National Security, LLC, and shall not be used for advertising or +product endorsement purposes. diff --git a/src/schur-complement-hip/Makefile b/src/schur-complement-hip/Makefile new file mode 100644 index 000000000..68e4db52f --- /dev/null +++ b/src/schur-complement-hip/Makefile @@ -0,0 +1,56 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = hipcc +OPTIMIZE = yes +DEBUG = no +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cu + +obj = $(source:.cu=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +# Standard Flags +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../schur-complement-cuda + +# Linker Flags +LDFLAGS = + +# Debug Flags +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +# Optimization Flags +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cu ../schur-complement-cuda/reference.h + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) + +run: $(program) + $(LAUNCHER) ./$(program) 4096 32 100 diff --git a/src/schur-complement-hip/main.cu b/src/schur-complement-hip/main.cu new file mode 100644 index 000000000..db0b1e0b8 --- /dev/null +++ b/src/schur-complement-hip/main.cu @@ -0,0 +1,202 @@ +#include +#include +#include +#include +#include +#include +#include "reference.h" + +#define CHECK(call) \ + do { \ + const hipError_t err = (call); \ + if (err != hipSuccess) { \ + fprintf(stderr, "HIP error %s:%d '%s': %s\n", __FILE__, __LINE__, \ + #call, hipGetErrorString(err)); \ + exit(EXIT_FAILURE); \ + } \ + } while (0) + +// Original HiOp mapping is one thread per row i, looping over all partner rows +// j -- only O(m) work items, which leaves the GPU underused. Here each thread +// instead owns a single (i,j) row pair and performs one merge, raising the +// parallelism to O(m^2). Threads are laid out 2D with the fast (x) dimension on +// j so the writes to W[i][j] are coalesced; BLOCK_X is warp/wavefront-aligned. +#define BLOCK_X 32 +#define BLOCK_Y 8 + +// W += alpha * J * D^{-1} * J^T : thread (i,j) computes one upper-triangle entry +// (j >= i) of the symmetric diagonal block. +__global__ void mdinvmtrans_diag(int nrows, + const int* __restrict__ row_start, + const int* __restrict__ jcol, + const double* __restrict__ values, + const double* __restrict__ D, + int row_dest_start, int col_dest_start, + double alpha, double* __restrict__ W, int m_W) +{ + const int j = blockIdx.x * blockDim.x + threadIdx.x; + const int i = blockIdx.y * blockDim.y + threadIdx.y; + if (i >= nrows || j >= nrows || j < i) return; + + double acc = 0.0; + if (i == j) { + for (int k = row_start[i]; k < row_start[i + 1]; k++) + acc += values[k] / D[jcol[k]] * values[k]; + } else { + int ki = row_start[i], kj = row_start[j]; + const int kie = row_start[i + 1], kje = row_start[j + 1]; + while (ki < kie && kj < kje) { + const int ci = jcol[ki], cj = jcol[kj]; + if (ci == cj) { acc += values[ki] / D[ci] * values[kj]; ki++; kj++; } + else if (ci < cj) ki++; + else kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; +} + +// W += alpha * J1 * D^{-1} * J2^T : thread (i,j) computes one entry of the full +// off-diagonal block by merging row i of J1 against row j of J2. +__global__ void mdinvntrans(int m1, int m2, + const int* __restrict__ rs1, const int* __restrict__ jc1, + const double* __restrict__ v1, + const int* __restrict__ rs2, const int* __restrict__ jc2, + const double* __restrict__ v2, + const double* __restrict__ D, + int row_dest_start, int col_dest_start, + double alpha, double* __restrict__ W, int m_W) +{ + const int j = blockIdx.x * blockDim.x + threadIdx.x; + const int i = blockIdx.y * blockDim.y + threadIdx.y; + if (i >= m1 || j >= m2) return; + + double acc = 0.0; + int ki = rs1[i], kj = rs2[j]; + const int kie = rs1[i + 1], kje = rs2[j + 1]; + while (ki < kie && kj < kje) { + const int ci = jc1[ki], cj = jc2[kj]; + if (ci == cj) { acc += v1[ki] / D[ci] * v2[kj]; ki++; kj++; } + else if (ci < cj) ki++; + else kj++; + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; +} + +int main(int argc, char* argv[]) +{ + if (argc != 4) { + printf("Usage: %s \n", argv[0]); + return 1; + } + + const int m = atoi(argv[1]); + const int nnz_row = atoi(argv[2]); + const int repeat = atoi(argv[3]); + + const int nx = 8 * nnz_row + 1024; + const double alpha = -1.0; + + std::vector h_rs1, h_jc1, h_rs2, h_jc2; + std::vector h_v1, h_v2, h_D; + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + + const int nnz = (int)h_v1.size(); + const size_t w_elems = (size_t)m * m; + const size_t w_bytes = w_elems * sizeof(double); + + int *d_rs1, *d_jc1, *d_rs2, *d_jc2; + double *d_v1, *d_v2, *d_D, *d_W; + CHECK(hipMalloc((void**)&d_rs1, (m + 1) * sizeof(int))); + CHECK(hipMalloc((void**)&d_jc1, nnz * sizeof(int))); + CHECK(hipMalloc((void**)&d_v1, nnz * sizeof(double))); + CHECK(hipMalloc((void**)&d_rs2, (m + 1) * sizeof(int))); + CHECK(hipMalloc((void**)&d_jc2, nnz * sizeof(int))); + CHECK(hipMalloc((void**)&d_v2, nnz * sizeof(double))); + CHECK(hipMalloc((void**)&d_D, nx * sizeof(double))); + CHECK(hipMalloc((void**)&d_W, w_bytes)); + + CHECK(hipMemcpy(d_rs1, h_rs1.data(), (m + 1) * sizeof(int), hipMemcpyHostToDevice)); + CHECK(hipMemcpy(d_jc1, h_jc1.data(), nnz * sizeof(int), hipMemcpyHostToDevice)); + CHECK(hipMemcpy(d_v1, h_v1.data(), nnz * sizeof(double), hipMemcpyHostToDevice)); + CHECK(hipMemcpy(d_rs2, h_rs2.data(), (m + 1) * sizeof(int), hipMemcpyHostToDevice)); + CHECK(hipMemcpy(d_jc2, h_jc2.data(), nnz * sizeof(int), hipMemcpyHostToDevice)); + CHECK(hipMemcpy(d_v2, h_v2.data(), nnz * sizeof(double), hipMemcpyHostToDevice)); + CHECK(hipMemcpy(d_D, h_D.data(), nx * sizeof(double), hipMemcpyHostToDevice)); + + const dim3 block(BLOCK_X, BLOCK_Y); + const dim3 grid((m + BLOCK_X - 1) / BLOCK_X, (m + BLOCK_Y - 1) / BLOCK_Y); + + std::vector h_W(w_elems); + std::vector h_ref(w_elems); + int errors = 0; + + // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- + + // host/device correctness check (run once, verify against reference) before timing + CHECK(hipMemset(d_W, 0, w_bytes)); + mdinvmtrans_diag<<>>(m, d_rs1, d_jc1, d_v1, d_D, 0, 0, alpha, d_W, m); + CHECK(hipMemcpy(h_W.data(), d_W, w_bytes, hipMemcpyDeviceToHost)); + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvmtrans_diag(m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (int i = 0; i < m; i++) + for (int j = i; j < m; j++) + if (!close_enough(h_W[(size_t)i * m + j], h_ref[(size_t)i * m + j], 1e-10)) { + errors++; i = m; break; + } + + // benchmark + CHECK(hipMemset(d_W, 0, w_bytes)); + CHECK(hipDeviceSynchronize()); + auto start = std::chrono::steady_clock::now(); + + for (int i = 0; i < repeat; i++) + mdinvmtrans_diag<<>>( + m, d_rs1, d_jc1, d_v1, d_D, 0, 0, alpha, d_W, m); + + CHECK(hipDeviceSynchronize()); + auto end = std::chrono::steady_clock::now(); + auto time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvMtrans (diag block): %f (us)\n", + time * 1e-3 / repeat); + + // --- off-diagonal block: W += alpha * J1 D^{-1} J2^T --------------------- + + // host/device correctness check before timing + CHECK(hipMemset(d_W, 0, w_bytes)); + mdinvntrans<<>>(m, m, d_rs1, d_jc1, d_v1, d_rs2, d_jc2, d_v2, d_D, 0, 0, alpha, d_W, m); + CHECK(hipMemcpy(h_W.data(), d_W, w_bytes, hipMemcpyDeviceToHost)); + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvntrans(m, m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), + h_rs2.data(), h_jc2.data(), h_v2.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (size_t k = 0; k < w_elems; k++) + if (!close_enough(h_W[k], h_ref[k], 1e-10)) { errors++; break; } + + // benchmark + CHECK(hipMemset(d_W, 0, w_bytes)); + CHECK(hipDeviceSynchronize()); + start = std::chrono::steady_clock::now(); + + for (int i = 0; i < repeat; i++) + mdinvntrans<<>>( + m, m, d_rs1, d_jc1, d_v1, d_rs2, d_jc2, d_v2, d_D, 0, 0, alpha, d_W, m); + + CHECK(hipDeviceSynchronize()); + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvNtrans (off-diag): %f (us)\n", + time * 1e-3 / repeat); + + printf("%s\n", errors == 0 ? "PASS" : "FAIL"); + + CHECK(hipFree(d_rs1)); CHECK(hipFree(d_jc1)); CHECK(hipFree(d_v1)); + CHECK(hipFree(d_rs2)); CHECK(hipFree(d_jc2)); CHECK(hipFree(d_v2)); + CHECK(hipFree(d_D)); CHECK(hipFree(d_W)); + + return 0; +} diff --git a/src/schur-complement-omp/CMakeLists.txt b/src/schur-complement-omp/CMakeLists.txt new file mode 100644 index 000000000..72d9eab75 --- /dev/null +++ b/src/schur-complement-omp/CMakeLists.txt @@ -0,0 +1,11 @@ +# schur-complement-omp/CMakeLists.txt + +set(INC_DIR "${CMAKE_CURRENT_LIST_DIR}/../schur-complement-cuda") + +add_hecbench_benchmark( + NAME schur-complement + MODEL omp + SOURCES main.cpp + CATEGORIES math + INCLUDE_DIRS ${INC_DIR} +) diff --git a/src/schur-complement-omp/LICENSE b/src/schur-complement-omp/LICENSE new file mode 100644 index 000000000..78a030c8d --- /dev/null +++ b/src/schur-complement-omp/LICENSE @@ -0,0 +1,43 @@ +BSD 3-Clause License. Please also read “Additional BSD Notice” below. + +Copyright (c) 2017, Lawrence Livermore National Laboratory +LLNL-CODE-742473.All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: +i. Redistributions of source code must retain the above copyright notice, this list +of conditions and the disclaimer below. +ii. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the disclaimer (as noted below) in the documentation and/or +other materials provided with the distribution. +iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +endorse or promote products derived from this software without specific prior written +permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Additional BSD Notice +1. This notice is required to be provided under our contract with the U.S. Department +of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +Contract No. DE-AC52-07NA27344 with the DOE. +2. Neither the United States Government nor Lawrence Livermore National Security, LLC +nor any of their employees, makes any warranty, express or implied, or assumes any +liability or responsibility for the accuracy, completeness, or usefulness of any +information, apparatus, product, or process disclosed, or represents that its use would +not infringe privately-owned rights. +3. Also, reference herein to any specific commercial products, process, or services by +trade name, trademark, manufacturer or otherwise does not necessarily constitute or +imply its endorsement, recommendation, or favoring by the United States Government or +Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +herein do not necessarily state or reflect those of the United States Government or +Lawrence Livermore National Security, LLC, and shall not be used for advertising or +product endorsement purposes. diff --git a/src/schur-complement-omp/Makefile b/src/schur-complement-omp/Makefile new file mode 100644 index 000000000..098af89bc --- /dev/null +++ b/src/schur-complement-omp/Makefile @@ -0,0 +1,62 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = icpx +OPTIMIZE = yes +DEBUG = no +DEVICE = gpu +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +# Standard Flags +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../schur-complement-cuda + +# Linker Flags +LDFLAGS = + +# Debug Flags +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +# Optimization Flags +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(DEVICE),gpu) + CFLAGS +=-fiopenmp -fopenmp-targets=spir64 -D__STRICT_ANSI__ +else + CFLAGS +=-qopenmp +endif +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp ../schur-complement-cuda/reference.h Makefile + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) + +run: $(program) + $(LAUNCHER) ./$(program) 4096 32 100 diff --git a/src/schur-complement-omp/Makefile.aomp b/src/schur-complement-omp/Makefile.aomp new file mode 100644 index 000000000..2d721b386 --- /dev/null +++ b/src/schur-complement-omp/Makefile.aomp @@ -0,0 +1,66 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = clang++ +OPTIMIZE = yes +DEBUG = no +DEVICE = gpu +ARCH = gfx906 +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +# Standard Flags +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../schur-complement-cuda + +# Linker Flags +LDFLAGS = + +# Debug Flags +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +# Optimization Flags +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(DEVICE),gpu) + CFLAGS += -target x86_64-pc-linux-gnu \ + -fopenmp -fopenmp-targets=amdgcn-amd-amdhsa \ + -Xopenmp-target=amdgcn-amd-amdhsa \ + -march=$(ARCH) +else + CFLAGS +=-fopenmp +endif +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp ../schur-complement-cuda/reference.h Makefile.aomp + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) + +run: $(program) + $(LAUNCHER) ./$(program) 4096 32 100 diff --git a/src/schur-complement-omp/Makefile.nvc b/src/schur-complement-omp/Makefile.nvc new file mode 100644 index 000000000..13ed482e4 --- /dev/null +++ b/src/schur-complement-omp/Makefile.nvc @@ -0,0 +1,63 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = nvc++ +OPTIMIZE = yes +DEBUG = no +DEVICE = gpu +SM = cc70 +LAUNCHER = + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +# Standard Flags +CFLAGS := $(EXTRA_CFLAGS) -Wall -I../schur-complement-cuda + +# Linker Flags +LDFLAGS = + +# Debug Flags +ifeq ($(DEBUG),yes) + CFLAGS += -g + LDFLAGS += -g +endif + +# Optimization Flags +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(DEVICE),gpu) + CFLAGS +=-Minfo -mp=gpu -gpu=$(SM) +else + CFLAGS += +endif +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp ../schur-complement-cuda/reference.h Makefile.nvc + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) + +run: $(program) + $(LAUNCHER) ./$(program) 4096 32 100 diff --git a/src/schur-complement-omp/main.cpp b/src/schur-complement-omp/main.cpp new file mode 100644 index 000000000..63806e901 --- /dev/null +++ b/src/schur-complement-omp/main.cpp @@ -0,0 +1,181 @@ +#include +#include +#include +#include +#include +#include +#include "reference.h" + +int main(int argc, char* argv[]) +{ + if (argc != 4) { + printf("Usage: %s \n", argv[0]); + return 1; + } + + const int m = atoi(argv[1]); + const int nnz_row = atoi(argv[2]); + const int repeat = atoi(argv[3]); + + const int nx = 8 * nnz_row + 1024; + const double alpha = -1.0; + + std::vector h_rs1, h_jc1, h_rs2, h_jc2; + std::vector h_v1, h_v2, h_D; + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + + const int nnz = (int)h_v1.size(); + const size_t w_elems = (size_t)m * m; + + const int m_W = m; + const int row_dest_start = 0, col_dest_start = 0; + + int* rs1 = h_rs1.data(); int* jc1 = h_jc1.data(); double* v1 = h_v1.data(); + int* rs2 = h_rs2.data(); int* jc2 = h_jc2.data(); double* v2 = h_v2.data(); + double* D = h_D.data(); + std::vector h_W(w_elems, 0.0); + double* W = h_W.data(); + std::vector h_ref(w_elems); + int errors = 0; + + #pragma omp target enter data map(to: rs1[0:m+1], jc1[0:nnz], v1[0:nnz], \ + rs2[0:m+1], jc2[0:nnz], v2[0:nnz], \ + D[0:nx], W[0:w_elems]) + + // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- + + // host/device correctness check (run once, verify against reference) before timing + #pragma omp target teams distribute parallel for thread_limit(128) + for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; + + #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + for (int i = 0; i < m; i++) { + for (int j = 0; j < m; j++) { + if (j < i) continue; + double acc = 0.0; + if (i == j) { + for (int k = rs1[i]; k < rs1[i + 1]; k++) + acc += v1[k] / D[jc1[k]] * v1[k]; + } else { + int ki = rs1[i], kj = rs1[j]; + while (ki < rs1[i + 1] && kj < rs1[j + 1]) { + if (jc1[ki] == jc1[kj]) { acc += v1[ki] / D[jc1[ki]] * v1[kj]; ki++; kj++; } + else if (jc1[ki] < jc1[kj]) ki++; else kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + } + } + #pragma omp target update from(W[0:w_elems]) + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvmtrans_diag(m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (int i = 0; i < m; i++) + for (int j = i; j < m; j++) + if (!close_enough(h_W[(size_t)i * m + j], h_ref[(size_t)i * m + j], 1e-10)) { + errors++; i = m; break; + } + + // benchmark + #pragma omp target teams distribute parallel for thread_limit(128) + for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; + + auto start = std::chrono::steady_clock::now(); + + for (int r = 0; r < repeat; r++) { + // one iteration per (i,j) row pair (collapsed) instead of per row i + #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + for (int i = 0; i < m; i++) { + for (int j = 0; j < m; j++) { + if (j < i) continue; + double acc = 0.0; + if (i == j) { + for (int k = rs1[i]; k < rs1[i + 1]; k++) + acc += v1[k] / D[jc1[k]] * v1[k]; + } else { + int ki = rs1[i], kj = rs1[j]; + while (ki < rs1[i + 1] && kj < rs1[j + 1]) { + if (jc1[ki] == jc1[kj]) { acc += v1[ki] / D[jc1[ki]] * v1[kj]; ki++; kj++; } + else if (jc1[ki] < jc1[kj]) ki++; else kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + } + } + } + + auto end = std::chrono::steady_clock::now(); + auto time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvMtrans (diag block): %f (us)\n", + time * 1e-3 / repeat); + + // --- off-diagonal block: W += alpha * J1 D^{-1} J2^T --------------------- + + // host/device correctness check before timing + #pragma omp target teams distribute parallel for thread_limit(128) + for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; + + #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + for (int i = 0; i < m; i++) { + for (int j = 0; j < m; j++) { + double acc = 0.0; + int ki = rs1[i], kj = rs2[j]; + while (ki < rs1[i + 1] && kj < rs2[j + 1]) { + if (jc1[ki] == jc2[kj]) { acc += v1[ki] / D[jc1[ki]] * v2[kj]; ki++; kj++; } + else if (jc1[ki] < jc2[kj]) ki++; else kj++; + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + } + } + #pragma omp target update from(W[0:w_elems]) + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvntrans(m, m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), + h_rs2.data(), h_jc2.data(), h_v2.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (size_t k = 0; k < w_elems; k++) + if (!close_enough(h_W[k], h_ref[k], 1e-10)) { errors++; break; } + + // benchmark + #pragma omp target teams distribute parallel for thread_limit(128) + for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; + + start = std::chrono::steady_clock::now(); + + for (int r = 0; r < repeat; r++) { + #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + for (int i = 0; i < m; i++) { + for (int j = 0; j < m; j++) { + double acc = 0.0; + int ki = rs1[i], kj = rs2[j]; + while (ki < rs1[i + 1] && kj < rs2[j + 1]) { + if (jc1[ki] == jc2[kj]) { + acc += v1[ki] / D[jc1[ki]] * v2[kj]; + ki++; kj++; + } else if (jc1[ki] < jc2[kj]) { + ki++; + } else { + kj++; + } + } + W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + } + } + } + + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvNtrans (off-diag): %f (us)\n", + time * 1e-3 / repeat); + + #pragma omp target exit data map(delete: rs1[0:m+1], jc1[0:nnz], v1[0:nnz], \ + rs2[0:m+1], jc2[0:nnz], v2[0:nnz], \ + D[0:nx], W[0:w_elems]) + + printf("%s\n", errors == 0 ? "PASS" : "FAIL"); + + return 0; +} diff --git a/src/schur-complement-sycl/CMakeLists.txt b/src/schur-complement-sycl/CMakeLists.txt new file mode 100644 index 000000000..3109e22ff --- /dev/null +++ b/src/schur-complement-sycl/CMakeLists.txt @@ -0,0 +1,11 @@ +# schur-complement-sycl/CMakeLists.txt + +set(INC_DIR "${CMAKE_CURRENT_LIST_DIR}/../schur-complement-cuda") + +add_hecbench_benchmark( + NAME schur-complement + MODEL sycl + SOURCES main.cpp + CATEGORIES math + INCLUDE_DIRS ${INC_DIR} +) diff --git a/src/schur-complement-sycl/LICENSE b/src/schur-complement-sycl/LICENSE new file mode 100644 index 000000000..78a030c8d --- /dev/null +++ b/src/schur-complement-sycl/LICENSE @@ -0,0 +1,43 @@ +BSD 3-Clause License. Please also read “Additional BSD Notice” below. + +Copyright (c) 2017, Lawrence Livermore National Laboratory +LLNL-CODE-742473.All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: +i. Redistributions of source code must retain the above copyright notice, this list +of conditions and the disclaimer below. +ii. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the disclaimer (as noted below) in the documentation and/or +other materials provided with the distribution. +iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +endorse or promote products derived from this software without specific prior written +permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Additional BSD Notice +1. This notice is required to be provided under our contract with the U.S. Department +of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +Contract No. DE-AC52-07NA27344 with the DOE. +2. Neither the United States Government nor Lawrence Livermore National Security, LLC +nor any of their employees, makes any warranty, express or implied, or assumes any +liability or responsibility for the accuracy, completeness, or usefulness of any +information, apparatus, product, or process disclosed, or represents that its use would +not infringe privately-owned rights. +3. Also, reference herein to any specific commercial products, process, or services by +trade name, trademark, manufacturer or otherwise does not necessarily constitute or +imply its endorsement, recommendation, or favoring by the United States Government or +Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +herein do not necessarily state or reflect those of the United States Government or +Lawrence Livermore National Security, LLC, and shall not be used for advertising or +product endorsement purposes. diff --git a/src/schur-complement-sycl/Makefile b/src/schur-complement-sycl/Makefile new file mode 100644 index 000000000..f3b5f5e50 --- /dev/null +++ b/src/schur-complement-sycl/Makefile @@ -0,0 +1,82 @@ +#=============================================================================== +# User Options +#=============================================================================== + +# Compiler can be set below, or via environment variable +CC = icpx +OPTIMIZE = yes +DEBUG = no +LAUNCHER = + +GPU = yes +CUDA = no +CUDA_ARCH = sm_70 +HIP = no +HIP_ARCH = gfx908 +#GCC_TOOLCHAIN = "/auto/software/gcc/x86_64/gcc-9.1.0/" + +#=============================================================================== +# Program name & source code list +#=============================================================================== + +program = main + +source = main.cpp + +obj = $(source:.cpp=.o) + +#=============================================================================== +# Sets Flags +#=============================================================================== + +# Standard Flags +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -fsycl -I../schur-complement-cuda \ + --gcc-toolchain=$(GCC_TOOLCHAIN) + +ifeq ($(VENDOR), AdaptiveCpp) + CFLAGS_TMP := $(CFLAGS) + CFLAGS = $(filter-out -fsycl, $(CFLAGS_TMP)) +endif + +# Linker Flags +LDFLAGS = + +ifeq ($(CUDA), yes) + CFLAGS += -fsycl-targets=nvptx64-nvidia-cuda \ + -Xsycl-target-backend --cuda-gpu-arch=$(CUDA_ARCH) +endif + +ifeq ($(HIP), yes) + CFLAGS += -fsycl-targets=amdgcn-amd-amdhsa \ + -Xsycl-target-backend --offload-arch=$(HIP_ARCH) +endif + +# Debug Flags +ifeq ($(DEBUG),yes) + CFLAGS += -g -DDEBUG + LDFLAGS += -g +endif + +# Optimization Flags +ifeq ($(OPTIMIZE),yes) + CFLAGS += -O3 +endif + +ifeq ($(GPU),yes) + CFLAGS +=-DUSE_GPU +endif +#=============================================================================== +# Targets to Build +#=============================================================================== + +$(program): $(obj) + $(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS) + +%.o: %.cpp ../schur-complement-cuda/reference.h Makefile + $(CC) $(CFLAGS) -c $< -o $@ + +clean: + rm -rf $(program) $(obj) + +run: $(program) + $(LAUNCHER) ./$(program) 4096 32 100 diff --git a/src/schur-complement-sycl/main.cpp b/src/schur-complement-sycl/main.cpp new file mode 100644 index 000000000..6103a2f3e --- /dev/null +++ b/src/schur-complement-sycl/main.cpp @@ -0,0 +1,180 @@ +#include +#include +#include +#include +#include +#include +#include "reference.h" + +int main(int argc, char* argv[]) +{ + if (argc != 4) { + printf("Usage: %s \n", argv[0]); + return 1; + } + + const int m = atoi(argv[1]); + const int nnz_row = atoi(argv[2]); + const int repeat = atoi(argv[3]); + + const int nx = 8 * nnz_row + 1024; + const double alpha = -1.0; + + std::vector h_rs1, h_jc1, h_rs2, h_jc2; + std::vector h_v1, h_v2, h_D; + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + + const int nnz = (int)h_v1.size(); + const size_t w_elems = (size_t)m * m; + const size_t w_bytes = w_elems * sizeof(double); + +#ifdef USE_GPU + sycl::queue q(sycl::gpu_selector_v, sycl::property::queue::in_order()); +#else + sycl::queue q(sycl::cpu_selector_v, sycl::property::queue::in_order()); +#endif + + int *d_rs1 = sycl::malloc_device(m + 1, q); + int *d_jc1 = sycl::malloc_device(nnz, q); + double *d_v1 = sycl::malloc_device(nnz, q); + int *d_rs2 = sycl::malloc_device(m + 1, q); + int *d_jc2 = sycl::malloc_device(nnz, q); + double *d_v2 = sycl::malloc_device(nnz, q); + double *d_D = sycl::malloc_device(nx, q); + double *d_W = sycl::malloc_device(w_elems, q); + + q.memcpy(d_rs1, h_rs1.data(), (m + 1) * sizeof(int)); + q.memcpy(d_jc1, h_jc1.data(), nnz * sizeof(int)); + q.memcpy(d_v1, h_v1.data(), nnz * sizeof(double)); + q.memcpy(d_rs2, h_rs2.data(), (m + 1) * sizeof(int)); + q.memcpy(d_jc2, h_jc2.data(), nnz * sizeof(int)); + q.memcpy(d_v2, h_v2.data(), nnz * sizeof(double)); + q.memcpy(d_D, h_D.data(), nx * sizeof(double)); + q.wait(); + + const int m_W = m; + const int row_dest_start = 0, col_dest_start = 0; + + // Each work-item owns a single (i,j) row pair (O(m^2) parallelism) rather + // than a whole row (O(m)). The 2D nd_range puts j on the fast dimension so + // writes to W[i][j] are contiguous within a work-group; the x extent (32) is + // wavefront-aligned. + const int BX = 32, BY = 8; + const size_t gx = ((size_t)(m + BX - 1) / BX) * BX; + const size_t gy = ((size_t)(m + BY - 1) / BY) * BY; + const sycl::range<2> local(BY, BX); + const sycl::range<2> glob(gy, gx); + + std::vector h_W(w_elems); + std::vector h_ref(w_elems); + int errors = 0; + + // W += alpha * J * D^{-1} * J^T : one work-item owns row i, merges it against + // every later row j and writes the upper triangle of the diagonal block. + auto run_mdinvmtrans_diag = [&]() { + q.parallel_for(sycl::nd_range<2>(glob, local), [=](sycl::nd_item<2> item) { + const int i = item.get_global_id(0); + const int j = item.get_global_id(1); + if (i >= m || j >= m || j < i) return; + double acc = 0.0; + if (i == j) { + for (int k = d_rs1[i]; k < d_rs1[i + 1]; k++) + acc += d_v1[k] / d_D[d_jc1[k]] * d_v1[k]; + } else { + int ki = d_rs1[i], kj = d_rs1[j]; + const int kie = d_rs1[i + 1], kje = d_rs1[j + 1]; + while (ki < kie && kj < kje) { + const int ci = d_jc1[ki], cj = d_jc1[kj]; + if (ci == cj) { acc += d_v1[ki] / d_D[ci] * d_v1[kj]; ki++; kj++; } + else if (ci < cj) ki++; + else kj++; + } + } + d_W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + }); + }; + + // W += alpha * J1 * D^{-1} * J2^T : one work-item owns row i of J1, merges it + // against every row j of J2 and writes the full off-diagonal block. + auto run_mdinvntrans = [&]() { + q.parallel_for(sycl::nd_range<2>(glob, local), [=](sycl::nd_item<2> item) { + const int i = item.get_global_id(0); + const int j = item.get_global_id(1); + if (i >= m || j >= m) return; + double acc = 0.0; + int ki = d_rs1[i], kj = d_rs2[j]; + const int kie = d_rs1[i + 1], kje = d_rs2[j + 1]; + while (ki < kie && kj < kje) { + const int ci = d_jc1[ki], cj = d_jc2[kj]; + if (ci == cj) { acc += d_v1[ki] / d_D[ci] * d_v2[kj]; ki++; kj++; } + else if (ci < cj) ki++; + else kj++; + } + d_W[(size_t)(i + row_dest_start) * m_W + j + col_dest_start] += alpha * acc; + }); + }; + + // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- + + // host/device correctness check (run once, verify against reference) before timing + q.memset(d_W, 0, w_bytes).wait(); + run_mdinvmtrans_diag(); + q.memcpy(h_W.data(), d_W, w_bytes).wait(); + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvmtrans_diag(m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (int i = 0; i < m; i++) + for (int j = i; j < m; j++) + if (!close_enough(h_W[(size_t)i * m + j], h_ref[(size_t)i * m + j], 1e-10)) { + errors++; i = m; break; + } + + // benchmark + q.memset(d_W, 0, w_bytes).wait(); + auto start = std::chrono::steady_clock::now(); + + for (int i = 0; i < repeat; i++) run_mdinvmtrans_diag(); + q.wait(); + + auto end = std::chrono::steady_clock::now(); + auto time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvMtrans (diag block): %f (us)\n", + time * 1e-3 / repeat); + + // --- off-diagonal block: W += alpha * J1 D^{-1} J2^T --------------------- + + // host/device correctness check before timing + q.memset(d_W, 0, w_bytes).wait(); + run_mdinvntrans(); + q.memcpy(h_W.data(), d_W, w_bytes).wait(); + + std::fill(h_ref.begin(), h_ref.end(), 0.0); + reference_mdinvntrans(m, m, m, h_rs1.data(), h_jc1.data(), h_v1.data(), + h_rs2.data(), h_jc2.data(), h_v2.data(), h_D.data(), + 0, 0, alpha, h_ref.data()); + for (size_t k = 0; k < w_elems; k++) + if (!close_enough(h_W[k], h_ref[k], 1e-10)) { errors++; break; } + + // benchmark + q.memset(d_W, 0, w_bytes).wait(); + start = std::chrono::steady_clock::now(); + + for (int i = 0; i < repeat; i++) run_mdinvntrans(); + q.wait(); + + end = std::chrono::steady_clock::now(); + time = std::chrono::duration_cast(end - start).count(); + printf("Average execution time of addMDinvNtrans (off-diag): %f (us)\n", + time * 1e-3 / repeat); + + printf("%s\n", errors == 0 ? "PASS" : "FAIL"); + + sycl::free(d_rs1, q); sycl::free(d_jc1, q); sycl::free(d_v1, q); + sycl::free(d_rs2, q); sycl::free(d_jc2, q); sycl::free(d_v2, q); + sycl::free(d_D, q); sycl::free(d_W, q); + + return 0; +} From b09b6ba87dd942274768926834dc0b5c169ebf9a Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Sat, 15 Aug 2026 14:48:50 -0700 Subject: [PATCH 2/4] [schur-complement] validate the command-line arguments and align the SYCL and OpenMP kernel launches The CSR row pointers and column indices are 32-bit, so reject non-positive sizes and sizes whose nonzero count would overflow an int instead of silently generating a corrupt matrix. Also share the block size macros in the SYCL version and set the team size explicitly in the OpenMP version. Co-authored-by: Cursor --- src/schur-complement-cuda/main.cu | 10 ++++++++++ src/schur-complement-hip/main.cu | 10 ++++++++++ src/schur-complement-omp/main.cpp | 26 ++++++++++++++++++-------- src/schur-complement-sycl/main.cpp | 28 +++++++++++++++++++++------- 4 files changed, 59 insertions(+), 15 deletions(-) diff --git a/src/schur-complement-cuda/main.cu b/src/schur-complement-cuda/main.cu index 953f165c8..8aa7c6d4c 100644 --- a/src/schur-complement-cuda/main.cu +++ b/src/schur-complement-cuda/main.cu @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -93,6 +94,15 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); + // the CSR row pointers and column indices are 32-bit, as in HiOp + if (m <= 0 || nnz_row <= 0 || repeat <= 0 || + (long long)m * nnz_row > INT_MAX || + (long long)8 * nnz_row + 1024 > INT_MAX) { + printf("Invalid arguments: , and must be " + "positive, and the number of nonzeros must fit in a 32-bit int\n"); + return 1; + } + // number of variables (columns of J); a few times wider than nnz_row so the // per-row column lists overlap only partially, exercising the merge loop const int nx = 8 * nnz_row + 1024; diff --git a/src/schur-complement-hip/main.cu b/src/schur-complement-hip/main.cu index db0b1e0b8..84697b4c1 100644 --- a/src/schur-complement-hip/main.cu +++ b/src/schur-complement-hip/main.cu @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -93,6 +94,15 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); + // the CSR row pointers and column indices are 32-bit, as in HiOp + if (m <= 0 || nnz_row <= 0 || repeat <= 0 || + (long long)m * nnz_row > INT_MAX || + (long long)8 * nnz_row + 1024 > INT_MAX) { + printf("Invalid arguments: , and must be " + "positive, and the number of nonzeros must fit in a 32-bit int\n"); + return 1; + } + const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; diff --git a/src/schur-complement-omp/main.cpp b/src/schur-complement-omp/main.cpp index 63806e901..916b1ff20 100644 --- a/src/schur-complement-omp/main.cpp +++ b/src/schur-complement-omp/main.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -17,6 +18,15 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); + // the CSR row pointers and column indices are 32-bit, as in HiOp + if (m <= 0 || nnz_row <= 0 || repeat <= 0 || + (long long)m * nnz_row > INT_MAX || + (long long)8 * nnz_row + 1024 > INT_MAX) { + printf("Invalid arguments: , and must be " + "positive, and the number of nonzeros must fit in a 32-bit int\n"); + return 1; + } + const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; @@ -47,10 +57,10 @@ int main(int argc, char* argv[]) // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- // host/device correctness check (run once, verify against reference) before timing - #pragma omp target teams distribute parallel for thread_limit(128) + #pragma omp target teams distribute parallel for for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; - #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + #pragma omp target teams distribute parallel for collapse(2) num_threads(256) for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { if (j < i) continue; @@ -80,14 +90,14 @@ int main(int argc, char* argv[]) } // benchmark - #pragma omp target teams distribute parallel for thread_limit(128) + #pragma omp target teams distribute parallel for for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; auto start = std::chrono::steady_clock::now(); for (int r = 0; r < repeat; r++) { // one iteration per (i,j) row pair (collapsed) instead of per row i - #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + #pragma omp target teams distribute parallel for collapse(2) num_threads(256) for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { if (j < i) continue; @@ -115,10 +125,10 @@ int main(int argc, char* argv[]) // --- off-diagonal block: W += alpha * J1 D^{-1} J2^T --------------------- // host/device correctness check before timing - #pragma omp target teams distribute parallel for thread_limit(128) + #pragma omp target teams distribute parallel for for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; - #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + #pragma omp target teams distribute parallel for collapse(2) num_threads(256) for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { double acc = 0.0; @@ -140,13 +150,13 @@ int main(int argc, char* argv[]) if (!close_enough(h_W[k], h_ref[k], 1e-10)) { errors++; break; } // benchmark - #pragma omp target teams distribute parallel for thread_limit(128) + #pragma omp target teams distribute parallel for for (size_t k = 0; k < w_elems; k++) W[k] = 0.0; start = std::chrono::steady_clock::now(); for (int r = 0; r < repeat; r++) { - #pragma omp target teams distribute parallel for collapse(2) thread_limit(128) + #pragma omp target teams distribute parallel for collapse(2) num_threads(256) for (int i = 0; i < m; i++) { for (int j = 0; j < m; j++) { double acc = 0.0; diff --git a/src/schur-complement-sycl/main.cpp b/src/schur-complement-sycl/main.cpp index 6103a2f3e..ca3cd93ba 100644 --- a/src/schur-complement-sycl/main.cpp +++ b/src/schur-complement-sycl/main.cpp @@ -1,11 +1,20 @@ #include #include +#include #include #include #include #include #include "reference.h" +// Original HiOp mapping is one thread per row i, looping over all partner rows +// j -- only O(m) work items, which leaves the GPU underused. Here each thread +// instead owns a single (i,j) row pair and performs one merge, raising the +// parallelism to O(m^2). Threads are laid out 2D with the fast (x) dimension on +// j so the writes to W[i][j] are coalesced; BLOCK_X is warp-aligned. +#define BLOCK_X 32 +#define BLOCK_Y 8 + int main(int argc, char* argv[]) { if (argc != 4) { @@ -17,6 +26,15 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); + // the CSR row pointers and column indices are 32-bit, as in HiOp + if (m <= 0 || nnz_row <= 0 || repeat <= 0 || + (long long)m * nnz_row > INT_MAX || + (long long)8 * nnz_row + 1024 > INT_MAX) { + printf("Invalid arguments: , and must be " + "positive, and the number of nonzeros must fit in a 32-bit int\n"); + return 1; + } + const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; @@ -57,13 +75,9 @@ int main(int argc, char* argv[]) const int m_W = m; const int row_dest_start = 0, col_dest_start = 0; - // Each work-item owns a single (i,j) row pair (O(m^2) parallelism) rather - // than a whole row (O(m)). The 2D nd_range puts j on the fast dimension so - // writes to W[i][j] are contiguous within a work-group; the x extent (32) is - // wavefront-aligned. - const int BX = 32, BY = 8; - const size_t gx = ((size_t)(m + BX - 1) / BX) * BX; - const size_t gy = ((size_t)(m + BY - 1) / BY) * BY; + const int BX = BLOCK_X, BY = BLOCK_Y; + const size_t gx = (m + BX - 1) / BX * BX; + const size_t gy = (m + BY - 1) / BY * BY; const sycl::range<2> local(BY, BX); const sycl::range<2> glob(gy, gx); From fa1da9235eb5c59266bdc43a3adecc506880a0ca Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Fri, 21 Aug 2026 13:36:37 -0700 Subject: [PATCH 3/4] [schur-complement] share the argument validation and report allocation and launch failures - Move the command-line checks into reference.h so the four versions cannot drift apart, and bound the number of rows as well: m + 1 row pointers and rounding m up to the block size have to stay within an int. - Report a failed host allocation instead of terminating on bad_alloc, and a failed device allocation in the SYCL version instead of dereferencing null in the kernel. - Check the CUDA and HIP kernel launches, so an invalid launch configuration is reported rather than showing up as a verification failure. - Fix the two SYCL kernel comments that still described the original HiOp one-work-item-per-row mapping. Co-authored-by: Cursor --- src/schur-complement-cuda/main.cu | 32 ++++++++++------- src/schur-complement-cuda/reference.h | 13 +++++++ src/schur-complement-hip/main.cu | 32 ++++++++++------- src/schur-complement-omp/main.cpp | 28 ++++++++------- src/schur-complement-sycl/main.cpp | 50 +++++++++++++++++---------- 5 files changed, 98 insertions(+), 57 deletions(-) diff --git a/src/schur-complement-cuda/main.cu b/src/schur-complement-cuda/main.cu index 8aa7c6d4c..a61ada1e1 100644 --- a/src/schur-complement-cuda/main.cu +++ b/src/schur-complement-cuda/main.cu @@ -1,8 +1,8 @@ #include #include -#include #include #include +#include #include #include #include "reference.h" @@ -94,10 +94,7 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); - // the CSR row pointers and column indices are 32-bit, as in HiOp - if (m <= 0 || nnz_row <= 0 || repeat <= 0 || - (long long)m * nnz_row > INT_MAX || - (long long)8 * nnz_row + 1024 > INT_MAX) { + if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " "positive, and the number of nonzeros must fit in a 32-bit int\n"); return 1; @@ -108,15 +105,24 @@ int main(int argc, char* argv[]) const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; // HiOp assembles the Schur complement with -1 + const size_t w_elems = (size_t)m * m; + const size_t w_bytes = w_elems * sizeof(double); + std::vector h_rs1, h_jc1, h_rs2, h_jc2; - std::vector h_v1, h_v2, h_D; - gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); - gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); - gen_diag(nx, 789, h_D); + std::vector h_v1, h_v2, h_D, h_W, h_ref; + try { + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + h_W.resize(w_elems); + h_ref.resize(w_elems); + } catch (const std::bad_alloc&) { + printf("Failed to allocate the host buffers: the dense block alone needs " + "%zu bytes\n", w_bytes); + return 1; + } const int nnz = (int)h_v1.size(); - const size_t w_elems = (size_t)m * m; - const size_t w_bytes = w_elems * sizeof(double); int *d_rs1, *d_jc1, *d_rs2, *d_jc2; double *d_v1, *d_v2, *d_D, *d_W; @@ -140,8 +146,6 @@ int main(int argc, char* argv[]) const dim3 block(BLOCK_X, BLOCK_Y); const dim3 grid((m + BLOCK_X - 1) / BLOCK_X, (m + BLOCK_Y - 1) / BLOCK_Y); - std::vector h_W(w_elems); - std::vector h_ref(w_elems); int errors = 0; // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- @@ -149,6 +153,7 @@ int main(int argc, char* argv[]) // host/device correctness check (run once, verify against reference) before timing CHECK(cudaMemset(d_W, 0, w_bytes)); mdinvmtrans_diag<<>>(m, d_rs1, d_jc1, d_v1, d_D, 0, 0, alpha, d_W, m); + CHECK(cudaGetLastError()); CHECK(cudaMemcpy(h_W.data(), d_W, w_bytes, cudaMemcpyDeviceToHost)); std::fill(h_ref.begin(), h_ref.end(), 0.0); @@ -180,6 +185,7 @@ int main(int argc, char* argv[]) // host/device correctness check before timing CHECK(cudaMemset(d_W, 0, w_bytes)); mdinvntrans<<>>(m, m, d_rs1, d_jc1, d_v1, d_rs2, d_jc2, d_v2, d_D, 0, 0, alpha, d_W, m); + CHECK(cudaGetLastError()); CHECK(cudaMemcpy(h_W.data(), d_W, w_bytes, cudaMemcpyDeviceToHost)); std::fill(h_ref.begin(), h_ref.end(), 0.0); diff --git a/src/schur-complement-cuda/reference.h b/src/schur-complement-cuda/reference.h index b55fde3fb..805f56b49 100644 --- a/src/schur-complement-cuda/reference.h +++ b/src/schur-complement-cuda/reference.h @@ -2,6 +2,7 @@ #define REFERENCE_H #include +#include #include #include #include @@ -15,6 +16,18 @@ // D is a diagonal (stored as a vector), and W is a dense matrix into whose // upper triangle the (symmetric) result is scattered. +// The CSR row pointers and column indices are 32-bit, as in HiOp, so reject +// sizes whose derived quantities would overflow rather than silently +// generating a corrupt matrix. The 1024 headroom on m covers both the m + 1 +// row pointers and rounding m up to a multiple of the block size. +static bool valid_problem_size(int m, int nnz_row, int repeat) +{ + return m > 0 && nnz_row > 0 && repeat > 0 && + m <= INT_MAX - 1024 && // m + 1, and the launch range + nnz_row <= (INT_MAX - 1024) / 8 && // nx = 8 * nnz_row + 1024 + (long long)m * nnz_row <= INT_MAX; // number of nonzeros +} + // Build a CSR matrix with `m` rows, `nx` columns and exactly `nnz_row` // nonzeros per row (column indices sorted ascending within each row). static void gen_csr(int m, int nx, int nnz_row, unsigned seed, diff --git a/src/schur-complement-hip/main.cu b/src/schur-complement-hip/main.cu index 84697b4c1..5cd54b0c1 100644 --- a/src/schur-complement-hip/main.cu +++ b/src/schur-complement-hip/main.cu @@ -1,8 +1,8 @@ #include #include -#include #include #include +#include #include #include #include "reference.h" @@ -94,10 +94,7 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); - // the CSR row pointers and column indices are 32-bit, as in HiOp - if (m <= 0 || nnz_row <= 0 || repeat <= 0 || - (long long)m * nnz_row > INT_MAX || - (long long)8 * nnz_row + 1024 > INT_MAX) { + if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " "positive, and the number of nonzeros must fit in a 32-bit int\n"); return 1; @@ -106,15 +103,24 @@ int main(int argc, char* argv[]) const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; + const size_t w_elems = (size_t)m * m; + const size_t w_bytes = w_elems * sizeof(double); + std::vector h_rs1, h_jc1, h_rs2, h_jc2; - std::vector h_v1, h_v2, h_D; - gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); - gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); - gen_diag(nx, 789, h_D); + std::vector h_v1, h_v2, h_D, h_W, h_ref; + try { + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + h_W.resize(w_elems); + h_ref.resize(w_elems); + } catch (const std::bad_alloc&) { + printf("Failed to allocate the host buffers: the dense block alone needs " + "%zu bytes\n", w_bytes); + return 1; + } const int nnz = (int)h_v1.size(); - const size_t w_elems = (size_t)m * m; - const size_t w_bytes = w_elems * sizeof(double); int *d_rs1, *d_jc1, *d_rs2, *d_jc2; double *d_v1, *d_v2, *d_D, *d_W; @@ -138,8 +144,6 @@ int main(int argc, char* argv[]) const dim3 block(BLOCK_X, BLOCK_Y); const dim3 grid((m + BLOCK_X - 1) / BLOCK_X, (m + BLOCK_Y - 1) / BLOCK_Y); - std::vector h_W(w_elems); - std::vector h_ref(w_elems); int errors = 0; // --- diagonal block: W += alpha * J D^{-1} J^T --------------------------- @@ -147,6 +151,7 @@ int main(int argc, char* argv[]) // host/device correctness check (run once, verify against reference) before timing CHECK(hipMemset(d_W, 0, w_bytes)); mdinvmtrans_diag<<>>(m, d_rs1, d_jc1, d_v1, d_D, 0, 0, alpha, d_W, m); + CHECK(hipGetLastError()); CHECK(hipMemcpy(h_W.data(), d_W, w_bytes, hipMemcpyDeviceToHost)); std::fill(h_ref.begin(), h_ref.end(), 0.0); @@ -178,6 +183,7 @@ int main(int argc, char* argv[]) // host/device correctness check before timing CHECK(hipMemset(d_W, 0, w_bytes)); mdinvntrans<<>>(m, m, d_rs1, d_jc1, d_v1, d_rs2, d_jc2, d_v2, d_D, 0, 0, alpha, d_W, m); + CHECK(hipGetLastError()); CHECK(hipMemcpy(h_W.data(), d_W, w_bytes, hipMemcpyDeviceToHost)); std::fill(h_ref.begin(), h_ref.end(), 0.0); diff --git a/src/schur-complement-omp/main.cpp b/src/schur-complement-omp/main.cpp index 916b1ff20..b80b3c27f 100644 --- a/src/schur-complement-omp/main.cpp +++ b/src/schur-complement-omp/main.cpp @@ -1,8 +1,8 @@ #include #include -#include #include #include +#include #include #include #include "reference.h" @@ -18,10 +18,7 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); - // the CSR row pointers and column indices are 32-bit, as in HiOp - if (m <= 0 || nnz_row <= 0 || repeat <= 0 || - (long long)m * nnz_row > INT_MAX || - (long long)8 * nnz_row + 1024 > INT_MAX) { + if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " "positive, and the number of nonzeros must fit in a 32-bit int\n"); return 1; @@ -30,14 +27,23 @@ int main(int argc, char* argv[]) const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; + const size_t w_elems = (size_t)m * m; + std::vector h_rs1, h_jc1, h_rs2, h_jc2; - std::vector h_v1, h_v2, h_D; - gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); - gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); - gen_diag(nx, 789, h_D); + std::vector h_v1, h_v2, h_D, h_W, h_ref; + try { + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + h_W.resize(w_elems); + h_ref.resize(w_elems); + } catch (const std::bad_alloc&) { + printf("Failed to allocate the host buffers: the dense block alone needs " + "%zu bytes\n", w_elems * sizeof(double)); + return 1; + } const int nnz = (int)h_v1.size(); - const size_t w_elems = (size_t)m * m; const int m_W = m; const int row_dest_start = 0, col_dest_start = 0; @@ -45,9 +51,7 @@ int main(int argc, char* argv[]) int* rs1 = h_rs1.data(); int* jc1 = h_jc1.data(); double* v1 = h_v1.data(); int* rs2 = h_rs2.data(); int* jc2 = h_jc2.data(); double* v2 = h_v2.data(); double* D = h_D.data(); - std::vector h_W(w_elems, 0.0); double* W = h_W.data(); - std::vector h_ref(w_elems); int errors = 0; #pragma omp target enter data map(to: rs1[0:m+1], jc1[0:nnz], v1[0:nnz], \ diff --git a/src/schur-complement-sycl/main.cpp b/src/schur-complement-sycl/main.cpp index ca3cd93ba..74845308e 100644 --- a/src/schur-complement-sycl/main.cpp +++ b/src/schur-complement-sycl/main.cpp @@ -1,8 +1,8 @@ #include #include -#include #include #include +#include #include #include #include "reference.h" @@ -26,10 +26,7 @@ int main(int argc, char* argv[]) const int nnz_row = atoi(argv[2]); const int repeat = atoi(argv[3]); - // the CSR row pointers and column indices are 32-bit, as in HiOp - if (m <= 0 || nnz_row <= 0 || repeat <= 0 || - (long long)m * nnz_row > INT_MAX || - (long long)8 * nnz_row + 1024 > INT_MAX) { + if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " "positive, and the number of nonzeros must fit in a 32-bit int\n"); return 1; @@ -38,15 +35,24 @@ int main(int argc, char* argv[]) const int nx = 8 * nnz_row + 1024; const double alpha = -1.0; + const size_t w_elems = (size_t)m * m; + const size_t w_bytes = w_elems * sizeof(double); + std::vector h_rs1, h_jc1, h_rs2, h_jc2; - std::vector h_v1, h_v2, h_D; - gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); - gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); - gen_diag(nx, 789, h_D); + std::vector h_v1, h_v2, h_D, h_W, h_ref; + try { + gen_csr(m, nx, nnz_row, 123, h_rs1, h_jc1, h_v1); + gen_csr(m, nx, nnz_row, 456, h_rs2, h_jc2, h_v2); + gen_diag(nx, 789, h_D); + h_W.resize(w_elems); + h_ref.resize(w_elems); + } catch (const std::bad_alloc&) { + printf("Failed to allocate the host buffers: the dense block alone needs " + "%zu bytes\n", w_bytes); + return 1; + } const int nnz = (int)h_v1.size(); - const size_t w_elems = (size_t)m * m; - const size_t w_bytes = w_elems * sizeof(double); #ifdef USE_GPU sycl::queue q(sycl::gpu_selector_v, sycl::property::queue::in_order()); @@ -63,6 +69,14 @@ int main(int argc, char* argv[]) double *d_D = sycl::malloc_device(nx, q); double *d_W = sycl::malloc_device(w_elems, q); + if (d_rs1 == nullptr || d_jc1 == nullptr || d_v1 == nullptr || + d_rs2 == nullptr || d_jc2 == nullptr || d_v2 == nullptr || + d_D == nullptr || d_W == nullptr) { + printf("Failed to allocate the device buffers: the dense block alone needs " + "%zu bytes\n", w_bytes); + return 1; + } + q.memcpy(d_rs1, h_rs1.data(), (m + 1) * sizeof(int)); q.memcpy(d_jc1, h_jc1.data(), nnz * sizeof(int)); q.memcpy(d_v1, h_v1.data(), nnz * sizeof(double)); @@ -76,17 +90,15 @@ int main(int argc, char* argv[]) const int row_dest_start = 0, col_dest_start = 0; const int BX = BLOCK_X, BY = BLOCK_Y; - const size_t gx = (m + BX - 1) / BX * BX; - const size_t gy = (m + BY - 1) / BY * BY; + const size_t gx = ((size_t)m + BX - 1) / BX * BX; + const size_t gy = ((size_t)m + BY - 1) / BY * BY; const sycl::range<2> local(BY, BX); const sycl::range<2> glob(gy, gx); - std::vector h_W(w_elems); - std::vector h_ref(w_elems); int errors = 0; - // W += alpha * J * D^{-1} * J^T : one work-item owns row i, merges it against - // every later row j and writes the upper triangle of the diagonal block. + // W += alpha * J * D^{-1} * J^T : work-item (i,j) merges row i against row j + // and writes one entry of the upper triangle of the diagonal block. auto run_mdinvmtrans_diag = [&]() { q.parallel_for(sycl::nd_range<2>(glob, local), [=](sycl::nd_item<2> item) { const int i = item.get_global_id(0); @@ -110,8 +122,8 @@ int main(int argc, char* argv[]) }); }; - // W += alpha * J1 * D^{-1} * J2^T : one work-item owns row i of J1, merges it - // against every row j of J2 and writes the full off-diagonal block. + // W += alpha * J1 * D^{-1} * J2^T : work-item (i,j) merges row i of J1 + // against row j of J2 and writes one entry of the off-diagonal block. auto run_mdinvntrans = [&]() { q.parallel_for(sycl::nd_range<2>(glob, local), [=](sycl::nd_item<2> item) { const int i = item.get_global_id(0); From fb78012cf5d104ce32344ad339e93a3bc78230c2 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Tue, 8 Sep 2026 06:46:33 -0700 Subject: [PATCH 4/4] [schur-complement] validate dense block size Co-authored-by: Cursor --- src/schur-complement-cuda/main.cu | 2 +- src/schur-complement-cuda/reference.h | 7 +++++-- src/schur-complement-hip/main.cu | 2 +- src/schur-complement-omp/main.cpp | 2 +- src/schur-complement-sycl/main.cpp | 2 +- 5 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/schur-complement-cuda/main.cu b/src/schur-complement-cuda/main.cu index a61ada1e1..9258936be 100644 --- a/src/schur-complement-cuda/main.cu +++ b/src/schur-complement-cuda/main.cu @@ -96,7 +96,7 @@ int main(int argc, char* argv[]) if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " - "positive, and the number of nonzeros must fit in a 32-bit int\n"); + "positive, and all derived sizes must fit their integer types\n"); return 1; } diff --git a/src/schur-complement-cuda/reference.h b/src/schur-complement-cuda/reference.h index 805f56b49..293a18ec2 100644 --- a/src/schur-complement-cuda/reference.h +++ b/src/schur-complement-cuda/reference.h @@ -2,6 +2,7 @@ #define REFERENCE_H #include +#include #include #include #include @@ -19,13 +20,15 @@ // The CSR row pointers and column indices are 32-bit, as in HiOp, so reject // sizes whose derived quantities would overflow rather than silently // generating a corrupt matrix. The 1024 headroom on m covers both the m + 1 -// row pointers and rounding m up to a multiple of the block size. +// row pointers and rounding m up to a multiple of the block size. Also ensure +// that the byte size of the dense m-by-m double block fits in size_t. static bool valid_problem_size(int m, int nnz_row, int repeat) { return m > 0 && nnz_row > 0 && repeat > 0 && m <= INT_MAX - 1024 && // m + 1, and the launch range nnz_row <= (INT_MAX - 1024) / 8 && // nx = 8 * nnz_row + 1024 - (long long)m * nnz_row <= INT_MAX; // number of nonzeros + (long long)m * nnz_row <= INT_MAX && // number of nonzeros + (size_t)m <= SIZE_MAX / sizeof(double) / (size_t)m; // dense bytes } // Build a CSR matrix with `m` rows, `nx` columns and exactly `nnz_row` diff --git a/src/schur-complement-hip/main.cu b/src/schur-complement-hip/main.cu index 5cd54b0c1..6b1b872d8 100644 --- a/src/schur-complement-hip/main.cu +++ b/src/schur-complement-hip/main.cu @@ -96,7 +96,7 @@ int main(int argc, char* argv[]) if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " - "positive, and the number of nonzeros must fit in a 32-bit int\n"); + "positive, and all derived sizes must fit their integer types\n"); return 1; } diff --git a/src/schur-complement-omp/main.cpp b/src/schur-complement-omp/main.cpp index b80b3c27f..af9b2efbe 100644 --- a/src/schur-complement-omp/main.cpp +++ b/src/schur-complement-omp/main.cpp @@ -20,7 +20,7 @@ int main(int argc, char* argv[]) if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " - "positive, and the number of nonzeros must fit in a 32-bit int\n"); + "positive, and all derived sizes must fit their integer types\n"); return 1; } diff --git a/src/schur-complement-sycl/main.cpp b/src/schur-complement-sycl/main.cpp index 74845308e..8439aab24 100644 --- a/src/schur-complement-sycl/main.cpp +++ b/src/schur-complement-sycl/main.cpp @@ -28,7 +28,7 @@ int main(int argc, char* argv[]) if (!valid_problem_size(m, nnz_row, repeat)) { printf("Invalid arguments: , and must be " - "positive, and the number of nonzeros must fit in a 32-bit int\n"); + "positive, and all derived sizes must fit their integer types\n"); return 1; }