diff --git a/CMakeLists.txt b/CMakeLists.txt index 65dfdba..c22f731 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,6 +28,13 @@ cmake_minimum_required(VERSION 3.25 FATAL_ERROR) project(linalglib LANGUAGES CXX) +if(MSVC) + # Keep the project and all fetched dependencies on the DLL-based MSVC runtime + # so that targets like yaml-cpp do not end up linked against /MTd while the + # executable is built against /MDd. + set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$:Debug>DLL") +endif() + set(CMAKE_CXX_STANDARD 17) option(ENABLE_COVERAGE "Build with code coverage support" OFF) diff --git a/clients/testing/test_arguments.h b/clients/testing/test_arguments.h index cbbd599..a2459ec 100644 --- a/clients/testing/test_arguments.h +++ b/clients/testing/test_arguments.h @@ -61,10 +61,7 @@ namespace testing std::string generate_test_name() const { std::string name = group; - if(this->backend != backend::CPU) - { - name += "_" + backend_to_string(this->backend); - } + name += "_" + backend_to_string(this->backend); if(this->uplo != uplo::lower) { name += "_" + uplo_to_string(this->uplo); diff --git a/clients/testing/tests/test_symmetric_ruiz_scaling.yaml b/clients/testing/tests/test_symmetric_ruiz_scaling.yaml index 877723a..1c32e3e 100644 --- a/clients/testing/tests/test_symmetric_ruiz_scaling.yaml +++ b/clients/testing/tests/test_symmetric_ruiz_scaling.yaml @@ -9,7 +9,7 @@ Tests: "matrices/SPD/nos7/nos7.mtx"] max_iters: [30] tol: [1e-6] - backend: [GPU] + backend: [CPU] small: matrix_file: ["matrices/SPD/bcsstm02/bcsstm02.mtx", diff --git a/library/include/csr_matrix.h b/library/include/csr_matrix.h index f732c80..5a19402 100644 --- a/library/include/csr_matrix.h +++ b/library/include/csr_matrix.h @@ -32,6 +32,14 @@ #include "vector.h" +#include "linalg_types.h" + +// struct csrmv_descr; +// struct csrgemm_descr; +// struct csrsv_descr; +// struct csric0_descr; +// struct csrilu0_descr; + /*! \file * \brief csr_matrx.h provides class for CSR sparse matrices */ @@ -95,6 +103,12 @@ namespace linalg /*! \brief Flag indicating if the matrix data is currently on the host (CPU) or device (GPU). */ bool on_host; + csrmv_descr* descr_mv; + csrgemm_descr* descr_gemm; + csrtrsv_descr* descr_sv; + csric0_descr* descr_ic; + csrilu0_descr* descr_ilu; + public: /*! \brief Default constructor. * Initializes an empty CSR matrix with zero dimensions and no non-zero elements. diff --git a/library/include/linalg_enums.h b/library/include/linalg_enums.h index 1904e6c..9a33e1e 100644 --- a/library/include/linalg_enums.h +++ b/library/include/linalg_enums.h @@ -65,7 +65,8 @@ namespace linalg default_algorithm, /*!< Default algorithm */ merge_path, rowsplit, - nnzsplit + nnzsplit, + lrb }; /*! \brief Enumeration for CSR matrix-matrix addition algorithms. diff --git a/library/src/backend/device/cuda/csrmv_kernels.cuh b/library/src/backend/device/cuda/csrmv_kernels.cuh index 1d86a75..478c6c8 100644 --- a/library/src/backend/device/cuda/csrmv_kernels.cuh +++ b/library/src/backend/device/cuda/csrmv_kernels.cuh @@ -30,16 +30,16 @@ #include "common.cuh" template -__global__ void csrmv_vector_kernel(int m, - int n, - int nnz, - const T alpha, - const int* __restrict__ csr_row_ptr, - const int* __restrict__ csr_col_ind, - const T* __restrict__ csr_val, - const T* __restrict__ x, - const T beta, - T* __restrict__ y) +__global__ void csrmv_row_split_kernel(int m, + int n, + int nnz, + const T alpha, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + const T* __restrict__ csr_val, + const T* __restrict__ x, + const T beta, + T* __restrict__ y) { const int tid = threadIdx.x; const int bid = blockIdx.x; @@ -103,16 +103,16 @@ __device__ inline int } template -__global__ void csrmv_stream_kernel(int m, - int n, - int nnz, - const T alpha, - const int* __restrict__ csr_row_ptr, - const int* __restrict__ csr_col_ind, - const T* __restrict__ csr_val, - const T* __restrict__ x, - const T beta, - T* __restrict__ y) +__global__ void csrmv_nnz_split_kernel(int m, + int n, + int nnz, + const T alpha, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + const T* __restrict__ csr_val, + const T* __restrict__ x, + const T beta, + T* __restrict__ y) { const int tid = threadIdx.x; const int bid = blockIdx.x; @@ -232,4 +232,336 @@ __global__ void csrmv_stream_kernel(int m, } } + +__device__ int ilog2(unsigned int x) { + // 31 minus leading zeros equals floor(log2(x)) for x > 0 + return 31 - __clz(x); +} + +__device__ __forceinline__ int ceil_log2_32(unsigned int x) { + if (x <= 1) return 0; + return 32 - __clz(x - 1); +} + +template +__global__ void compute_analysis_pass1(int m, + const int* __restrict__ csr_row_ptr, + int* __restrict__ bin_count, + int* __restrict__ row_index_in_bin) +{ + const int tid = threadIdx.x; + const int bid = blockIdx.x; + const int gid = tid + BLOCKSIZE * bid; + + if(gid < m) + { + const int row_length = csr_row_ptr[gid + 1] - csr_row_ptr[gid]; + const int bin = (row_length != 0) ? ceil_log2_32(row_length) : 0; + + row_index_in_bin[gid] = atomicAdd(&bin_count[bin], 1); + } +} + +template +__global__ void compute_analysis_pass2(int m, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ bin_count, + const int* __restrict__ row_index_in_bin, + int* __restrict__ bin_start_ptr, + int* __restrict__ row_index_in_bin_sorted) +{ + const int tid = threadIdx.x; + const int bid = blockIdx.x; + const int gid = tid + BLOCKSIZE * bid; + + __shared__ int bin_start_ptr_shared[32]; + + if(tid == 0) + { + int count = 0; + for(int i = 0; i < 32; i++) + { + const int tmp = bin_count[i]; + bin_start_ptr[i] = count; + bin_start_ptr_shared[i] = count; + count += tmp; + } + } + + __syncthreads(); + + if(gid < m) + { + const int row_length = csr_row_ptr[gid + 1] - csr_row_ptr[gid]; + const int bin = (row_length != 0) ? ceil_log2_32(row_length) : 0; + + row_index_in_bin_sorted[bin_start_ptr_shared[bin] + row_index_in_bin[gid]] = gid; + } +} + +template +__global__ void csrmv_lrb_small_kernel(int m, + int n, + int nnz, + int bin, + int bin_count, + const T alpha, + const int* __restrict__ bin_start_ptr, + const int* __restrict__ row_index_in_bin_sorted, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + const T* __restrict__ csr_val, + const T* __restrict__ x, + const T beta, + T* __restrict__ y) +{ + const int tid = threadIdx.x; + const int bid = blockIdx.x; + const int gid = tid + BLOCKSIZE * bid; + + __shared__ T shared[BLOCKSIZE]; + + if(gid < bin_count) + { + const int row = row_index_in_bin_sorted[bin_start_ptr[bin] + gid]; + + assert(row < m); + + const int start = csr_row_ptr[row]; + const int end = csr_row_ptr[row + 1]; + + T sum = static_cast(0); + for(int j = start; j < end; j++) + { + const int col = csr_col_ind[j]; + const T val = csr_val[j]; + + sum = std::fma(x[col], val, sum); + } + + if(beta == static_cast(0)) + { + y[row] = alpha * sum; + } + else + { + y[row] = std::fma(alpha, sum, beta * y[row]); + } + } +} + +template +__global__ void csrmv_lrb_medium_kernel(int m, + int n, + int nnz, + int bin, + int bin_count, + const T alpha, + const int* __restrict__ bin_start_ptr, + const int* __restrict__ row_index_in_bin_sorted, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + const T* __restrict__ csr_val, + const T* __restrict__ x, + const T beta, + T* __restrict__ y) +{ + const int tid = threadIdx.x; + const int bid = blockIdx.x; + const int gid = tid + BLOCKSIZE * bid; + + const int lid = tid & WARPSIZE - 1; + //const int wid = tid / WARPSIZE; + + for(int i = gid / WARPSIZE; i < bin_count; i += (BLOCKSIZE / WARPSIZE) * gridDim.x) + { + const int row = row_index_in_bin_sorted[bin_start_ptr[bin] + i]; + + const int row_start = csr_row_ptr[row]; + const int row_end = csr_row_ptr[row + 1]; + + T sum = static_cast(0); + for(int j = row_start + lid; j < row_end; j += WARPSIZE) + { + const int col = csr_col_ind[j]; + const T val = csr_val[j]; + + sum = std::fma(x[col], val, sum); + } + + warp_reduction_sum(&sum); + + if(lid == 0) + { + if(beta == static_cast(0)) + { + y[row] = alpha * sum; + } + else + { + y[row] = std::fma(alpha, sum, beta * y[row]); + } + } + } +} + +template +__global__ void csrmv_lrb_medium_large_kernel(int m, + int n, + int nnz, + int bin, + int bin_count, + const T alpha, + const int* __restrict__ bin_start_ptr, + const int* __restrict__ row_index_in_bin_sorted, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + const T* __restrict__ csr_val, + const T* __restrict__ x, + const T beta, + T* __restrict__ y) +{ + const int tid = threadIdx.x; + const int bid = blockIdx.x; + const int gid = tid + BLOCKSIZE * bid; + + __shared__ T shared[BLOCKSIZE]; + + for(int i = bid; i < bin_count; i += gridDim.x) + { + const int row = row_index_in_bin_sorted[bin_start_ptr[bin] + i]; + + const int row_start = csr_row_ptr[row]; + const int row_end = csr_row_ptr[row + 1]; + + T sum = static_cast(0); + for(int j = row_start + tid; j < row_end; j += BLOCKSIZE) + { + const int col = csr_col_ind[j]; + const T val = csr_val[j]; + + sum = std::fma(x[col], val, sum); + } + + shared[tid] = sum; + __syncthreads(); + + block_reduction_sum(shared, tid); + + if(tid == 0) + { + if(beta == static_cast(0)) + { + y[row] = alpha * shared[0]; + } + else + { + y[row] = std::fma(alpha, shared[0], beta * y[row]); + } + } + } +} + +template +__global__ void csrmv_lrb_large_kernel(int m, + int n, + int nnz, + int bin, + int bin_count, + const T alpha, + const int* __restrict__ bin_start_ptr, + const int* __restrict__ row_index_in_bin_sorted, + const int* __restrict__ csr_row_ptr, + const int* __restrict__ csr_col_ind, + const T* __restrict__ csr_val, + const T* __restrict__ x, + const T beta, + T* __restrict__ y) +{ + const int tid = threadIdx.x; + const int bid = blockIdx.x; + const int gid = tid + BLOCKSIZE * bid; + + __shared__ T shared[BLOCKSIZE]; + + const int bin_size = 1 << bin; + + const int blocks_per_row = bin_size / BLOCKSIZE; + + const int row + = row_index_in_bin_sorted[bin_start_ptr[bid / blocks_per_row] + bid % blocks_per_row]; + + const int row_start = csr_row_ptr[row] + BLOCKSIZE * (bid % blocks_per_row); + const int row_end = csr_row_ptr[row + 1]; + + T sum = static_cast(0); + for(int j = row_start + tid; j < row_end; j += BLOCKSIZE * blocks_per_row) + { + const int col = csr_col_ind[j]; + const T val = csr_val[j]; + + sum = std::fma(x[col], val, sum); + } + + shared[tid] = sum; + __syncthreads(); + + block_reduction_sum(shared, tid); + + if(tid == 0) + { + if(beta == static_cast(0)) + { + atomicAdd(&y[row], alpha * shared[0]); + } + else + { + // if(bid / blocks_per_row == 0) + // { + // // atomicAdd(&y[row], alpha * shared[0] + beta); + // } + // else + // { + // atomicAdd(&y[row], alpha * shared[0]); + // } + } + } + + // __shared__ T shared[BLOCKSIZE]; + + // for(int i = bid; i < bin_count; i += gridDim.x) + // { + // const int row = row_index_in_bin_sorted[bin_start_ptr[bin] + i]; + + // const int row_start = csr_row_ptr[row]; + // const int row_end = csr_row_ptr[row + 1]; + + // T sum = static_cast(0); + // for(int j = row_start + tid; j < row_end; j += BLOCKSIZE) + // { + // const int col = csr_col_ind[j]; + // const T val = csr_val[j]; + + // sum = std::fma(x[col], val, sum); + // } + + // shared[tid] = sum; + // __syncthreads(); + + // block_reduction_sum(shared, tid); + + // if(tid == 0) + // { + // if(beta == static_cast(0)) + // { + // y[row] = alpha * shared[0]; + // } + // else + // { + // y[row] = std::fma(alpha, shared[0], beta * y[row]); + // } + // } + // } +} + #endif diff --git a/library/src/backend/device/cuda/cuda_csric0.cu b/library/src/backend/device/cuda/cuda_csric0.cu index 214b9bc..aa652a3 100644 --- a/library/src/backend/device/cuda/cuda_csric0.cu +++ b/library/src/backend/device/cuda/cuda_csric0.cu @@ -74,6 +74,9 @@ void linalg::cuda_csric0_analysis(int m, { std::cout << "cuda_csric0_analysis m: " << m << " n: " << n << " nnz: " << nnz << std::endl; + // Free cuda memory that may have been allocated from previous calls to analysis + free_csric0_cuda_data(descr); + // Free any previous allocations? assert(descr->done_array == nullptr); assert(descr->row_perm == nullptr); diff --git a/library/src/backend/device/cuda/cuda_csrilu0.cu b/library/src/backend/device/cuda/cuda_csrilu0.cu index c3c39c6..9408a24 100644 --- a/library/src/backend/device/cuda/cuda_csrilu0.cu +++ b/library/src/backend/device/cuda/cuda_csrilu0.cu @@ -70,6 +70,9 @@ void linalg::cuda_csrilu0_analysis(int m, const T* csr_val, csrilu0_descr* descr) { + // Free cuda memory that may have been allocated from previous calls to analysis + free_csrilu0_cuda_data(descr); + // Free any previous allocations? assert(descr->done_array == nullptr); assert(descr->row_perm == nullptr); diff --git a/library/src/backend/device/cuda/cuda_csrtrsv.cu b/library/src/backend/device/cuda/cuda_csrtrsv.cu index 964ccc8..dcfd940 100644 --- a/library/src/backend/device/cuda/cuda_csrtrsv.cu +++ b/library/src/backend/device/cuda/cuda_csrtrsv.cu @@ -69,6 +69,9 @@ void linalg::cuda_csrtrsv_analysis(int m, diagonal_type diag_type, csrtrsv_descr* descr) { + // Free cuda memory that may have been allocated from previous calls to analysis + free_csrtrsv_cuda_data(descr); + // Free any previous allocations? assert(descr->done_array == nullptr); assert(descr->row_perm == nullptr); diff --git a/library/src/backend/device/cuda/cuda_matrix_vector.cu b/library/src/backend/device/cuda/cuda_matrix_vector.cu index bfa1ceb..17855cf 100644 --- a/library/src/backend/device/cuda/cuda_matrix_vector.cu +++ b/library/src/backend/device/cuda/cuda_matrix_vector.cu @@ -24,6 +24,9 @@ // //******************************************************************************** +#include +#include +#include #include #include "cuda_matrix_vector.h" @@ -35,44 +38,434 @@ #include "../../../trace.h" -//------------------------------------------------------------------------------- -// sparse matrix-vector product y = A*x -//------------------------------------------------------------------------------- -template -void linalg::cuda_matrix_vector_product(int m, - int n, - int nnz, - const int* csr_row_ptr, - const int* csr_col_ind, - const T* csr_val, - const T* x, - T* y) +namespace linalg { - ROUTINE_TRACE("linalg::cuda_matrix_vector_product_impl"); + static std::string csrmv_alg_to_string(csrmv_algorithm alg) + { + switch(alg) + { + case csrmv_algorithm::default_algorithm: + return "default_algorithm"; + case csrmv_algorithm::rowsplit: + return "rowsplit"; + case csrmv_algorithm::nnzsplit: + return "nnzsplit"; + case csrmv_algorithm::merge_path: + return "marge_path"; + case csrmv_algorithm::lrb: + return "lrb"; + } - int avg_nnz_per_row = nnz / m; + return "invalid"; + } - if(avg_nnz_per_row <= 8) + static void csrmv_analysis_lrb_dispatch( + int m, int n, int nnz, const int* csr_row_ptr, const int* csr_col_ind, csrmv_descr* descr) { - csrmv_vector_kernel<256, 4><<<((m - 1) / (256 / 4) + 1), 256>>>( - m, n, nnz, T(1), csr_row_ptr, csr_col_ind, csr_val, x, T(0), y); + // Free any previous allocations? + assert(descr->bin_start_ptr == nullptr); + assert(descr->row_index_in_bin == nullptr); + assert(descr->row_index_in_bin_sorted == nullptr); + + CHECK_CUDA(cudaMalloc((void**)&(descr->bin_count), sizeof(int) * 32)); + CHECK_CUDA(cudaMalloc((void**)&(descr->bin_start_ptr), sizeof(int) * (32 + 1))); + CHECK_CUDA(cudaMalloc((void**)&(descr->row_index_in_bin), sizeof(int) * m)); + CHECK_CUDA(cudaMalloc((void**)&(descr->row_index_in_bin_sorted), sizeof(int) * m)); + + CHECK_CUDA(cudaMemset(descr->bin_count, 0, sizeof(int) * 32)); + + compute_analysis_pass1<256><<<((m - 1) / 256 + 1), 256>>>( + m, + csr_row_ptr, + descr->bin_count, + descr->row_index_in_bin); + + CHECK_CUDA(cudaMemcpy(descr->hbin_count.data(), + descr->bin_count, + sizeof(int) * 32, + cudaMemcpyDeviceToHost)); + + compute_analysis_pass2<256><<<((m - 1) / 256 + 1), 256>>>( + m, + csr_row_ptr, + descr->bin_count, + descr->row_index_in_bin, + descr->bin_start_ptr, + descr->row_index_in_bin_sorted); + + + + + + + + + // // Free any previous allocations? + // assert(descr->bin_start_ptr == nullptr); + // assert(descr->row_index_in_bin == nullptr); + // assert(descr->row_index_in_bin_sorted == nullptr); + + // CHECK_CUDA(cudaMalloc((void**)&(descr->bin_start_ptr), sizeof(int) * (32 + 1))); + // CHECK_CUDA(cudaMalloc((void**)&(descr->row_index_in_bin), sizeof(int) * m)); + // CHECK_CUDA(cudaMalloc((void**)&(descr->row_index_in_bin_sorted), sizeof(int) * m)); + + // CHECK_CUDA(cudaMemset(descr->row_index_in_bin_sorted, 0, sizeof(int) * m)); + + // for(int i = 0; i < 32; i++) + // { + // descr->hbin_count[i] = 0; + // } + + // std::vector hcsr_row_ptr(m + 1); + // CHECK_CUDA(cudaMemcpy( + // hcsr_row_ptr.data(), csr_row_ptr, sizeof(int) * (m + 1), cudaMemcpyDeviceToHost)); + + // std::vector row_index_in_bin(m, 0); + // std::array bin_start_ptr = {}; + + // for(int i = 0; i < m; i++) + // { + // const int row_length = hcsr_row_ptr[i + 1] - hcsr_row_ptr[i]; + // const int bin = (row_length != 0) ? std::ceil(std::log2(row_length)) : 0; + + // row_index_in_bin[i] = descr->hbin_count[bin]; + // descr->hbin_count[bin]++; + // bin_start_ptr[bin]++; + // } + + // // std::cout << "descr->hbin_count" << std::endl; + // // for(int i = 0; i < 32; i++) + // // { + // // std::cout << descr->hbin_count[i] << " "; + // // } + // // std::cout << "" << std::endl; + + // // std::cout << "bin_start_ptr" << std::endl; + // // for(int i = 0; i < 32 + 1; i++) + // // { + // // std::cout << bin_start_ptr[i] << " "; + // // } + // // std::cout << "" << std::endl; + + // // std::cout << "row_index_in_bin" << std::endl; + // // for(int i = 0; i < m; i++) + // // { + // // std::cout << row_index_in_bin[i] << " "; + // // } + // // std::cout << "" << std::endl; + + // int count = 0; + // for(int i = 0; i < 32; i++) + // { + // const int tmp = bin_start_ptr[i]; + // bin_start_ptr[i] = count; + // count += tmp; + // } + + // // std::cout << "bin_start_ptr" << std::endl; + // // for(int i = 0; i < 32 + 1; i++) + // // { + // // std::cout << bin_start_ptr[i] << " "; + // // } + // // std::cout << "" << std::endl; + + // std::vector row_index_in_bin_sorted(m, 0); + // for(int i = 0; i < m; i++) + // { + // const int row_length = hcsr_row_ptr[i + 1] - hcsr_row_ptr[i]; + // const int bin = (row_length != 0) ? std::ceil(std::log2(row_length)) : 0; + + // row_index_in_bin_sorted[bin_start_ptr[bin] + row_index_in_bin[i]] = i; + // } + + // // std::cout << "bin_start_ptr" << std::endl; + // // for(int i = 0; i < 32 + 1; i++) + // // { + // // std::cout << bin_start_ptr[i] << " "; + // // } + // // std::cout << "" << std::endl; + + // // std::cout << "row_index_in_bin_sorted" << std::endl; + // // for(int i = 0; i < m; i++) + // // { + // // std::cout << row_index_in_bin_sorted[i] << " "; + // // } + // // std::cout << "" << std::endl; + + // CHECK_CUDA(cudaMemcpy(descr->bin_start_ptr, + // bin_start_ptr.data(), + // sizeof(int) * (32 + 1), + // cudaMemcpyHostToDevice)); + // CHECK_CUDA(cudaMemcpy(descr->row_index_in_bin, + // row_index_in_bin.data(), + // sizeof(int) * m, + // cudaMemcpyHostToDevice)); + // CHECK_CUDA(cudaMemcpy(descr->row_index_in_bin_sorted, + // row_index_in_bin_sorted.data(), + // sizeof(int) * m, + // cudaMemcpyHostToDevice)); } - else if(avg_nnz_per_row <= 16) + + static void csrmv_analysis_merge_path_dispatch( + int m, int n, int nnz, const int* csr_row_ptr, const int* csr_col_ind, csrmv_descr* descr) { - csrmv_vector_kernel<256, 8><<<((m - 1) / (256 / 8) + 1), 256>>>( - m, n, nnz, T(1), csr_row_ptr, csr_col_ind, csr_val, x, T(0), y); } - else if(avg_nnz_per_row <= 32) + + static void csrmv_analysis_algorithm_dispatch(int m, + int n, + int nnz, + const int* csr_row_ptr, + const int* csr_col_ind, + csrmv_algorithm alg, + csrmv_descr* descr) { - csrmv_vector_kernel<256, 16><<<((m - 1) / (256 / 16) + 1), 256>>>( - m, n, nnz, T(1), csr_row_ptr, csr_col_ind, csr_val, x, T(0), y); + switch(alg) + { + case csrmv_algorithm::default_algorithm: + case csrmv_algorithm::rowsplit: + case csrmv_algorithm::nnzsplit: + break; + case csrmv_algorithm::merge_path: + csrmv_analysis_merge_path_dispatch(m, n, nnz, csr_row_ptr, csr_col_ind, descr); + break; + case csrmv_algorithm::lrb: + csrmv_analysis_lrb_dispatch(m, n, nnz, csr_row_ptr, csr_col_ind, descr); + break; + default: + throw std::runtime_error("Unknown csrmv_algorithm"); + } } - else + + template + static void csrmv_row_split_dispatch(int m, + int n, + int nnz, + T alpha, + const int* csr_row_ptr, + const int* csr_col_ind, + const T* csr_val, + const T* x, + T beta, + T* y, + const csrmv_descr* descr) { - csrmv_vector_kernel<256, 32><<<((m - 1) / (256 / 32) + 1), 256>>>( - m, n, nnz, T(1), csr_row_ptr, csr_col_ind, csr_val, x, T(0), y); + const int avg_nnz_per_row = nnz / m; + + if(avg_nnz_per_row <= 8) + { + csrmv_row_split_kernel<256, 4><<<((m - 1) / (256 / 4) + 1), 256>>>( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); + } + else if(avg_nnz_per_row <= 16) + { + csrmv_row_split_kernel<256, 8><<<((m - 1) / (256 / 8) + 1), 256>>>( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); + } + else if(avg_nnz_per_row <= 32) + { + csrmv_row_split_kernel<256, 16><<<((m - 1) / (256 / 16) + 1), 256>>>( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); + } + else + { + csrmv_row_split_kernel<256, 32><<<((m - 1) / (256 / 32) + 1), 256>>>( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); + } + } + + template + static void csrmv_nnz_split_dispatch(int m, + int n, + int nnz, + T alpha, + const int* csr_row_ptr, + const int* csr_col_ind, + const T* csr_val, + const T* x, + T beta, + T* y, + const csrmv_descr* descr) + { + CHECK_CUDA(cudaMemset(y, 0, sizeof(T) * m)); // need to call kernel to handle beta + + csrmv_nnz_split_kernel<256, 32, 8><<<((nnz - 1) / (8 * 256) + 1), 256>>>( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); + } + + template + static void csrmv_lrb_dispatch(int m, + int n, + int nnz, + T alpha, + const int* csr_row_ptr, + const int* csr_col_ind, + const T* csr_val, + const T* x, + T beta, + T* y, + const csrmv_descr* descr) + { + // Short rows (bin sizes 2-16) + for(int bin = 0; bin < 5; bin++) + { + if(descr->hbin_count[bin] > 0) + { + // std::cout << "small bin: " << bin << " hbin_count[bin]: " << descr->hbin_count[bin] + // << std::endl; + csrmv_lrb_small_kernel<256> + <<<((descr->hbin_count[bin] - 1) / 256 + 1), 256>>>(m, + n, + nnz, + bin, + descr->hbin_count[bin], + alpha, + descr->bin_start_ptr, + descr->row_index_in_bin_sorted, + csr_row_ptr, + csr_col_ind, + csr_val, + x, + beta, + y); + } + } + + // Medium rows (one warp per row, bin sizes 32-128) + for(int bin = 5; bin < 8; bin++) + { + if(descr->hbin_count[bin] > 0) + { + // int bin_size = 1 << bin; + // std::cout << "medium bin: " << bin << " hbin_count[bin]: " << descr->hbin_count[bin] + // << " bin_size: " << bin_size << std::endl; + csrmv_lrb_medium_kernel<256, 32> + <<<((descr->hbin_count[bin] - 1) / (256 / 32) + 1), 256>>>( + m, + n, + nnz, + bin, + descr->hbin_count[bin], + alpha, + descr->bin_start_ptr, + descr->row_index_in_bin_sorted, + csr_row_ptr, + csr_col_ind, + csr_val, + x, + beta, + y); + } + } + + // Medium-large rows (one block per row, bin sizes 256-2^32) + for(int bin = 8; bin < 32 /*14*/; bin++) + { + if(descr->hbin_count[bin] > 0) + { + // std::cout << "medium-large bin: " << bin + // << " hbin_count[bin]: " << descr->hbin_count[bin] << std::endl; + csrmv_lrb_medium_large_kernel<256> + <<<((descr->hbin_count[bin] - 1) / 256 + 1), 256>>>( + m, + n, + nnz, + bin, + descr->hbin_count[bin], + alpha, + descr->bin_start_ptr, + descr->row_index_in_bin_sorted, + csr_row_ptr, + csr_col_ind, + csr_val, + x, + beta, + y); + } + } + + // for(int bin = 0/*14*/; bin < 32; bin++) + // { + // if(descr->hbin_count[bin] > 0) + // { + // // How many blocks do I need? + // const int bin_size = 1 << bin; + + // const int num_blocks = descr->hbin_count[bin] * (bin_size / 256); + + // csrmv_lrb_large_kernel<256> + // <<>>( + // m, + // n, + // nnz, + // bin, + // descr->hbin_count[bin], + // alpha, + // descr->bin_start_ptr, + // descr->row_index_in_bin_sorted, + // csr_row_ptr, + // csr_col_ind, + // csr_val, + // x, + // beta, + // y); + // } + // } + } + + template + static void csrmv_merge_path_dispatch(int m, + int n, + int nnz, + T alpha, + const int* csr_row_ptr, + const int* csr_col_ind, + const T* csr_val, + const T* x, + T beta, + T* y, + const csrmv_descr* descr) + { + } + + template + static void csrmv_algorithm_dispatch(int m, + int n, + int nnz, + T alpha, + const int* csr_row_ptr, + const int* csr_col_ind, + const T* csr_val, + const T* x, + T beta, + T* y, + csrmv_algorithm alg, + const csrmv_descr* descr) + { + // std::cout << "alg: " << csrmv_alg_to_string(alg) << std::endl; + + switch(alg) + { + case csrmv_algorithm::default_algorithm: + case csrmv_algorithm::rowsplit: + csrmv_row_split_dispatch( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y, descr); + break; + case csrmv_algorithm::nnzsplit: + csrmv_nnz_split_dispatch( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y, descr); + break; + case csrmv_algorithm::merge_path: + csrmv_merge_path_dispatch( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y, descr); + break; + case csrmv_algorithm::lrb: + csrmv_lrb_dispatch( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y, descr); + break; + default: + throw std::runtime_error("Unknown csrmv_algorithm"); + } } - CHECK_CUDA_LAUNCH_ERROR(); } //------------------------------------------------------------------------------- @@ -99,6 +492,21 @@ void linalg::free_csrmv_cuda_data(csrmv_descr* descr) { if(descr != nullptr) { + if(descr->bin_start_ptr != nullptr) + { + CHECK_CUDA(cudaFree(descr->bin_start_ptr)); + descr->bin_start_ptr = nullptr; + } + if(descr->row_index_in_bin != nullptr) + { + CHECK_CUDA(cudaFree(descr->row_index_in_bin)); + descr->row_index_in_bin = nullptr; + } + if(descr->row_index_in_bin_sorted != nullptr) + { + CHECK_CUDA(cudaFree(descr->row_index_in_bin_sorted)); + descr->row_index_in_bin_sorted = nullptr; + } } } @@ -112,6 +520,10 @@ void linalg::cuda_csrmv_analysis(int m, csrmv_algorithm alg, csrmv_descr* descr) { + // Free cuda memory that may have been allocated from previous calls to analysis + free_csrmv_cuda_data(descr); + + csrmv_analysis_algorithm_dispatch(m, n, nnz, csr_row_ptr, csr_col_ind, alg, descr); } template @@ -130,47 +542,10 @@ void linalg::cuda_csrmv_solve(int m, { ROUTINE_TRACE("linalg::cuda_csrmv_solve"); - // const int avg_nnz_per_row = nnz / m; - - // if(avg_nnz_per_row <= 8) - // { - // csrmv_vector_kernel<256, 4><<<((m - 1) / (256 / 4) + 1), 256>>>( - // m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); - // } - // else if(avg_nnz_per_row <= 16) - // { - // csrmv_vector_kernel<256, 8><<<((m - 1) / (256 / 8) + 1), 256>>>( - // m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); - // } - // else if(avg_nnz_per_row <= 32) - // { - // csrmv_vector_kernel<256, 16><<<((m - 1) / (256 / 16) + 1), 256>>>( - // m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); - // } - // else - // { - // csrmv_vector_kernel<256, 32><<<((m - 1) / (256 / 32) + 1), 256>>>( - // m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); - // } - - //int grid_size = ((nnz - 1) / (8 * 256) + 1); - //std::cout << "AAAAAA nnz: " << nnz << " grid_size: " << grid_size << std::endl; - - - - - - CHECK_CUDA(cudaMemset(y, 0, sizeof(T) * m)); // need to call kernel to handle beta - - csrmv_stream_kernel<256, 32, 8><<<((nnz - 1) / (8 * 256) + 1), 256>>>( - m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y); - + csrmv_algorithm_dispatch( + m, n, nnz, alpha, csr_row_ptr, csr_col_ind, csr_val, x, beta, y, alg, descr); } -template void linalg::cuda_matrix_vector_product( - int, int, int, const int*, const int*, const double*, const double*, double*); -template void linalg::cuda_matrix_vector_product( - int, int, int, const int*, const int*, const float*, const float*, float*); template void linalg::cuda_compute_residual( int, int, int, const int*, const int*, const double*, const double*, const double*, double*); template void linalg::cuda_compute_residual( diff --git a/library/src/backend/device/cuda/cuda_matrix_vector.h b/library/src/backend/device/cuda/cuda_matrix_vector.h index e55cb86..1c8a8be 100644 --- a/library/src/backend/device/cuda/cuda_matrix_vector.h +++ b/library/src/backend/device/cuda/cuda_matrix_vector.h @@ -31,15 +31,6 @@ namespace linalg { template - void cuda_matrix_vector_product(int m, - int n, - int nnz, - const int* csr_row_ptr, - const int* csr_col_ind, - const T* csr_val, - const T* x, - T* y); - template void cuda_compute_residual(int m, int n, int nnz, diff --git a/library/src/backend/host/host_matrix_vector.cpp b/library/src/backend/host/host_matrix_vector.cpp index ec829a7..e1d8af4 100644 --- a/library/src/backend/host/host_matrix_vector.cpp +++ b/library/src/backend/host/host_matrix_vector.cpp @@ -121,6 +121,7 @@ void linalg::host_csrmv_solve(double alpha, case csrmv_algorithm::merge_path: case csrmv_algorithm::rowsplit: case csrmv_algorithm::nnzsplit: + case csrmv_algorithm::lrb: host_csrmv_impl(A.get_m(), A.get_n(), A.get_nnz(), diff --git a/library/src/csr_matrix.cpp b/library/src/csr_matrix.cpp index 81a2a77..b08df0c 100644 --- a/library/src/csr_matrix.cpp +++ b/library/src/csr_matrix.cpp @@ -49,6 +49,11 @@ csr_matrix::csr_matrix() , nnz(0) , on_host(true) { + create_csrmv_descr(&this->descr_mv); + create_csrgemm_descr(&this->descr_gemm); + create_csrtrsv_descr(&this->descr_sv); + create_csric0_descr(&this->descr_ic); + create_csrilu0_descr(&this->descr_ilu); } template @@ -71,11 +76,22 @@ csr_matrix::csr_matrix(const std::vector& csr_row_ptr, this->nnz = nnz; this->on_host = true; + + create_csrmv_descr(&this->descr_mv); + create_csrgemm_descr(&this->descr_gemm); + create_csrtrsv_descr(&this->descr_sv); + create_csric0_descr(&this->descr_ic); + create_csrilu0_descr(&this->descr_ilu); } template csr_matrix::~csr_matrix() { + destroy_csrmv_descr(this->descr_mv); + destroy_csrgemm_descr(this->descr_gemm); + destroy_csrtrsv_descr(this->descr_sv); + destroy_csric0_descr(this->descr_ic); + destroy_csrilu0_descr(this->descr_ilu); } template @@ -304,14 +320,16 @@ void csr_matrix::multiply_by_vector(vector& y, const vector& x) const { ROUTINE_TRACE("csr_matrix::multiply_by_vector"); - csrmv_descr* descr = nullptr; - create_csrmv_descr(&descr); + // csrmv_descr* descr = nullptr; + // create_csrmv_descr(&descr); - csrmv_analysis(*this, csrmv_algorithm::default_algorithm, descr); + csrmv_analysis(*this, csrmv_algorithm::lrb, descr_mv); + // csrmv_analysis(*this, csrmv_algorithm::default_algorithm, descr); - csrmv_solve(1.0, *this, x, 0.0, y, csrmv_algorithm::default_algorithm, descr); + csrmv_solve(1.0, *this, x, 0.0, y, csrmv_algorithm::lrb, descr_mv); + // csrmv_solve(1.0, *this, x, 0.0, y, csrmv_algorithm::default_algorithm, descr); - destroy_csrmv_descr(descr); + // destroy_csrmv_descr(descr); } template @@ -319,14 +337,14 @@ void csr_matrix::multiply_by_vector_and_add(vector& y, const vector& x) { ROUTINE_TRACE("csr_matrix::multiply_by_vector_and_add"); - csrmv_descr* descr = nullptr; - create_csrmv_descr(&descr); + // csrmv_descr* descr = nullptr; + // create_csrmv_descr(&descr); - csrmv_analysis(*this, csrmv_algorithm::default_algorithm, descr); + csrmv_analysis(*this, csrmv_algorithm::default_algorithm, descr_mv); - csrmv_solve(1.0, *this, x, 1.0, y, csrmv_algorithm::default_algorithm, descr); + csrmv_solve(1.0, *this, x, 1.0, y, csrmv_algorithm::default_algorithm, descr_mv); - destroy_csrmv_descr(descr); + // destroy_csrmv_descr(descr); } template @@ -342,14 +360,14 @@ void csr_matrix::multiply_by_matrix(csr_matrix& C, const csr_matrix& B) D.move_to_device(); } - csrgemm_descr* descr = nullptr; - create_csrgemm_descr(&descr); + // csrgemm_descr* descr = nullptr; + // create_csrgemm_descr(&descr); - csrgemm_nnz(*this, B, D, C, csrgemm_algorithm::default_algorithm, descr); + csrgemm_nnz(*this, B, D, C, csrgemm_algorithm::default_algorithm, descr_gemm); - csrgemm_solve(1.0, *this, B, 0.0, D, C, csrgemm_algorithm::default_algorithm, descr); + csrgemm_solve(1.0, *this, B, 0.0, D, C, csrgemm_algorithm::default_algorithm, descr_gemm); - destroy_csrgemm_descr(descr); + // destroy_csrgemm_descr(descr); } template @@ -357,13 +375,13 @@ void csr_matrix::triangular_solve_lower(vector& x, const vector& y, boo { ROUTINE_TRACE("csr_matrix::triangular_solve_lower"); - csrtrsv_descr* descr = nullptr; - create_csrtrsv_descr(&descr); + // csrtrsv_descr* descr = nullptr; + // create_csrtrsv_descr(&descr); csrtrsv_analysis(*this, triangular_type::lower, unit_diag ? diagonal_type::unit : diagonal_type::non_unit, - descr); + descr_sv); // Perform the triangular solve csrtrsv_solve(*this, @@ -372,9 +390,9 @@ void csr_matrix::triangular_solve_lower(vector& x, const vector& y, boo static_cast(1.0), triangular_type::lower, unit_diag ? diagonal_type::unit : diagonal_type::non_unit, - descr); + descr_sv); - destroy_csrtrsv_descr(descr); + // destroy_csrtrsv_descr(descr); } template @@ -382,13 +400,13 @@ void csr_matrix::triangular_solve_upper(vector& x, const vector& y, boo { ROUTINE_TRACE("csr_matrix::triangular_solve_upper"); - csrtrsv_descr* descr = nullptr; - create_csrtrsv_descr(&descr); + // csrtrsv_descr* descr = nullptr; + // create_csrtrsv_descr(&descr); csrtrsv_analysis(*this, triangular_type::upper, unit_diag ? diagonal_type::unit : diagonal_type::non_unit, - descr); + descr_sv); // Perform the triangular solve csrtrsv_solve(*this, @@ -397,9 +415,9 @@ void csr_matrix::triangular_solve_upper(vector& x, const vector& y, boo static_cast(1.0), triangular_type::upper, unit_diag ? diagonal_type::unit : diagonal_type::non_unit, - descr); + descr_sv); - destroy_csrtrsv_descr(descr); + // destroy_csrtrsv_descr(descr); } template @@ -407,15 +425,15 @@ void csr_matrix::compute_incomplete_cholesky_factorization() { ROUTINE_TRACE("csr_matrix::compute_incomplete_cholesky_factorization"); - csric0_descr* descr = nullptr; - create_csric0_descr(&descr); + // csric0_descr* descr = nullptr; + // create_csric0_descr(&descr); - csric0_analysis(*this, descr); + csric0_analysis(*this, descr_ic); // Compute Cholesky factorization inplace - csric0_compute(*this, descr); + csric0_compute(*this, descr_ic); - destroy_csric0_descr(descr); + // destroy_csric0_descr(descr); } template @@ -423,15 +441,15 @@ void csr_matrix::compute_incomplete_LU_factorization() { ROUTINE_TRACE("csr_matrix::compute_incomplete_LU_factorization"); - csrilu0_descr* descr = nullptr; - create_csrilu0_descr(&descr); + // csrilu0_descr* descr = nullptr; + // create_csrilu0_descr(&descr); - csrilu0_analysis(*this, descr); + csrilu0_analysis(*this, descr_ilu); // Compute ILU factorization inplace - csrilu0_compute(*this, descr); + csrilu0_compute(*this, descr_ilu); - destroy_csrilu0_descr(descr); + // destroy_csrilu0_descr(descr); } template diff --git a/library/src/descriptors/csrmv_descr_internal.h b/library/src/descriptors/csrmv_descr_internal.h index 087c07c..ba58c3c 100644 --- a/library/src/descriptors/csrmv_descr_internal.h +++ b/library/src/descriptors/csrmv_descr_internal.h @@ -27,10 +27,19 @@ #ifndef CSRMV_DESCR_INTERNAL_H #define CSRMV_DESCR_INTERNAL_H +#include + namespace linalg { struct csrmv_descr { + // LRB algorithm + std::array hbin_count; // how many rows belong to each bin + + int* bin_count; // desvice array of size 32, how many rows belong to each bin + int* bin_start_ptr; // device array of size (32 + 1) + int* row_index_in_bin; // device array of size m row_indices_perm? + int* row_index_in_bin_sorted; // device array of size m row_indices? }; } diff --git a/library/src/linalg_math.cpp b/library/src/linalg_math.cpp index c28c150..622e2ee 100644 --- a/library/src/linalg_math.cpp +++ b/library/src/linalg_math.cpp @@ -206,7 +206,10 @@ void linalg::create_csrmv_descr(csrmv_descr** descr) { ROUTINE_TRACE("linalg::create_csrmv_descr"); - *descr = new csrmv_descr; + *descr = new csrmv_descr; + (*descr)->bin_start_ptr = nullptr; + (*descr)->row_index_in_bin = nullptr; + (*descr)->row_index_in_bin_sorted = nullptr; } void linalg::destroy_csrmv_descr(csrmv_descr* descr) diff --git a/tools/Benchmark_Metrics_Tools.md b/tools/Benchmark_Metrics_Tools.md new file mode 100644 index 0000000..e0670a2 --- /dev/null +++ b/tools/Benchmark_Metrics_Tools.md @@ -0,0 +1,115 @@ +# GoogleTest Benchmark Metrics Tools + +A collection of Python scripts to run GoogleTest matrix benchmarks, parse performance output into structured YAML files, and visualize effective bandwidth. + +--- + +## Prerequisites + +Install required Python dependencies before using the scripts: + +```bash +pip install pyyaml matplotlib numpy +``` + +--- + +## Tool Overview + +### 1. `run_metrics.py` (Benchmark Runner & Parser) + +Runs your C++ GoogleTest executable with a specified filter, extracts performance metrics (`Solve time` and `Effective Bandwidth`), normalizes matrix identifiers (e.g., `quick_ci_GPU_nos1_mtx` -> `nos1.mtx`), and exports the results to a timestamped YAML file. + +**Usage:** + +```bash +python run_metrics.py --filter [--exe ] +``` + +**Options:** +* `--filter` *(Required)*: The test filter substring (e.g., `matrix_vector_product`). +* `--exe` *(Optional)*: Path to test executable (Default: `./test_main.exe`). + +**Example:** + +```bash +python run_metrics.py --filter matrix_vector_product --exe ./test_main.exe +``` + +**Output:** +Generates a file named `_.yaml` containing structured metrics: + +```yaml +- matrix_file: nos1.mtx + solve_time_ms: 7.6949 + effective_bandwidth_gbps: 0.22025 +- matrix_file: nos2.mtx + solve_time_ms: 7.4423 + effective_bandwidth_gbps: 0.924284 +``` + +--- + +### 2. `plot_metrics.py` (Single Run Visualization) + +Reads a generated benchmark YAML file and creates a bar chart displaying effective bandwidth (GB/s) per matrix file, complete with exact numerical values rendered above each bar. + +**Usage:** + +```bash +python plot_metrics.py [--output ] +``` + +**Options:** +* `YAML_FILE` *(Required)*: Path to the parsed benchmark YAML file. +* `-o, --output` *(Optional)*: Target path for the output PNG (Default: replaces `.yaml` extension with `.png`). + +**Example:** + +```bash +python plot_metrics.py matrix_vector_product_20260830_141922.yaml -o bandwidth_single_run.png +``` + +--- + +### 3. `compare_metrics.py` (Dual Run Comparison) + +Compares two benchmark YAML runs (e.g., Baseline vs. Optimized, or GPU vs. CPU) on a side-by-side grouped bar chart. Automatically aligns matching matrices and handles missing entries cleanly. + +**Usage:** + +```bash +python compare_metrics.py [--label1 ] [--label2 ] [--output ] +``` + +**Options:** +* `FILE1` *(Required)*: Path to the first benchmark YAML file. +* `FILE2` *(Required)*: Path to the second benchmark YAML file. +* `--label1` *(Optional)*: Legend label for Dataset 1 (Default: filename of Dataset 1). +* `--label2` *(Optional)*: Legend label for Dataset 2 (Default: filename of Dataset 2). +* `-o, --output` *(Optional)*: Target path for the comparison plot PNG. + +**Example:** + +```bash +python compare_metrics.py run_baseline.yaml run_optimized.yaml --label1 "Baseline Kernel" --label2 "Optimized Kernel" -o kernel_comparison.png +``` + +--- + +## End-to-End Workflow Example + +1. **Execute benchmarks and parse output:** + ```bash + python run_metrics.py --filter matrix_vector_product --exe ./bin/sparse_tests.exe + ``` + +2. **Visualize single benchmark run:** + ```bash + python plot_metrics.py matrix_vector_product_20260830_141922.yaml + ``` + +3. **Compare two separate optimization runs:** + ```bash + python compare_metrics.py matrix_vector_product_20260830_100000.yaml matrix_vector_product_20260830_141922.yaml --label1 "v1.0" --label2 "v1.1" -o release_comparison.png + ``` diff --git a/tools/comparison_matrix_vector_product_20260830_145014_vs_matrix_vector_product_20260830_145014.png b/tools/comparison_matrix_vector_product_20260830_145014_vs_matrix_vector_product_20260830_145014.png new file mode 100644 index 0000000..4aabb3d Binary files /dev/null and b/tools/comparison_matrix_vector_product_20260830_145014_vs_matrix_vector_product_20260830_145014.png differ diff --git a/tools/plot_bandwidth.py b/tools/plot_bandwidth.py new file mode 100644 index 0000000..5d6e575 --- /dev/null +++ b/tools/plot_bandwidth.py @@ -0,0 +1,62 @@ +import argparse +from pathlib import Path +import yaml +import matplotlib.pyplot as plt + +def plot_bandwidth(yaml_path, output_image=None): + with open(yaml_path, 'r') as f: + data = yaml.safe_load(f) + + if not data: + print(f"No data found in {yaml_path}") + return + + matrices = [item['matrix_file'] for item in data] + bandwidths = [item['effective_bandwidth_gbps'] for item in data] + + fig, ax = plt.subplots(figsize=(10, 6)) + bars = ax.bar(matrices, bandwidths, color='#1f77b4', edgecolor='#003366', width=0.6) + + ax.set_xlabel("Matrix File", fontsize=11, fontweight='bold', labelpad=10) + ax.set_ylabel("Effective Bandwidth (GB/s)", fontsize=11, fontweight='bold') + ax.set_title(f"Bandwidth Comparison: {Path(yaml_path).stem}", fontsize=13, fontweight='bold', pad=15) + + ax.grid(axis='y', linestyle='--', alpha=0.6) + ax.set_axisbelow(True) + + plt.xticks(rotation=45, ha='right', fontsize=10) + plt.yticks(fontsize=10) + + # Add numeric labels above each bar + for bar in bars: + height = bar.get_height() + ax.annotate( + f"{height:.2f}", + xy=(bar.get_x() + bar.get_width() / 2, height), + xytext=(0, 3), + textcoords="offset points", + ha='center', + va='bottom', + fontsize=9 + ) + + plt.tight_layout() + + # Save output plot + if not output_image: + output_image = Path(yaml_path).with_suffix('.png') + + plt.savefig(output_image, dpi=300) + print(f"Plot successfully saved to: {output_image}") + plt.show() + +def main(): + parser = argparse.ArgumentParser(description="Plot matrix effective bandwidth from YAML benchmark output.") + parser.add_argument("yaml_file", help="Path to the input YAML file") + parser.add_argument("--output", "-o", help="Path for saving the output PNG image (optional)") + args = parser.parse_args() + + plot_bandwidth(args.yaml_file, args.output) + +if __name__ == "__main__": + main() diff --git a/tools/plot_bandwidth_compare.py b/tools/plot_bandwidth_compare.py new file mode 100644 index 0000000..d9408e0 --- /dev/null +++ b/tools/plot_bandwidth_compare.py @@ -0,0 +1,90 @@ +import argparse +from pathlib import Path +import yaml +import numpy as np +import matplotlib.pyplot as plt + +def load_yaml_data(filepath): + """Loads YAML file and returns a dictionary mapping matrix_file -> effective_bandwidth_gbps.""" + with open(filepath, 'r') as f: + data = yaml.safe_load(f) + if not data: + return {} + return {item['matrix_file']: item['effective_bandwidth_gbps'] for item in data} + +def plot_comparison(file1_path, file2_path, label1=None, label2=None, output_image=None): + data1 = load_yaml_data(file1_path) + data2 = load_yaml_data(file2_path) + + if not data1 and not data2: + print("No valid data found in either YAML file.") + return + + # Use file names as legend labels if custom labels aren't provided + label1 = label1 or Path(file1_path).stem + label2 = label2 or Path(file2_path).stem + + # Preserve order from file 1 and append any new matrices found in file 2 + all_matrices = list(data1.keys()) + for matrix in data2.keys(): + if matrix not in all_matrices: + all_matrices.append(matrix) + + bw1 = [data1.get(m, 0.0) for m in all_matrices] + bw2 = [data2.get(m, 0.0) for m in all_matrices] + + x = np.arange(len(all_matrices)) + width = 0.35 + + fig, ax = plt.subplots(figsize=(12, 6)) + rects1 = ax.bar(x - width/2, bw1, width, label=label1, color='#1f77b4', edgecolor='#003366') + rects2 = ax.bar(x + width/2, bw2, width, label=label2, color='#ff7f0e', edgecolor='#b33c00') + + ax.set_xlabel("Matrix File", fontsize=11, fontweight='bold', labelpad=10) + ax.set_ylabel("Effective Bandwidth (GB/s)", fontsize=11, fontweight='bold') + ax.set_title(f"Bandwidth Comparison: {label1} vs {label2}", fontsize=13, fontweight='bold', pad=15) + + ax.set_xticks(x) + ax.set_xticklabels(all_matrices, rotation=45, ha='right', fontsize=10) + ax.legend(fontsize=10) + ax.grid(axis='y', linestyle='--', alpha=0.6) + ax.set_axisbelow(True) + + # Add numerical values above each bar + def annotate_bars(rects): + for rect in rects: + height = rect.get_height() + if height > 0: + ax.annotate( + f"{height:.2f}", + xy=(rect.get_x() + rect.get_width() / 2, height), + xytext=(0, 3), + textcoords="offset points", + ha='center', va='bottom', fontsize=8 + ) + + annotate_bars(rects1) + annotate_bars(rects2) + + plt.tight_layout() + + if not output_image: + output_image = f"comparison_{label1}_vs_{label2}.png" + + plt.savefig(output_image, dpi=300) + print(f"Comparison plot saved to: {output_image}") + plt.show() + +def main(): + parser = argparse.ArgumentParser(description="Compare matrix bandwidth metrics from two YAML benchmark outputs.") + parser.add_argument("file1", help="Path to the first YAML benchmark file") + parser.add_argument("file2", help="Path to the second YAML benchmark file") + parser.add_argument("--label1", help="Custom legend label for dataset 1 (e.g., 'GPU Run 1')") + parser.add_argument("--label2", help="Custom legend label for dataset 2 (e.g., 'GPU Run 2')") + parser.add_argument("--output", "-o", help="Path for saving output plot PNG") + args = parser.parse_args() + + plot_comparison(args.file1, args.file2, args.label1, args.label2, args.output) + +if __name__ == "__main__": + main() diff --git a/tools/record_test_performance.py b/tools/record_test_performance.py new file mode 100644 index 0000000..c3d7341 --- /dev/null +++ b/tools/record_test_performance.py @@ -0,0 +1,92 @@ +import subprocess +import re +import yaml +import argparse +from datetime import datetime + +def extract_matrix_name(param_string): + """ + Extracts the matrix filename by looking after 'GPU_' or 'CPU_' + and converting the trailing '_mtx' to '.mtx'. + + Example: 'quick_ci_GPU_nos1_mtx' -> 'nos1.mtx' + """ + match = re.search(r'(?:GPU|CPU)_(.*)', param_string) + if match: + extracted = match.group(1) + # Replace the final trailing '_mtx' with '.mtx' + if extracted.endswith('_mtx'): + return extracted[:-4] + '.mtx' + return extracted + return param_string + +def run_tests_and_parse(executable, test_filter): + command = [executable, f"--gtest_filter=*{test_filter}*"] + print(f"Running: {' '.join(command)}") + + process = subprocess.run(command, capture_output=True, text=True) + + results = [] + current_test = None + + re_run = re.compile(r'\[\s*RUN\s*\]\s+(.*)') + re_time = re.compile(r'Solve time:\s*([0-9.]+)\s*ms') + re_bw = re.compile(r'Effective Bandwidth:\s*([0-9.]+)\s*GB/s') + re_end = re.compile(r'\[\s*(OK|FAILED)\s*\]') + + for line in process.stdout.splitlines(): + run_match = re_run.match(line) + if run_match: + full_test_name = run_match.group(1) + raw_param = full_test_name.split('/')[-1] + + # Clean up raw_param to get 'nos1.mtx', 'nos2.mtx', etc. + matrix_filename = extract_matrix_name(raw_param) + + current_test = { + 'matrix_file': matrix_filename, + 'solve_time_ms': None, + 'effective_bandwidth_gbps': None + } + continue + + if current_test: + time_match = re_time.search(line) + if time_match: + current_test['solve_time_ms'] = float(time_match.group(1)) + + bw_match = re_bw.search(line) + if bw_match: + current_test['effective_bandwidth_gbps'] = float(bw_match.group(1)) + + end_match = re_end.match(line) + if end_match: + if current_test['solve_time_ms'] is not None and current_test['effective_bandwidth_gbps'] is not None: + results.append(current_test) + current_test = None + + return results + +def main(): + parser = argparse.ArgumentParser(description="Parse GoogleTest performance metrics into YAML.") + parser.add_argument("--exe", default="./test_main.exe", help="Path to the googletest executable") + parser.add_argument("--filter", required=True, help="Test filter string (e.g., matrix_vector_product)") + args = parser.parse_args() + + results = run_tests_and_parse(args.exe, args.filter) + + if not results: + print("No test metrics found. Verify your executable path and filter.") + return + + datestamp = datetime.now().strftime("%Y%m%d_%H%M%S") + output_filename = f"{args.filter}_{datestamp}.yaml" + + with open(output_filename, 'w') as f: + yaml.dump(results, f, sort_keys=False, default_flow_style=False) + + print(f"Successfully parsed {len(results)} tests.") + print(f"Results saved to: {output_filename}") + +if __name__ == "__main__": + main()