From a417e0b50716199008a5026c61030780f4d69321 Mon Sep 17 00:00:00 2001 From: govindchari Date: Tue, 26 May 2026 11:39:43 -0700 Subject: [PATCH 1/2] Batch solve --- CMakeLists.txt | 14 +- algebra/cuda/cuda_linalg.cu | 22 +- algebra/cuda/cudss_backend.cu | 933 +++++++++++++++++++++++++++++- algebra/cuda/cudss_backend.h | 20 +- batch_cudss_benchmark_findings.md | 163 ++++++ benchmarks/pdg_batch_compare.cpp | 341 +++++++++++ include/qoco.h | 3 +- include/qoco_batch.h | 90 +++ include/structs.h | 23 +- src/qoco_batch.c | 164 ++++++ tests/unit_tests/batch_test.cpp | 128 ++++ 11 files changed, 1871 insertions(+), 30 deletions(-) create mode 100644 batch_cudss_benchmark_findings.md create mode 100644 benchmarks/pdg_batch_compare.cpp create mode 100644 include/qoco_batch.h create mode 100644 src/qoco_batch.c create mode 100644 tests/unit_tests/batch_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 84a21cec..531dc24f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -108,6 +108,7 @@ if(${MATLAB}) endif() set_property(GLOBAL APPEND PROPERTY QOCO_SOURCES "${CMAKE_CURRENT_SOURCE_DIR}/src/qoco_api.c" + ${CMAKE_CURRENT_SOURCE_DIR}/src/qoco_batch.c ${CMAKE_CURRENT_SOURCE_DIR}/src/input_validation.c ${CMAKE_CURRENT_SOURCE_DIR}/src/common_linalg.c ${CMAKE_CURRENT_SOURCE_DIR}/src/kkt.c @@ -118,6 +119,7 @@ set_property(GLOBAL APPEND PROPERTY QOCO_SOURCES "${CMAKE_CURRENT_SOURCE_DIR}/sr set_property(GLOBAL APPEND PROPERTY QOCO_HEADERS "${CMAKE_CURRENT_SOURCE_DIR}/include/qoco.h" ${CMAKE_CURRENT_SOURCE_DIR}/include/qoco_api.h + ${CMAKE_CURRENT_SOURCE_DIR}/include/qoco_batch.h ${CMAKE_CURRENT_SOURCE_DIR}/include/input_validation.h ${CMAKE_CURRENT_SOURCE_DIR}/include/qoco_linalg.h ${CMAKE_CURRENT_SOURCE_DIR}/include/kkt.h @@ -190,6 +192,11 @@ if(IS_LINUX OR IS_MACOS) endif() target_include_directories(qoco PUBLIC ${qoco_include}) +if(${QOCO_ALGEBRA_BACKEND} STREQUAL "cuda") + target_compile_definitions(qoco PUBLIC QOCO_ALGEBRA_BACKEND_CUDA) +else() + target_compile_definitions(qoco PUBLIC QOCO_ALGEBRA_BACKEND_BUILTIN) +endif() target_sources(qoco PRIVATE ${qoco_sources}) # Build qoco static library. @@ -205,6 +212,11 @@ if(IS_LINUX OR IS_MACOS) endif() target_include_directories(qocostatic PUBLIC ${qoco_include}) +if(${QOCO_ALGEBRA_BACKEND} STREQUAL "cuda") + target_compile_definitions(qocostatic PUBLIC QOCO_ALGEBRA_BACKEND_CUDA) +else() + target_compile_definitions(qocostatic PUBLIC QOCO_ALGEBRA_BACKEND_BUILTIN) +endif() target_sources(qocostatic PRIVATE ${qoco_sources}) # Build qoco demo. @@ -251,4 +263,4 @@ install( ) install(FILES ${qoco_headers} DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/qoco" -) \ No newline at end of file +) diff --git a/algebra/cuda/cuda_linalg.cu b/algebra/cuda/cuda_linalg.cu index 73ed8c10..51d0e9a7 100644 --- a/algebra/cuda/cuda_linalg.cu +++ b/algebra/cuda/cuda_linalg.cu @@ -352,15 +352,19 @@ void sync_vector_to_device(QOCOVectorf* v) void sync_matrix_to_device(QOCOMatrix* M) { if (M->csc && M->d_csc) { - CUDA_CHECK(cudaMemcpy(M->d_csc_host->x, M->csc->x, - M->csc->nnz * sizeof(QOCOFloat), - cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(M->d_csc_host->i, M->csc->i, - M->csc->nnz * sizeof(QOCOInt), - cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(M->d_csc_host->p, M->csc->p, - (M->csc->n + 1) * sizeof(QOCOInt), - cudaMemcpyHostToDevice)); + if (M->csc->nnz > 0) { + CUDA_CHECK(cudaMemcpy(M->d_csc_host->x, M->csc->x, + M->csc->nnz * sizeof(QOCOFloat), + cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(M->d_csc_host->i, M->csc->i, + M->csc->nnz * sizeof(QOCOInt), + cudaMemcpyHostToDevice)); + } + if (M->csc->p && M->d_csc_host->p) { + CUDA_CHECK(cudaMemcpy(M->d_csc_host->p, M->csc->p, + (M->csc->n + 1) * sizeof(QOCOInt), + cudaMemcpyHostToDevice)); + } } } diff --git a/algebra/cuda/cudss_backend.cu b/algebra/cuda/cudss_backend.cu index 2ae9a564..e12a3086 100644 --- a/algebra/cuda/cudss_backend.cu +++ b/algebra/cuda/cudss_backend.cu @@ -9,8 +9,19 @@ */ #include "cudss_backend.h" +#ifdef __cplusplus +extern "C" { +#endif +#include "equilibration.h" +#include "input_validation.h" +#include "qoco_utils.h" +#ifdef __cplusplus +} +#endif #include #include +#include +#include #define CUDA_CHECK(call) \ do { \ @@ -47,6 +58,86 @@ static void* g_cusparse_handle = NULL; static void* g_cublas_handle = NULL; static int g_libs_loaded = 0; +typedef struct { + double serial_factor_sec; + double serial_solve_sec; + long long serial_factor_calls; + long long serial_solve_calls; + double batch_factor_sec; + double batch_solve_sec; + long long batch_factor_calls; + long long batch_solve_calls; +} QOCOCudaLinsysTiming; + +static QOCOCudaLinsysTiming g_linsys_timing = {0}; +static int g_linsys_timing_enabled = 0; + +static double qoco_cuda_now_sec(void) +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec / 1e9; +} + +extern "C" void qoco_cuda_linsys_timing_reset(void) +{ + memset(&g_linsys_timing, 0, sizeof(g_linsys_timing)); + g_linsys_timing_enabled = 1; +} + +extern "C" void qoco_cuda_linsys_timing_set_enabled(int enabled) +{ + g_linsys_timing_enabled = enabled != 0; +} + +extern "C" void qoco_cuda_linsys_timing_get(QOCOCudaLinsysTiming* timing) +{ + if (timing) { + *timing = g_linsys_timing; + } +} + +static void qoco_cuda_timed_cudss_execute( + cudssHandle_t handle, cudssPhase_t phase, cudssConfig_t config, + cudssData_t data, cudssMatrix_t matA, cudssMatrix_t matB, + cudssMatrix_t matC, unsigned char is_batch) +{ + if (!g_linsys_timing_enabled || + (phase != CUDSS_PHASE_FACTORIZATION && phase != CUDSS_PHASE_SOLVE)) { + CUDSS_CHECK(g_cuda_funcs.cudssExecute(handle, phase, config, data, matA, + matB, matC)); + return; + } + + CUDA_CHECK(cudaDeviceSynchronize()); + const double start_sec = qoco_cuda_now_sec(); + CUDSS_CHECK(g_cuda_funcs.cudssExecute(handle, phase, config, data, matA, matB, + matC)); + CUDA_CHECK(cudaDeviceSynchronize()); + const double elapsed_sec = qoco_cuda_now_sec() - start_sec; + + if (is_batch) { + if (phase == CUDSS_PHASE_FACTORIZATION) { + g_linsys_timing.batch_factor_sec += elapsed_sec; + g_linsys_timing.batch_factor_calls++; + } + else { + g_linsys_timing.batch_solve_sec += elapsed_sec; + g_linsys_timing.batch_solve_calls++; + } + } + else { + if (phase == CUDSS_PHASE_FACTORIZATION) { + g_linsys_timing.serial_factor_sec += elapsed_sec; + g_linsys_timing.serial_factor_calls++; + } + else { + g_linsys_timing.serial_solve_sec += elapsed_sec; + g_linsys_timing.serial_solve_calls++; + } + } +} + // Global accessor for function pointers (for use in cuda_linalg.cu) CudaLibFuncs* get_cuda_funcs(void) { return &g_cuda_funcs; } @@ -109,14 +200,26 @@ int load_cuda_libraries(void) g_cuda_funcs.cudssMatrixCreateCsr = (typeof(g_cuda_funcs.cudssMatrixCreateCsr))dlsym(g_cudss_handle, "cudssMatrixCreateCsr"); + g_cuda_funcs.cudssMatrixCreateBatchCsr = + (typeof(g_cuda_funcs.cudssMatrixCreateBatchCsr))dlsym( + g_cudss_handle, "cudssMatrixCreateBatchCsr"); g_cuda_funcs.cudssExecute = (typeof(g_cuda_funcs.cudssExecute))dlsym(g_cudss_handle, "cudssExecute"); g_cuda_funcs.cudssMatrixCreateDn = (typeof(g_cuda_funcs.cudssMatrixCreateDn))dlsym(g_cudss_handle, "cudssMatrixCreateDn"); + g_cuda_funcs.cudssMatrixCreateBatchDn = + (typeof(g_cuda_funcs.cudssMatrixCreateBatchDn))dlsym( + g_cudss_handle, "cudssMatrixCreateBatchDn"); g_cuda_funcs.cudssMatrixSetValues = (typeof(g_cuda_funcs.cudssMatrixSetValues))dlsym(g_cudss_handle, "cudssMatrixSetValues"); + g_cuda_funcs.cudssMatrixSetBatchValues = + (typeof(g_cuda_funcs.cudssMatrixSetBatchValues))dlsym( + g_cudss_handle, "cudssMatrixSetBatchValues"); + g_cuda_funcs.cudssMatrixSetBatchCsrPointers = + (typeof(g_cuda_funcs.cudssMatrixSetBatchCsrPointers))dlsym( + g_cudss_handle, "cudssMatrixSetBatchCsrPointers"); g_cuda_funcs.cudssMatrixDestroy = (typeof(g_cuda_funcs.cudssMatrixDestroy))dlsym(g_cudss_handle, "cudssMatrixDestroy"); @@ -297,6 +400,43 @@ struct LinSysData { QOCOInt* d_GttoKKTcsr; }; +struct CudaBatchLinSysData { + QOCOInt batch_count; + QOCOInt Kn; + QOCOInt nnz; + + cudssHandle_t handle; + cudssConfig_t config; + cudssData_t data; + + cudssMatrix_t K_batch; + cudssMatrix_t rhs_batch; + cudssMatrix_t xyz_batch; + + QOCOInt* d_csr_row_ptr; + QOCOInt* d_csr_col_ind; + + void** h_csr_row_ptrs; + void** h_csr_col_inds; + void** h_csr_values; + void** h_rhs_values; + void** h_xyz_values; + + void** d_csr_row_ptrs; + void** d_csr_col_inds; + void** d_csr_values; + void** d_rhs_values; + void** d_xyz_values; + + QOCOInt* h_nrows; + QOCOInt* h_ncols; + QOCOInt* h_nnz; + QOCOInt* h_nrhs; + QOCOInt* h_ld; +}; + +extern "C" void qoco_cuda_batch_cleanup(QOCOBatchSolver* batch); + // Convert CSC to CSR on CPU and copy to GPU static void csc_to_csr_device(const QOCOCscMatrix* csc, QOCOInt** csr_row_ptr, QOCOInt** csr_col_ind, QOCOFloat** csr_val, @@ -519,6 +659,16 @@ static LinSysData* cudss_setup(QOCOProblemData* data, QOCOSettings* settings, valueType_setup, CUDSS_MTYPE_SYMMETRIC, CUDSS_MVIEW_UPPER, CUDSS_BASE_ZERO)); + // Create dense matrix wrappers for solution and RHS vectors (column vectors). + CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateDn( + &linsys_data->d_rhs_matrix, (int64_t)linsys_data->Kn, 1, + (int64_t)linsys_data->Kn, linsys_data->d_rhs_matrix_data, valueType_setup, + CUDSS_LAYOUT_COL_MAJOR)); + CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateDn( + &linsys_data->d_xyz_matrix, (int64_t)linsys_data->Kn, 1, + (int64_t)linsys_data->Kn, linsys_data->d_xyz_matrix_data, valueType_setup, + CUDSS_LAYOUT_COL_MAJOR)); + // Run analysis phase. CUDSS_CHECK(g_cuda_funcs.cudssExecute( linsys_data->handle, CUDSS_PHASE_ANALYSIS, linsys_data->config, @@ -582,18 +732,6 @@ static LinSysData* cudss_setup(QOCOProblemData* data, QOCOSettings* settings, qoco_free(h_csr_col_ind); qoco_free(csc2csr); - // Create dense matrix wrappers for solution and RHS vectors (column vectors) - // Note: d_rhs_matrix wraps d_rhs_matrix_data, d_xyz_matrix wraps - // d_xyz_matrix_data - CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateDn( - &linsys_data->d_rhs_matrix, (int64_t)linsys_data->Kn, 1, - (int64_t)linsys_data->Kn, linsys_data->d_rhs_matrix_data, valueType_setup, - CUDSS_LAYOUT_COL_MAJOR)); - CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateDn( - &linsys_data->d_xyz_matrix, (int64_t)linsys_data->Kn, 1, - (int64_t)linsys_data->Kn, linsys_data->d_xyz_matrix_data, valueType_setup, - CUDSS_LAYOUT_COL_MAJOR)); - return linsys_data; } @@ -664,10 +802,10 @@ static void cudss_factor(LinSysData* linsys_data, QOCOInt n, (void)n; (void)kkt_dynamic_reg; - CUDSS_CHECK(g_cuda_funcs.cudssExecute( + qoco_cuda_timed_cudss_execute( linsys_data->handle, CUDSS_PHASE_FACTORIZATION, linsys_data->config, linsys_data->data, linsys_data->K_csr, linsys_data->d_xyz_matrix, - linsys_data->d_rhs_matrix)); + linsys_data->d_rhs_matrix, 0); } static void cudss_solve_system(LinSysData* linsys_data, const QOCOFloat* rhs, @@ -680,10 +818,10 @@ static void cudss_solve_system(LinSysData* linsys_data, const QOCOFloat* rhs, CUDA_CHECK(cudaMemset(linsys_data->d_xyz_matrix_data, 0, linsys_data->Kn * sizeof(QOCOFloat))); - CUDSS_CHECK(g_cuda_funcs.cudssExecute( + qoco_cuda_timed_cudss_execute( linsys_data->handle, CUDSS_PHASE_SOLVE, linsys_data->config, linsys_data->data, linsys_data->K_csr, linsys_data->d_xyz_matrix, - linsys_data->d_rhs_matrix)); + linsys_data->d_rhs_matrix, 0); if (sol != linsys_data->d_xyz_matrix_data) { CUDA_CHECK(cudaMemcpy(sol, linsys_data->d_xyz_matrix_data, @@ -947,6 +1085,769 @@ static void cudss_cleanup(LinSysData* linsys_data) qoco_free(linsys_data); } +static unsigned char cuda_batch_symbols_available(void) +{ + return g_cuda_funcs.cudssMatrixCreateBatchCsr && + g_cuda_funcs.cudssMatrixCreateBatchDn; +} + +static void cuda_batch_refresh_value_pointers(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data) +{ + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + LinSysData* linsys_data = (LinSysData*)batch->solvers[item]->linsys_data; + batch_data->h_csr_values[item] = (void*)linsys_data->d_csr_val; + batch_data->h_rhs_values[item] = (void*)linsys_data->d_rhs_matrix_data; + batch_data->h_xyz_values[item] = (void*)linsys_data->d_xyz_matrix_data; + } + + CUDA_CHECK(cudaMemcpy(batch_data->d_csr_values, batch_data->h_csr_values, + batch_data->batch_count * sizeof(void*), + cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(batch_data->d_rhs_values, batch_data->h_rhs_values, + batch_data->batch_count * sizeof(void*), + cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(batch_data->d_xyz_values, batch_data->h_xyz_values, + batch_data->batch_count * sizeof(void*), + cudaMemcpyHostToDevice)); +} + +static QOCOInt cuda_batch_validate_solvers(QOCOBatchSolver* batch) +{ + QOCOSolver* first = batch->solvers[0]; + QOCOProblemData* first_data = first->work->data; + QOCOInt first_Pnnz = get_nnz(first_data->P); + QOCOInt first_Annz = get_nnz(first_data->A); + QOCOInt first_Gnnz = get_nnz(first_data->G); + + for (QOCOInt item = 0; item < batch->batch_count; ++item) { + QOCOSolver* solver = batch->solvers[item]; + if (!solver || !solver->work || !solver->work->data || + !solver->linsys_data) { + return QOCO_DATA_VALIDATION_ERROR; + } + + QOCOProblemData* data = solver->work->data; + if (data->n != first_data->n || data->m != first_data->m || + data->p != first_data->p || data->l != first_data->l || + data->nsoc != first_data->nsoc || + solver->work->Wnnz != first->work->Wnnz || + get_nnz(data->P) != first_Pnnz || get_nnz(data->A) != first_Annz || + get_nnz(data->G) != first_Gnnz) { + return QOCO_DATA_VALIDATION_ERROR; + } + } + return QOCO_NO_ERROR; +} + +static QOCOCscMatrix* cuda_batch_construct_reference_kkt(QOCOSolver* solver) +{ + QOCOProblemData* data = solver->work->data; + QOCOSettings* settings = solver->settings; + QOCOCscMatrix* Kcsc = NULL; + + set_cpu_mode(1); + Kcsc = construct_kkt( + get_csc_matrix(data->P), get_csc_matrix(data->A), get_csc_matrix(data->G), + get_csc_matrix(data->At), get_csc_matrix(data->Gt), + settings->kkt_static_reg_A, data->n, data->m, data->p, data->l, + data->nsoc, get_data_vectori(data->q), NULL, NULL, NULL, NULL, NULL, + solver->work->Wnnz); + set_cpu_mode(0); + return Kcsc; +} + +extern "C" QOCOInt qoco_cuda_batch_setup(QOCOBatchSolver* batch) +{ + if (!batch || !batch->solvers || batch->batch_count <= 0) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + + if (!load_cuda_libraries()) { + fprintf(stderr, "Failed to load CUDA libraries\n"); + return QOCO_SETUP_ERROR; + } + + if (!cuda_batch_symbols_available()) { + fprintf(stderr, + "Loaded cuDSS does not provide cudssMatrixCreateBatchCsr and " + "cudssMatrixCreateBatchDn\n"); + return QOCO_SETUP_ERROR; + } + + QOCOInt validation = cuda_batch_validate_solvers(batch); + if (validation != QOCO_NO_ERROR) { + return qoco_error((enum qoco_error_code)validation); + } + + qoco_cuda_batch_cleanup(batch); + + CudaBatchLinSysData* batch_data = + (CudaBatchLinSysData*)qoco_calloc(1, sizeof(CudaBatchLinSysData)); + if (!batch_data) { + return qoco_error(QOCO_MALLOC_ERROR); + } + + QOCOSolver* first = batch->solvers[0]; + QOCOProblemData* first_problem = first->work->data; + batch_data->batch_count = batch->batch_count; + batch_data->Kn = first_problem->n + first_problem->m + first_problem->p; + + QOCOCscMatrix* Kcsc = cuda_batch_construct_reference_kkt(first); + if (!Kcsc) { + qoco_free(batch_data); + return QOCO_SETUP_ERROR; + } + batch_data->nnz = Kcsc->nnz; + + QOCOFloat* unused_csr_val = NULL; + csc_to_csr_device(Kcsc, &batch_data->d_csr_row_ptr, + &batch_data->d_csr_col_ind, &unused_csr_val, NULL, NULL, + NULL); + cudaFree(unused_csr_val); + free_qoco_csc_matrix(Kcsc); + + batch_data->h_csr_row_ptrs = + (void**)qoco_malloc(batch_data->batch_count * sizeof(void*)); + batch_data->h_csr_col_inds = + (void**)qoco_malloc(batch_data->batch_count * sizeof(void*)); + batch_data->h_csr_values = + (void**)qoco_malloc(batch_data->batch_count * sizeof(void*)); + batch_data->h_rhs_values = + (void**)qoco_malloc(batch_data->batch_count * sizeof(void*)); + batch_data->h_xyz_values = + (void**)qoco_malloc(batch_data->batch_count * sizeof(void*)); + batch_data->h_nrows = + (QOCOInt*)qoco_malloc(batch_data->batch_count * sizeof(QOCOInt)); + batch_data->h_ncols = + (QOCOInt*)qoco_malloc(batch_data->batch_count * sizeof(QOCOInt)); + batch_data->h_nnz = + (QOCOInt*)qoco_malloc(batch_data->batch_count * sizeof(QOCOInt)); + batch_data->h_nrhs = + (QOCOInt*)qoco_malloc(batch_data->batch_count * sizeof(QOCOInt)); + batch_data->h_ld = + (QOCOInt*)qoco_malloc(batch_data->batch_count * sizeof(QOCOInt)); + + if (!batch_data->h_csr_row_ptrs || !batch_data->h_csr_col_inds || + !batch_data->h_csr_values || !batch_data->h_rhs_values || + !batch_data->h_xyz_values || !batch_data->h_nrows || + !batch_data->h_ncols || !batch_data->h_nnz || !batch_data->h_nrhs || + !batch_data->h_ld) { + batch->batch_linsys_data = batch_data; + qoco_cuda_batch_cleanup(batch); + return qoco_error(QOCO_MALLOC_ERROR); + } + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + batch_data->h_csr_row_ptrs[item] = (void*)batch_data->d_csr_row_ptr; + batch_data->h_csr_col_inds[item] = (void*)batch_data->d_csr_col_ind; + batch_data->h_nrows[item] = batch_data->Kn; + batch_data->h_ncols[item] = batch_data->Kn; + batch_data->h_nnz[item] = batch_data->nnz; + batch_data->h_nrhs[item] = 1; + batch_data->h_ld[item] = batch_data->Kn; + } + + CUDA_CHECK(cudaMalloc(&batch_data->d_csr_row_ptrs, + batch_data->batch_count * sizeof(void*))); + CUDA_CHECK(cudaMalloc(&batch_data->d_csr_col_inds, + batch_data->batch_count * sizeof(void*))); + CUDA_CHECK(cudaMalloc(&batch_data->d_csr_values, + batch_data->batch_count * sizeof(void*))); + CUDA_CHECK(cudaMalloc(&batch_data->d_rhs_values, + batch_data->batch_count * sizeof(void*))); + CUDA_CHECK(cudaMalloc(&batch_data->d_xyz_values, + batch_data->batch_count * sizeof(void*))); + + CUDA_CHECK(cudaMemcpy(batch_data->d_csr_row_ptrs, + batch_data->h_csr_row_ptrs, + batch_data->batch_count * sizeof(void*), + cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(batch_data->d_csr_col_inds, + batch_data->h_csr_col_inds, + batch_data->batch_count * sizeof(void*), + cudaMemcpyHostToDevice)); + + cuda_batch_refresh_value_pointers(batch, batch_data); + + CUDSS_CHECK(g_cuda_funcs.cudssCreate(&batch_data->handle)); + CUDSS_CHECK(g_cuda_funcs.cudssConfigCreate(&batch_data->config)); + CUDSS_CHECK( + g_cuda_funcs.cudssDataCreate(batch_data->handle, &batch_data->data)); + int value = 0; + CUDSS_CHECK(g_cuda_funcs.cudssConfigSet(batch_data->config, + CUDSS_CONFIG_USE_SUPERPANELS, + (void*)&value, sizeof(int))); + + cudaDataType_t value_type = + (sizeof(QOCOFloat) == 8) ? CUDA_R_64F : CUDA_R_32F; + + CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateBatchDn( + &batch_data->rhs_batch, (int64_t)batch_data->batch_count, + batch_data->h_nrows, batch_data->h_nrhs, batch_data->h_ld, + batch_data->d_rhs_values, CUDA_R_32I, value_type, + CUDSS_LAYOUT_COL_MAJOR)); + CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateBatchDn( + &batch_data->xyz_batch, (int64_t)batch_data->batch_count, + batch_data->h_nrows, batch_data->h_nrhs, batch_data->h_ld, + batch_data->d_xyz_values, CUDA_R_32I, value_type, + CUDSS_LAYOUT_COL_MAJOR)); + CUDSS_CHECK(g_cuda_funcs.cudssMatrixCreateBatchCsr( + &batch_data->K_batch, (int64_t)batch_data->batch_count, + batch_data->h_nrows, batch_data->h_ncols, batch_data->h_nnz, + batch_data->d_csr_row_ptrs, NULL, batch_data->d_csr_col_inds, + batch_data->d_csr_values, CUDA_R_32I, value_type, + CUDSS_MTYPE_SYMMETRIC, CUDSS_MVIEW_UPPER, CUDSS_BASE_ZERO)); + + CUDSS_CHECK(g_cuda_funcs.cudssExecute( + batch_data->handle, CUDSS_PHASE_ANALYSIS, batch_data->config, + batch_data->data, batch_data->K_batch, batch_data->xyz_batch, + batch_data->rhs_batch)); + + batch->batch_linsys_data = batch_data; + batch->batch_linsys_stale = 0; + return QOCO_NO_ERROR; +} + +extern "C" void qoco_cuda_batch_cleanup(QOCOBatchSolver* batch) +{ + if (!batch || !batch->batch_linsys_data) { + return; + } + + CudaBatchLinSysData* batch_data = + (CudaBatchLinSysData*)batch->batch_linsys_data; + + if (g_libs_loaded) { + if (batch_data->K_batch) { + g_cuda_funcs.cudssMatrixDestroy(batch_data->K_batch); + } + if (batch_data->rhs_batch) { + g_cuda_funcs.cudssMatrixDestroy(batch_data->rhs_batch); + } + if (batch_data->xyz_batch) { + g_cuda_funcs.cudssMatrixDestroy(batch_data->xyz_batch); + } + if (batch_data->data) { + g_cuda_funcs.cudssDataDestroy(batch_data->handle, batch_data->data); + } + if (batch_data->config) { + g_cuda_funcs.cudssConfigDestroy(batch_data->config); + } + if (batch_data->handle) { + g_cuda_funcs.cudssDestroy(batch_data->handle); + } + } + + cudaFree(batch_data->d_csr_row_ptr); + cudaFree(batch_data->d_csr_col_ind); + cudaFree(batch_data->d_csr_row_ptrs); + cudaFree(batch_data->d_csr_col_inds); + cudaFree(batch_data->d_csr_values); + cudaFree(batch_data->d_rhs_values); + cudaFree(batch_data->d_xyz_values); + + qoco_free(batch_data->h_csr_row_ptrs); + qoco_free(batch_data->h_csr_col_inds); + qoco_free(batch_data->h_csr_values); + qoco_free(batch_data->h_rhs_values); + qoco_free(batch_data->h_xyz_values); + qoco_free(batch_data->h_nrows); + qoco_free(batch_data->h_ncols); + qoco_free(batch_data->h_nnz); + qoco_free(batch_data->h_nrhs); + qoco_free(batch_data->h_ld); + qoco_free(batch_data); + + batch->batch_linsys_data = NULL; + batch->batch_linsys_stale = 1; +} + +static void cuda_batch_set_matrix_values(CudaBatchLinSysData* batch_data) +{ + if (g_cuda_funcs.cudssMatrixSetBatchValues) { + CUDSS_CHECK(g_cuda_funcs.cudssMatrixSetBatchValues( + batch_data->K_batch, batch_data->d_csr_values)); + } +} + +static void cuda_batch_set_nt_identity(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active) +{ + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (active && !active[item]) { + continue; + } + QOCOSolver* solver = batch->solvers[item]; + cudss_set_nt_identity((LinSysData*)solver->linsys_data, + solver->work->data->m); + } + cuda_batch_set_matrix_values(batch_data); +} + +static void cuda_batch_update_nt(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active) +{ + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (active && !active[item]) { + continue; + } + QOCOSolver* solver = batch->solvers[item]; + cudss_update_nt((LinSysData*)solver->linsys_data, solver->work->WtW, + solver->settings->kkt_static_reg_G, + solver->work->data->m); + } + cuda_batch_set_matrix_values(batch_data); +} + +static void cuda_batch_factor(CudaBatchLinSysData* batch_data) +{ + qoco_cuda_timed_cudss_execute( + batch_data->handle, CUDSS_PHASE_FACTORIZATION, batch_data->config, + batch_data->data, batch_data->K_batch, batch_data->xyz_batch, + batch_data->rhs_batch, 1); +} + +static void cuda_batch_execute_solve(CudaBatchLinSysData* batch_data) +{ + qoco_cuda_timed_cudss_execute( + batch_data->handle, CUDSS_PHASE_SOLVE, batch_data->config, + batch_data->data, batch_data->K_batch, batch_data->xyz_batch, + batch_data->rhs_batch, 1); +} + +static void cuda_batch_pack_rhs(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active, + unsigned char use_residual) +{ + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + LinSysData* linsys_data = (LinSysData*)batch->solvers[item]->linsys_data; + if (active && !active[item]) { + CUDA_CHECK(cudaMemset(linsys_data->d_rhs_matrix_data, 0, + batch_data->Kn * sizeof(QOCOFloat))); + } + else { + QOCOFloat* rhs = + use_residual ? linsys_data->d_xyz_matrix_data + : get_data_vectorf(batch->solvers[item]->work->rhs); + CUDA_CHECK(cudaMemcpy(linsys_data->d_rhs_matrix_data, rhs, + batch_data->Kn * sizeof(QOCOFloat), + cudaMemcpyDeviceToDevice)); + } + CUDA_CHECK(cudaMemset(linsys_data->d_xyz_matrix_data, 0, + batch_data->Kn * sizeof(QOCOFloat))); + } +} + +static void cuda_batch_copy_solve_output(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active) +{ + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (active && !active[item]) { + continue; + } + LinSysData* linsys_data = (LinSysData*)batch->solvers[item]->linsys_data; + QOCOFloat* x = get_data_vectorf(batch->solvers[item]->work->xyz); + CUDA_CHECK(cudaMemcpy(x, linsys_data->d_xyz_matrix_data, + batch_data->Kn * sizeof(QOCOFloat), + cudaMemcpyDeviceToDevice)); + } +} + +static unsigned char cuda_batch_any_active(const unsigned char* active, + QOCOInt n) +{ + for (QOCOInt i = 0; i < n; ++i) { + if (active[i]) { + return 1; + } + } + return 0; +} + +static void cuda_batch_solve_refined(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active) +{ + cuda_batch_pack_rhs(batch, batch_data, active, 0); + cuda_batch_execute_solve(batch_data); + cuda_batch_copy_solve_output(batch, batch_data, active); + + unsigned char* refine_active = + (unsigned char*)qoco_calloc(batch_data->batch_count, sizeof(unsigned char)); + QOCOInt* ir_used = + (QOCOInt*)qoco_calloc(batch_data->batch_count, sizeof(QOCOInt)); + QOCOFloat* best_res = + (QOCOFloat*)qoco_calloc(batch_data->batch_count, sizeof(QOCOFloat)); + if (!refine_active || !ir_used || !best_res) { + qoco_free(refine_active); + qoco_free(ir_used); + qoco_free(best_res); + return; + } + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (active && !active[item]) { + continue; + } + + QOCOSolver* solver = batch->solvers[item]; + LinSysData* linsys_data = (LinSysData*)solver->linsys_data; + QOCOWorkspace* work = solver->work; + QOCOFloat* x = get_data_vectorf(work->xyz); + QOCOFloat* b = get_data_vectorf(work->rhs); + QOCOFloat* best_sol = get_data_vectorf(work->xyzbuff1); + + best_res[item] = + compute_linsys_residual(linsys_data, work, b, x, + linsys_data->d_xyz_matrix_data); + copy_arrayf(x, best_sol, batch_data->Kn); + refine_active[item] = + (best_res[item] >= solver->settings->ir_tol && + solver->settings->max_ir_iters > 0); + } + + while (cuda_batch_any_active(refine_active, batch_data->batch_count)) { + cuda_batch_pack_rhs(batch, batch_data, refine_active, 1); + cuda_batch_execute_solve(batch_data); + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (!refine_active[item]) { + continue; + } + + QOCOSolver* solver = batch->solvers[item]; + LinSysData* linsys_data = (LinSysData*)solver->linsys_data; + QOCOWorkspace* work = solver->work; + QOCOFloat* x = get_data_vectorf(work->xyz); + QOCOFloat* b = get_data_vectorf(work->rhs); + QOCOFloat* best_sol = get_data_vectorf(work->xyzbuff1); + + qoco_axpy(linsys_data->d_xyz_matrix_data, x, x, 1.0, batch_data->Kn); + QOCOFloat new_res = + compute_linsys_residual(linsys_data, work, b, x, + linsys_data->d_xyz_matrix_data); + + if (new_res >= best_res[item]) { + copy_arrayf(best_sol, x, batch_data->Kn); + refine_active[item] = 0; + continue; + } + + ir_used[item]++; + best_res[item] = new_res; + copy_arrayf(x, best_sol, batch_data->Kn); + refine_active[item] = + (new_res >= solver->settings->ir_tol && + ir_used[item] < solver->settings->max_ir_iters); + } + } + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (!active || active[item]) { + batch->solvers[item]->work->ir_iters += ir_used[item]; + } + } + + qoco_free(refine_active); + qoco_free(ir_used); + qoco_free(best_res); +} + +static void cuda_batch_finish_item(QOCOSolver* solver) +{ + QOCOWorkspace* work = solver->work; + stop_timer(&(work->solve_timer)); + + unsigned char restored = 0; + if (solver->sol->status == QOCO_NUMERICAL_ERROR) { + restored = restore_best_iterate(solver); + } + unscale_variables(work); + copy_solution(solver); + if (solver->settings->verbose) { + if (restored) { + printf("Best iterate (%d) restored\n", work->best_iter); + } + print_footer(solver->sol, + (enum qoco_solve_status)solver->sol->status); + } +} + +static void cuda_batch_initialize_ipm(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + unsigned char* active) +{ + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + QOCOSolver* solver = batch->solvers[item]; + QOCOWorkspace* work = solver->work; + QOCOProblemData* data = work->data; + + set_Wfull_identity(work->Wfull, work->Wnnzfull, work->Wsoc_idx, data); + work->a = 1.0; + + QOCOFloat* rhs = get_data_vectorf(work->rhs); + copy_and_negate_arrayf(get_data_vectorf(data->c), rhs, data->n); + copy_arrayf(get_data_vectorf(data->b), &rhs[data->n], data->p); + copy_arrayf(get_data_vectorf(data->h), &rhs[data->n + data->p], data->m); + } + + cuda_batch_set_nt_identity(batch, batch_data, active); + cuda_batch_factor(batch_data); + cuda_batch_solve_refined(batch, batch_data, active); + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + QOCOSolver* solver = batch->solvers[item]; + QOCOWorkspace* work = solver->work; + QOCOProblemData* data = work->data; + QOCOFloat* xyz = get_data_vectorf(work->xyz); + + copy_arrayf(xyz, get_data_vectorf(work->x), data->n); + copy_arrayf(&xyz[data->n], get_data_vectorf(work->y), data->p); + copy_arrayf(&xyz[data->n + data->p], get_data_vectorf(work->z), data->m); + copy_and_negate_arrayf(&xyz[data->n + data->p], get_data_vectorf(work->s), + data->m); + + bring2cone(get_data_vectorf(work->s), get_data_vectori(work->soc_idx), + data); + bring2cone(get_data_vectorf(work->z), get_data_vectori(work->soc_idx), + data); + + if (work->use_x0) { + QOCOFloat* x0 = get_data_vectorf(work->x0); + QOCOFloat* x = get_data_vectorf(work->x); + QOCOFloat* Dinvruiz = get_data_vectorf(work->scaling->Dinvruiz); + ew_product(x0, Dinvruiz, x, data->n); + + if (data->m > 0) { + QOCOFloat* s = get_data_vectorf(work->s); + QOCOFloat* h = get_data_vectorf(data->h); + SpMv(data->G, x, s); + qoco_axpy(s, h, s, -1.0, data->m); + bring2cone(s, get_data_vectori(work->soc_idx), data); + } + } + } +} + +static void cuda_batch_after_affine_solve(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active) +{ + (void)batch_data; + for (QOCOInt item = 0; item < batch->batch_count; ++item) { + if (active && !active[item]) { + continue; + } + QOCOSolver* solver = batch->solvers[item]; + QOCOWorkspace* work = solver->work; + QOCOProblemData* data = work->data; + + QOCOFloat* Wfull = get_data_vectorf(work->Wfull); + QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx); + QOCOInt* soc_idx = get_data_vectori(work->soc_idx); + QOCOFloat* lambda = get_data_vectorf(work->lambda); + QOCOFloat* Ds = get_data_vectorf(work->Ds); + QOCOFloat* ubuff1 = get_data_vectorf(work->ubuff1); + QOCOInt* q = get_data_vectori(data->q); + QOCOFloat* xyz = get_data_vectorf(work->xyz); + QOCOFloat* Dzaff = &xyz[data->n + data->p]; + + nt_multiply(Wfull, Wsoc_idx, soc_idx, Dzaff, ubuff1, data->l, data->m, + data->nsoc, q); + copy_and_negate_arrayf(ubuff1, ubuff1, data->m); + qoco_axpy(lambda, ubuff1, ubuff1, -1.0, data->m); + nt_multiply(Wfull, Wsoc_idx, soc_idx, ubuff1, Ds, data->l, data->m, + data->nsoc, q); + + compute_centering(solver); + construct_kkt_comb_rhs(work); + } +} + +static void cuda_batch_after_combined_solve(QOCOBatchSolver* batch, + CudaBatchLinSysData* batch_data, + const unsigned char* active, + QOCOInt iter) +{ + (void)batch_data; + for (QOCOInt item = 0; item < batch->batch_count; ++item) { + if (active && !active[item]) { + continue; + } + + QOCOSolver* solver = batch->solvers[item]; + QOCOWorkspace* work = solver->work; + QOCOProblemData* data = work->data; + + if (check_nan(work->xyz)) { + work->a = 0.0; + solver->sol->iters = iter; + solver->sol->ir_iters += work->ir_iters; + continue; + } + + QOCOFloat* Wfull = get_data_vectorf(work->Wfull); + QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx); + QOCOInt* soc_idx = get_data_vectori(work->soc_idx); + QOCOFloat* lambda = get_data_vectorf(work->lambda); + QOCOFloat* Ds = get_data_vectorf(work->Ds); + QOCOFloat* ubuff1 = get_data_vectorf(work->ubuff1); + QOCOFloat* ubuff2 = get_data_vectorf(work->ubuff2); + QOCOFloat* ubuff3 = get_data_vectorf(work->ubuff3); + QOCOInt* q = get_data_vectori(data->q); + QOCOFloat* xyz = get_data_vectorf(work->xyz); + QOCOFloat* Dz = &xyz[data->n + data->p]; + + cone_division(lambda, Ds, ubuff1, data->l, data->nsoc, q, soc_idx); + nt_multiply(Wfull, Wsoc_idx, soc_idx, Dz, ubuff2, data->l, data->m, + data->nsoc, q); + qoco_axpy(ubuff2, ubuff1, ubuff3, -1.0, data->m); + nt_multiply(Wfull, Wsoc_idx, soc_idx, ubuff3, Ds, data->l, data->m, + data->nsoc, q); + + QOCOFloat a = + qoco_min(linesearch(get_data_vectorf(work->s), Ds, 0.99, solver), + linesearch(get_data_vectorf(work->z), Dz, 0.99, solver)); + work->a = a; + + QOCOFloat* Dx = xyz; + QOCOFloat* Dy = &xyz[data->n]; + qoco_axpy(Dx, get_data_vectorf(work->x), get_data_vectorf(work->x), a, + data->n); + qoco_axpy(Ds, get_data_vectorf(work->s), get_data_vectorf(work->s), a, + data->m); + qoco_axpy(Dy, get_data_vectorf(work->y), get_data_vectorf(work->y), a, + data->p); + qoco_axpy(Dz, get_data_vectorf(work->z), get_data_vectorf(work->z), a, + data->m); + + solver->sol->iters = iter; + solver->sol->ir_iters += work->ir_iters; + if (solver->settings->verbose) { + log_iter(solver); + } + } +} + +extern "C" QOCOInt qoco_cuda_batch_solve(QOCOBatchSolver* batch) +{ + if (!batch || !batch->solvers || !batch->statuses) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + + if (!batch->batch_linsys_data || batch->batch_linsys_stale) { + QOCOInt exit = qoco_cuda_batch_setup(batch); + if (exit != QOCO_NO_ERROR) { + return exit; + } + } + + CudaBatchLinSysData* batch_data = + (CudaBatchLinSysData*)batch->batch_linsys_data; + cuda_batch_refresh_value_pointers(batch, batch_data); + + unsigned char* active = + (unsigned char*)qoco_calloc(batch_data->batch_count, sizeof(unsigned char)); + if (!active) { + return qoco_error(QOCO_MALLOC_ERROR); + } + + QOCOInt max_iters = 0; + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + QOCOSolver* solver = batch->solvers[item]; + if (qoco_validate_settings(solver->settings)) { + qoco_free(active); + return qoco_error(QOCO_SETTINGS_VALIDATION_ERROR); + } + active[item] = 1; + batch->statuses[item] = QOCO_UNSOLVED; + solver->sol->status = QOCO_UNSOLVED; + max_iters = qoco_max(max_iters, solver->settings->max_iters); + start_timer(&(solver->work->solve_timer)); + if (solver->settings->verbose) { + print_header(solver); + } + } + + log_ipm_iter(0); + cuda_batch_initialize_ipm(batch, batch_data, active); + + for (QOCOInt iter = 1; iter <= max_iters; ++iter) { + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (!active[item]) { + continue; + } + + QOCOSolver* solver = batch->solvers[item]; + QOCOWorkspace* work = solver->work; + QOCOProblemData* data = work->data; + + compute_kkt_residual(data, work->x, work->y, work->s, work->z, + work->kktres, solver->settings->kkt_static_reg_P, + work->xyzbuff1, work->xbuff, work->ubuff1); + solver->sol->obj = compute_objective( + data, work->x, work->xbuff, solver->settings->kkt_static_reg_P, + work->scaling->k); + work->mu = compute_mu(work->s, work->z, data->m); + + if (check_stopping(solver)) { + cuda_batch_finish_item(solver); + batch->statuses[item] = solver->sol->status; + active[item] = 0; + } + } + + if (!cuda_batch_any_active(active, batch_data->batch_count)) { + qoco_free(active); + return QOCO_NO_ERROR; + } + + log_ipm_iter(iter); + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (!active[item]) { + continue; + } + QOCOSolver* solver = batch->solvers[item]; + compute_nt_scaling(solver->work); + solver->work->ir_iters = 0; + } + + cuda_batch_update_nt(batch, batch_data, active); + cuda_batch_factor(batch_data); + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (active[item]) { + construct_kkt_aff_rhs(batch->solvers[item]->work); + } + } + cuda_batch_solve_refined(batch, batch_data, active); + + cuda_batch_after_affine_solve(batch, batch_data, active); + cuda_batch_solve_refined(batch, batch_data, active); + cuda_batch_after_combined_solve(batch, batch_data, active, iter); + } + + for (QOCOInt item = 0; item < batch_data->batch_count; ++item) { + if (!active[item]) { + continue; + } + QOCOSolver* solver = batch->solvers[item]; + solver->sol->status = QOCO_MAX_ITER; + restore_best_iterate(solver); + cuda_batch_finish_item(solver); + batch->statuses[item] = solver->sol->status; + } + + qoco_free(active); + return QOCO_NO_ERROR; +} + static const char* cudss_name() { return "cuda/cuDSS"; } LinSysBackend backend = {.linsys_name = cudss_name, diff --git a/algebra/cuda/cudss_backend.h b/algebra/cuda/cudss_backend.h index 8d17289f..9b52f726 100644 --- a/algebra/cuda/cudss_backend.h +++ b/algebra/cuda/cudss_backend.h @@ -41,15 +41,31 @@ typedef struct { cudssMatrix_t* mat, int64_t rows, int64_t cols, int64_t nnz, void* csrRowPtr, void* csrRowPtr_type, void* csrColInd, void* csrVal, cudaDataType_t idxType, cudaDataType_t valType, cudssMatrixType_t type, - int view, cudssIndexBase_t base); + cudssMatrixViewType_t view, cudssIndexBase_t base); + cudssStatus_t (*cudssMatrixCreateBatchCsr)( + cudssMatrix_t* mat, int64_t batchCount, void* nrows, void* ncols, + void* nnz, void** rowStart, void** rowEnd, void** colIndices, + void** values, cudaDataType_t indexType, cudaDataType_t valueType, + cudssMatrixType_t mtype, cudssMatrixViewType_t mview, + cudssIndexBase_t indexBase); cudssStatus_t (*cudssExecute)(cudssHandle_t handle, cudssPhase_t phase, cudssConfig_t config, cudssData_t data, cudssMatrix_t matA, cudssMatrix_t matB, cudssMatrix_t matC); cudssStatus_t (*cudssMatrixCreateDn)(cudssMatrix_t* mat, int64_t rows, int64_t cols, int64_t ld, void* data, - cudaDataType_t type, int layout); + cudaDataType_t type, + cudssLayout_t layout); + cudssStatus_t (*cudssMatrixCreateBatchDn)( + cudssMatrix_t* mat, int64_t batchCount, void* nrows, void* ncols, + void* ld, void** values, cudaDataType_t indexType, + cudaDataType_t valueType, cudssLayout_t layout); cudssStatus_t (*cudssMatrixSetValues)(cudssMatrix_t mat, void* data); + cudssStatus_t (*cudssMatrixSetBatchValues)(cudssMatrix_t mat, + void** values); + cudssStatus_t (*cudssMatrixSetBatchCsrPointers)( + cudssMatrix_t mat, void** rowOffsets, void** rowEnd, void** colIndices, + void** values); cudssStatus_t (*cudssMatrixDestroy)(cudssMatrix_t mat); cudssStatus_t (*cudssDataDestroy)(cudssHandle_t handle, cudssData_t data); cudssStatus_t (*cudssConfigDestroy)(cudssConfig_t config); diff --git a/batch_cudss_benchmark_findings.md b/batch_cudss_benchmark_findings.md new file mode 100644 index 00000000..2132ea92 --- /dev/null +++ b/batch_cudss_benchmark_findings.md @@ -0,0 +1,163 @@ +# Batch cuDSS Benchmark Findings + +This note summarizes the batch GPU backend timing checks performed while adding +cuDSS batch support. + +## Instrumentation + +The CUDA/cuDSS backend was instrumented with opt-in counters around only the +cuDSS factorization and solve `cudssExecute` calls. The timers synchronize the +device immediately before and after each measured cuDSS call. + +Measured phases: + +- Serial solver: `CUDSS_PHASE_FACTORIZATION` and `CUDSS_PHASE_SOLVE` calls from + the normal per-solver cuDSS backend. +- Batch solver: `CUDSS_PHASE_FACTORIZATION` and `CUDSS_PHASE_SOLVE` calls from + the cuDSS batch matrices created with `cudssMatrixCreateBatchCsr` and + `cudssMatrixCreateBatchDn`. + +Setup and cuDSS analysis are excluded from the cuDSS factor/solve totals. + +## Small SOCP Benchmark + +Problem: `cvxpy_qoco.socp_0` + +- `n = 3` +- `m = 3` +- `p = 2` +- one SOC of size `3` +- 100 variants with small deterministic perturbations to `b` and `h` +- CUDA backend +- `ruiz_iters = 0` +- `max_ir_iters = 0` + +Results: + +```text +serial setup: 0.008239 s +serial update: 0.001327 s +serial solve: 3.508448 s +serial update+solve: 3.509775 s +serial solved: 100 / 100 + +serial cuDSS factor: 0.019527 s (594 calls) +serial cuDSS solve: 0.038367 s (1088 calls) +serial cuDSS factor+solve: 0.057893 s +``` + +```text +batch setup: 0.757654 s +batch update: 0.001081 s +batch solve: 3.519075 s +batch update+solve: 3.520156 s +batch solved: 100 / 100 + +batch cuDSS factor: 0.000474 s (6 calls) +batch cuDSS solve: 0.000840 s (11 calls) +batch cuDSS factor+solve: 0.001314 s +``` + +Speedups: + +```text +solve speedup serial/batch: 0.997x +update+solve speedup serial/batch: 0.997x +including setup speedup serial/batch: 0.822x + +cuDSS solve speedup serial/batch: 45.677x +cuDSS factor+solve speedup: 44.060x +``` + +Interpretation: for this tiny SOCP, the batched cuDSS calls are much faster, but +the full solver runtime is dominated by per-item QOCO work outside cuDSS. + +## PDG Benchmark + +Problem: PDG test problem with 100 variants of the initial condition entries in +`b`. + +- `n = 2698` +- `nsoc = 598` +- initial condition offset in `b`: `1794` +- base initial condition: `100 50 50 -9 5 -9` +- relative perturbation: `1e-3` +- batch width: `100` +- CUDA backend + +Results: + +```text +serial setup: 0.031077 s +serial update: 0.009657 s +serial solve: 9.073872 s +serial update+solve: 9.083528 s +serial solved: 100 / 100 + +serial cuDSS factor: 0.325587 s (1245 calls) +serial cuDSS solve: 0.536786 s (3690 calls) +serial cuDSS factor+solve: 0.862373 s +``` + +```text +batch setup: 4.563072 s +batch update: 0.010039 s +batch solve: 8.693753 s +batch update+solve: 8.703792 s +batch solved: 100 / 100 + +batch cuDSS factor: 0.146206 s (13 calls) +batch cuDSS solve: 0.182242 s (38 calls) +batch cuDSS factor+solve: 0.328448 s +``` + +Speedups: + +```text +solve speedup serial/batch: 1.0437x +update+solve speedup serial/batch: 1.0436x +including setup speedup serial/batch: 0.6870x + +cuDSS solve speedup serial/batch: 2.9455x +cuDSS factor+solve speedup: 2.6256x +``` + +Interpretation: for PDG, batched cuDSS is clearly faster internally, but the +end-to-end solve only improves by about 4.4% excluding setup. Non-cuDSS QOCO +work still dominates total runtime. + +## Why Batch Setup Is Expensive + +The current `qoco_batch_setup` implementation allocates and fully initializes +one complete `QOCOSolver` per batch item: + +```text +for each item: + allocate QOCOSolver + call qoco_setup(...) +``` + +For the PDG benchmark with 100 items, batch setup does roughly: + +```text +100x normal solver setup ++ batch CSR pointer arrays and dense wrappers ++ batch cuDSS analysis +``` + +Each per-item `qoco_setup` builds/scales problem data, constructs KKT structures, +creates individual cuDSS handles/matrices/data, and runs per-solver cuDSS +analysis. The serial benchmark only sets up one solver and then updates `b` for +each instance, so its setup time is much lower. + +## Main Takeaways + +- The cuDSS batch API is working and reducing total cuDSS factor/solve time. +- End-to-end speedup is currently limited by QOCO-side per-item work around the + batched linear solves. +- Batch setup is expensive because the batch API currently owns many fully + initialized solver instances rather than sharing immutable problem data and + symbolic structure. +- Improving setup and end-to-end solve time likely requires a more shared batch + representation: one common problem/KKT structure, per-item vector/workspace + state, and per-item numeric KKT values only where needed. diff --git a/benchmarks/pdg_batch_compare.cpp b/benchmarks/pdg_batch_compare.cpp new file mode 100644 index 00000000..c861207a --- /dev/null +++ b/benchmarks/pdg_batch_compare.cpp @@ -0,0 +1,341 @@ +#include "pdg_data.h" +#include "qoco.h" + +#include +#include +#include +#include +#include + +namespace { + +constexpr QOCOInt kDefaultNumInstances = 1000; +constexpr QOCOInt kDefaultBatchWidth = 1000; +constexpr QOCOInt kPdgInitialConditionStart = 6 * (300 - 1); +constexpr QOCOFloat kRelativePerturbation = 1e-3; + +typedef struct { + double serial_factor_sec; + double serial_solve_sec; + long long serial_factor_calls; + long long serial_solve_calls; + double batch_factor_sec; + double batch_solve_sec; + long long batch_factor_calls; + long long batch_solve_calls; +} QOCOCudaLinsysTiming; + +extern "C" void qoco_cuda_linsys_timing_reset(void); +extern "C" void qoco_cuda_linsys_timing_set_enabled(int enabled); +extern "C" void qoco_cuda_linsys_timing_get(QOCOCudaLinsysTiming* timing); + +double now_seconds() +{ + using clock = std::chrono::steady_clock; + return std::chrono::duration(clock::now().time_since_epoch()).count(); +} + +void set_pdg_matrices(QOCOCscMatrix* P, QOCOCscMatrix* A, QOCOCscMatrix* G) +{ + qoco_set_csc(P, pdg_n, pdg_n, pdg_P_nnz, pdg_P_x, pdg_P_p, pdg_P_i); + qoco_set_csc(A, pdg_p, pdg_n, pdg_A_nnz, pdg_A_x, pdg_A_p, pdg_A_i); + qoco_set_csc(G, pdg_m, pdg_n, pdg_G_nnz, pdg_G_x, pdg_G_p, pdg_G_i); +} + +void make_b_values(std::vector& b_values, QOCOInt num_instances) +{ + b_values.resize((size_t)num_instances * (size_t)pdg_p); + + for (QOCOInt item = 0; item < num_instances; ++item) { + QOCOFloat* b = b_values.data() + (size_t)item * (size_t)pdg_p; + for (QOCOInt i = 0; i < pdg_p; ++i) { + b[i] = pdg_b[i]; + } + + for (QOCOInt i = 0; i < 6; ++i) { + const QOCOFloat base = pdg_b[kPdgInitialConditionStart + i]; + const QOCOFloat perturb = + kRelativePerturbation * base * + std::sin(0.017 * (double)(item + 1) * (double)(i + 1)); + b[kPdgInitialConditionStart + i] = base + perturb; + } + } +} + +const QOCOFloat* b_item(const std::vector& b_values, QOCOInt item) +{ + return b_values.data() + (size_t)item * (size_t)pdg_p; +} + +bool status_ok(QOCOInt status) +{ + return status == QOCO_SOLVED || status == QOCO_SOLVED_INACCURATE; +} + +void print_cudss_timing(const char* label, double factor_sec, + double solve_sec, long long factor_calls, + long long solve_calls) +{ + std::printf( + "%s cuDSS factor seconds: %.6f (%lld calls)\n" + "%s cuDSS solve seconds: %.6f (%lld calls)\n" + "%s cuDSS factor+solve seconds: %.6f\n", + label, factor_sec, factor_calls, label, solve_sec, solve_calls, label, + factor_sec + solve_sec); +} + +void warmup(QOCOCscMatrix* P, QOCOCscMatrix* A, QOCOCscMatrix* G, + QOCOSettings* settings, const std::vector& b_values) +{ + QOCOSolver* solver = (QOCOSolver*)malloc(sizeof(QOCOSolver)); + if (!solver) { + std::fprintf(stderr, "warmup solver allocation failed\n"); + std::exit(1); + } + + QOCOInt exit = qoco_setup(solver, pdg_n, pdg_m, pdg_p, P, pdg_c, A, pdg_b, G, + pdg_h, pdg_l, pdg_nsoc, pdg_q, settings); + if (exit != QOCO_NO_ERROR) { + std::fprintf(stderr, "warmup setup failed: %d\n", exit); + std::exit(1); + } + + qoco_update_vector_data(solver, nullptr, (QOCOFloat*)b_item(b_values, 0), + nullptr); + exit = qoco_solve(solver); + if (!status_ok(exit)) { + std::fprintf(stderr, "warmup solve failed: %d\n", exit); + std::exit(1); + } + qoco_cleanup(solver); +} + +} // namespace + +int main(int argc, char** argv) +{ + std::setvbuf(stdout, nullptr, _IOLBF, 0); + + QOCOInt num_instances = kDefaultNumInstances; + QOCOInt batch_width = kDefaultBatchWidth; + bool run_serial = true; + bool run_batch = true; + + if (argc > 1) { + num_instances = (QOCOInt)std::atoi(argv[1]); + } + if (argc > 2) { + batch_width = (QOCOInt)std::atoi(argv[2]); + } + if (argc > 3) { + run_serial = std::atoi(argv[3]) != 0; + } + if (argc > 4) { + run_batch = std::atoi(argv[4]) != 0; + } + if (num_instances <= 0 || batch_width <= 0) { + std::fprintf(stderr, + "usage: %s [num_instances=1000] [batch_width=1000] " + "[run_serial=1] [run_batch=1]\n", + argv[0]); + return 1; + } + + QOCOCscMatrix P; + QOCOCscMatrix A; + QOCOCscMatrix G; + set_pdg_matrices(&P, &A, &G); + + QOCOSettings settings; + set_default_settings(&settings); + settings.verbose = 0; + + std::vector b_values; + make_b_values(b_values, num_instances); + + std::printf("PDG batch comparison\n"); + std::printf("instances: %d\n", (int)num_instances); + std::printf("batch width: %d\n", (int)batch_width); + std::printf("initial condition b offset: %d\n", + (int)kPdgInitialConditionStart); + std::printf("relative perturbation: %.3e\n", + (double)kRelativePerturbation); + std::printf("base initial condition:"); + for (QOCOInt i = 0; i < 6; ++i) { + std::printf(" %.9g", (double)pdg_b[kPdgInitialConditionStart + i]); + } + std::printf("\n"); + + warmup(&P, &A, &G, &settings, b_values); + + double serial_setup_sec = 0.0; + double serial_update_sec = 0.0; + double serial_solve_sec = 0.0; + QOCOInt serial_solved = 0; + QOCOInt serial_failed = 0; + double serial_obj_accum = 0.0; + QOCOCudaLinsysTiming serial_cudss = {}; + + if (run_serial) { + QOCOSolver* serial_solver = (QOCOSolver*)malloc(sizeof(QOCOSolver)); + if (!serial_solver) { + std::fprintf(stderr, "serial solver allocation failed\n"); + return 1; + } + + double t0 = now_seconds(); + QOCOInt exit = qoco_setup(serial_solver, pdg_n, pdg_m, pdg_p, &P, pdg_c, + &A, pdg_b, &G, pdg_h, pdg_l, pdg_nsoc, pdg_q, + &settings); + serial_setup_sec = now_seconds() - t0; + if (exit != QOCO_NO_ERROR) { + std::fprintf(stderr, "serial setup failed: %d\n", exit); + return 1; + } + + qoco_cuda_linsys_timing_reset(); + for (QOCOInt item = 0; item < num_instances; ++item) { + t0 = now_seconds(); + qoco_update_vector_data(serial_solver, nullptr, + (QOCOFloat*)b_item(b_values, item), nullptr); + serial_update_sec += now_seconds() - t0; + + t0 = now_seconds(); + exit = qoco_solve(serial_solver); + serial_solve_sec += now_seconds() - t0; + + if (status_ok(exit)) { + ++serial_solved; + } + else { + ++serial_failed; + } + serial_obj_accum += serial_solver->sol->obj; + } + qoco_cuda_linsys_timing_get(&serial_cudss); + qoco_cleanup(serial_solver); + } + + const double serial_total_sec = serial_update_sec + serial_solve_sec; + + if (run_serial) { + std::printf("\nserial setup seconds: %.6f\n", serial_setup_sec); + std::printf("serial update seconds: %.6f\n", serial_update_sec); + std::printf("serial solve seconds: %.6f\n", serial_solve_sec); + std::printf("serial update+solve seconds: %.6f\n", serial_total_sec); + std::printf("serial solved: %d failed: %d obj checksum: %.12e\n", + (int)serial_solved, (int)serial_failed, serial_obj_accum); + print_cudss_timing("serial", serial_cudss.serial_factor_sec, + serial_cudss.serial_solve_sec, + serial_cudss.serial_factor_calls, + serial_cudss.serial_solve_calls); + } + else { + std::printf("\nserial skipped\n"); + } + + if (!run_batch) { + return serial_failed; + } + + double batch_setup_sec = 0.0; + double batch_update_sec = 0.0; + double batch_solve_sec = 0.0; + QOCOInt batch_solved = 0; + QOCOInt batch_failed = 0; + double batch_obj_accum = 0.0; + QOCOCudaLinsysTiming batch_cudss = {}; + + qoco_cuda_linsys_timing_reset(); + for (QOCOInt offset = 0; offset < num_instances; offset += batch_width) { + const QOCOInt this_batch = + (offset + batch_width <= num_instances) ? batch_width + : (num_instances - offset); + std::printf("batch chunk offset %d size %d\n", (int)offset, + (int)this_batch); + + QOCOBatchSolver batch; + double t0 = now_seconds(); + QOCOInt exit = qoco_batch_setup(&batch, this_batch, pdg_n, pdg_m, pdg_p, + &P, pdg_c, &A, pdg_b, &G, pdg_h, pdg_l, + pdg_nsoc, pdg_q, &settings); + batch_setup_sec += now_seconds() - t0; + if (exit != QOCO_NO_ERROR) { + std::fprintf(stderr, "batch setup failed at offset %d: %d\n", + (int)offset, exit); + return 1; + } + + t0 = now_seconds(); + for (QOCOInt item = 0; item < this_batch; ++item) { + exit = qoco_batch_update_vector_data( + &batch, item, nullptr, + (QOCOFloat*)b_item(b_values, offset + item), nullptr); + if (exit != QOCO_NO_ERROR) { + std::fprintf(stderr, "batch update failed for item %d: %d\n", + (int)(offset + item), exit); + qoco_batch_cleanup(&batch); + return 1; + } + } + batch_update_sec += now_seconds() - t0; + + t0 = now_seconds(); + exit = qoco_batch_solve(&batch); + batch_solve_sec += now_seconds() - t0; + if (exit != QOCO_NO_ERROR) { + std::fprintf(stderr, "batch solve dispatch failed at offset %d: %d\n", + (int)offset, exit); + qoco_batch_cleanup(&batch); + return 1; + } + + for (QOCOInt item = 0; item < this_batch; ++item) { + if (status_ok(batch.statuses[item])) { + ++batch_solved; + } + else { + ++batch_failed; + } + QOCOSolution* sol = qoco_batch_get_solution(&batch, item); + if (sol) { + batch_obj_accum += sol->obj; + } + } + qoco_batch_cleanup(&batch); + } + qoco_cuda_linsys_timing_get(&batch_cudss); + qoco_cuda_linsys_timing_set_enabled(0); + + const double batch_total_sec = batch_update_sec + batch_solve_sec; + + std::printf("\nbatch setup seconds: %.6f\n", batch_setup_sec); + std::printf("batch update seconds: %.6f\n", batch_update_sec); + std::printf("batch solve seconds: %.6f\n", batch_solve_sec); + std::printf("batch update+solve seconds: %.6f\n", batch_total_sec); + std::printf("batch solved: %d failed: %d obj checksum: %.12e\n", + (int)batch_solved, (int)batch_failed, batch_obj_accum); + print_cudss_timing("batch", batch_cudss.batch_factor_sec, + batch_cudss.batch_solve_sec, + batch_cudss.batch_factor_calls, + batch_cudss.batch_solve_calls); + + if (run_serial) { + std::printf("\nsolve speedup serial/batch: %.6f\n", + serial_solve_sec / batch_solve_sec); + std::printf("update+solve speedup serial/batch: %.6f\n", + serial_total_sec / batch_total_sec); + std::printf("including setup speedup serial/batch: %.6f\n", + (serial_setup_sec + serial_total_sec) / + (batch_setup_sec + batch_total_sec)); + std::printf("cuDSS solve speedup serial/batch: %.6f\n", + serial_cudss.serial_solve_sec / batch_cudss.batch_solve_sec); + std::printf("cuDSS factor+solve speedup serial/batch: %.6f\n", + (serial_cudss.serial_factor_sec + + serial_cudss.serial_solve_sec) / + (batch_cudss.batch_factor_sec + + batch_cudss.batch_solve_sec)); + } + + return batch_failed || serial_failed; +} diff --git a/include/qoco.h b/include/qoco.h index ea7f1c94..1953844f 100644 --- a/include/qoco.h +++ b/include/qoco.h @@ -22,6 +22,7 @@ extern "C" { #include "cone.h" #include "definitions.h" #include "qoco_api.h" +#include "qoco_batch.h" #include "qoco_utils.h" #include "structs.h" @@ -29,4 +30,4 @@ extern "C" { } #endif -#endif /* #ifndef QOCO_H */ \ No newline at end of file +#endif /* #ifndef QOCO_H */ diff --git a/include/qoco_batch.h b/include/qoco_batch.h new file mode 100644 index 00000000..1c9e5ed5 --- /dev/null +++ b/include/qoco_batch.h @@ -0,0 +1,90 @@ +/** + * @file qoco_batch.h + * @author Govind M. Chari + * + * @section LICENSE + * + * Copyright (c) 2026, Govind M. Chari + * This source code is licensed under the BSD 3-Clause License + * + * @section DESCRIPTION + * + * Exposes the batched QOCO API. + */ + +#ifndef QOCO_BATCH_H +#define QOCO_BATCH_H + +#include "definitions.h" +#include "qoco_api.h" +#include "structs.h" + +#ifdef __cplusplus +extern "C" { +#endif + +/** + * @brief Allocates a batch of QOCO solvers with common dimensions and sparsity. + * + * @param batch Pointer to batch solver. + * @param batch_count Number of solver instances. + * @param n Number of optimization variables. + * @param m Number of conic constraints. + * @param p Number of affine equality constraints. + * @param P Upper triangular part of quadratic cost Hessian in CSC form. + * @param c Linear cost vector. + * @param A Affine equality constraint matrix in CSC form. + * @param b Affine equality constraint offset vector. + * @param G Conic constraint matrix in CSC form. + * @param h Conic constraint offset vector. + * @param l Dimension of non-negative orthant. + * @param nsoc Number of second-order cones. + * @param q Dimension of each second-order cone. + * @param settings Settings struct. + * @return 0 if no error or flag containing error code. + */ +QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, + QOCOInt n, QOCOInt m, QOCOInt p, QOCOCscMatrix* P, + QOCOFloat* c, QOCOCscMatrix* A, QOCOFloat* b, + QOCOCscMatrix* G, QOCOFloat* h, QOCOInt l, + QOCOInt nsoc, QOCOInt* q, QOCOSettings* settings); + +/** + * @brief Updates vector data for one item in a batch. + */ +QOCOInt qoco_batch_update_vector_data(QOCOBatchSolver* batch, QOCOInt item, + QOCOFloat* cnew, QOCOFloat* bnew, + QOCOFloat* hnew); + +/** + * @brief Updates matrix values for one item in a batch. Sparsity must match + * the matrices supplied to qoco_batch_setup. + */ +QOCOInt qoco_batch_update_matrix_data(QOCOBatchSolver* batch, QOCOInt item, + QOCOFloat* Pxnew, QOCOFloat* Axnew, + QOCOFloat* Gxnew); + +/** + * @brief Solves all items in a batch. + * + * @return QOCO_NO_ERROR if dispatch succeeded. Per-item solve statuses are + * stored in batch->statuses. + */ +QOCOInt qoco_batch_solve(QOCOBatchSolver* batch); + +/** + * @brief Returns the solution for one item in a batch, or NULL for an invalid + * item. + */ +QOCOSolution* qoco_batch_get_solution(QOCOBatchSolver* batch, QOCOInt item); + +/** + * @brief Frees all memory owned by a batch solver. + */ +QOCOInt qoco_batch_cleanup(QOCOBatchSolver* batch); + +#ifdef __cplusplus +} +#endif + +#endif /* #ifndef QOCO_BATCH_H */ diff --git a/include/structs.h b/include/structs.h index 1f7a367f..89cd24ca 100644 --- a/include/structs.h +++ b/include/structs.h @@ -395,4 +395,25 @@ typedef struct { } QOCOSolver; -#endif /* #ifndef QOCO_STRUCTS_H */ \ No newline at end of file +/** + * @brief Batch solver struct. Owns a set of QOCO solver instances with common + * dimensions and sparsity structure. + */ +typedef struct { + /** Number of solver instances in the batch. */ + QOCOInt batch_count; + + /** Solver instances owned by the batch solver. */ + QOCOSolver** solvers; + + /** Solve status for each batch item. */ + QOCOInt* statuses; + + /** Backend-specific batch linear-system data. */ + void* batch_linsys_data; + + /** Whether backend-specific batch data needs to be rebuilt. */ + unsigned char batch_linsys_stale; +} QOCOBatchSolver; + +#endif /* #ifndef QOCO_STRUCTS_H */ diff --git a/src/qoco_batch.c b/src/qoco_batch.c new file mode 100644 index 00000000..f1fcbb89 --- /dev/null +++ b/src/qoco_batch.c @@ -0,0 +1,164 @@ +/** + * @file qoco_batch.c + * @author Govind M. Chari + * + * @section LICENSE + * + * Copyright (c) 2026, Govind M. Chari + * This source code is licensed under the BSD 3-Clause License + * + * @section DESCRIPTION + * + * Implements the batched QOCO API. + */ + +#include "qoco_batch.h" + +#ifdef QOCO_ALGEBRA_BACKEND_CUDA +QOCOInt qoco_cuda_batch_setup(QOCOBatchSolver* batch); +QOCOInt qoco_cuda_batch_solve(QOCOBatchSolver* batch); +void qoco_cuda_batch_cleanup(QOCOBatchSolver* batch); +#endif + +static unsigned char qoco_batch_valid_item(const QOCOBatchSolver* batch, + QOCOInt item) +{ + return batch && item >= 0 && item < batch->batch_count && batch->solvers && + batch->solvers[item]; +} + +static void qoco_batch_zero(QOCOBatchSolver* batch) +{ + batch->batch_count = 0; + batch->solvers = NULL; + batch->statuses = NULL; + batch->batch_linsys_data = NULL; + batch->batch_linsys_stale = 0; +} + +QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, + QOCOInt n, QOCOInt m, QOCOInt p, QOCOCscMatrix* P, + QOCOFloat* c, QOCOCscMatrix* A, QOCOFloat* b, + QOCOCscMatrix* G, QOCOFloat* h, QOCOInt l, + QOCOInt nsoc, QOCOInt* q, QOCOSettings* settings) +{ + if (!batch || batch_count <= 0) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + + qoco_batch_zero(batch); + batch->batch_count = batch_count; + batch->solvers = + (QOCOSolver**)qoco_calloc(batch_count, sizeof(QOCOSolver*)); + batch->statuses = (QOCOInt*)qoco_calloc(batch_count, sizeof(QOCOInt)); + if (!batch->solvers || !batch->statuses) { + qoco_batch_cleanup(batch); + return qoco_error(QOCO_MALLOC_ERROR); + } + + for (QOCOInt item = 0; item < batch_count; ++item) { + batch->solvers[item] = (QOCOSolver*)qoco_malloc(sizeof(QOCOSolver)); + if (!batch->solvers[item]) { + qoco_batch_cleanup(batch); + return qoco_error(QOCO_MALLOC_ERROR); + } + + QOCOInt exit = + qoco_setup(batch->solvers[item], n, m, p, P, c, A, b, G, h, l, nsoc, + q, settings); + if (exit != QOCO_NO_ERROR) { + qoco_free(batch->solvers[item]); + batch->solvers[item] = NULL; + qoco_batch_cleanup(batch); + return exit; + } + batch->statuses[item] = QOCO_UNSOLVED; + } + +#ifdef QOCO_ALGEBRA_BACKEND_CUDA + QOCOInt exit = qoco_cuda_batch_setup(batch); + if (exit != QOCO_NO_ERROR) { + qoco_batch_cleanup(batch); + return exit; + } +#endif + + return QOCO_NO_ERROR; +} + +QOCOInt qoco_batch_update_vector_data(QOCOBatchSolver* batch, QOCOInt item, + QOCOFloat* cnew, QOCOFloat* bnew, + QOCOFloat* hnew) +{ + if (!qoco_batch_valid_item(batch, item)) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + + qoco_update_vector_data(batch->solvers[item], cnew, bnew, hnew); + batch->statuses[item] = QOCO_UNSOLVED; + return QOCO_NO_ERROR; +} + +QOCOInt qoco_batch_update_matrix_data(QOCOBatchSolver* batch, QOCOInt item, + QOCOFloat* Pxnew, QOCOFloat* Axnew, + QOCOFloat* Gxnew) +{ + if (!qoco_batch_valid_item(batch, item)) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + + qoco_update_matrix_data(batch->solvers[item], Pxnew, Axnew, Gxnew); + batch->statuses[item] = QOCO_UNSOLVED; + return QOCO_NO_ERROR; +} + +QOCOInt qoco_batch_solve(QOCOBatchSolver* batch) +{ + if (!batch || !batch->solvers || !batch->statuses) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + +#ifdef QOCO_ALGEBRA_BACKEND_CUDA + return qoco_cuda_batch_solve(batch); +#else + for (QOCOInt item = 0; item < batch->batch_count; ++item) { + if (!batch->solvers[item]) { + return qoco_error(QOCO_DATA_VALIDATION_ERROR); + } + batch->statuses[item] = qoco_solve(batch->solvers[item]); + } + return QOCO_NO_ERROR; +#endif +} + +QOCOSolution* qoco_batch_get_solution(QOCOBatchSolver* batch, QOCOInt item) +{ + if (!qoco_batch_valid_item(batch, item)) { + return NULL; + } + return batch->solvers[item]->sol; +} + +QOCOInt qoco_batch_cleanup(QOCOBatchSolver* batch) +{ + if (!batch) { + return QOCO_NO_ERROR; + } + +#ifdef QOCO_ALGEBRA_BACKEND_CUDA + qoco_cuda_batch_cleanup(batch); +#endif + + if (batch->solvers) { + for (QOCOInt item = 0; item < batch->batch_count; ++item) { + if (batch->solvers[item]) { + qoco_cleanup(batch->solvers[item]); + } + } + } + + qoco_free(batch->solvers); + qoco_free(batch->statuses); + qoco_batch_zero(batch); + return QOCO_NO_ERROR; +} diff --git a/tests/unit_tests/batch_test.cpp b/tests/unit_tests/batch_test.cpp new file mode 100644 index 00000000..dbdaf60c --- /dev/null +++ b/tests/unit_tests/batch_test.cpp @@ -0,0 +1,128 @@ +#include "qoco.h" +#include "test_utils.h" +#include "gtest/gtest.h" + +#include + +namespace { + +void set_batch_test_settings(QOCOSettings* settings) +{ + set_default_settings(settings); + settings->ruiz_iters = 0; + settings->verbose = 0; +} + +void set_diag_qp_matrix(QOCOCscMatrix* P, QOCOFloat* Px) +{ + static QOCOInt Pp[] = {0, 1, 2}; + static QOCOInt Pi[] = {0, 1}; + qoco_set_csc(P, 2, 2, 2, Px, Pp, Pi); +} + +QOCOSolver* setup_serial(QOCOFloat* Px, QOCOFloat* c, + QOCOSettings* settings) +{ + QOCOCscMatrix P; + set_diag_qp_matrix(&P, Px); + + QOCOSolver* solver = (QOCOSolver*)std::malloc(sizeof(QOCOSolver)); + QOCOInt exit = qoco_setup(solver, 2, 0, 0, &P, c, nullptr, nullptr, nullptr, + nullptr, 0, 0, nullptr, settings); + EXPECT_EQ(exit, QOCO_NO_ERROR); + if (exit != QOCO_NO_ERROR) { + std::free(solver); + return nullptr; + } + return solver; +} + +} // namespace + +TEST(batch, vector_and_matrix_updates_match_serial_solves) +{ + QOCOFloat Px_base[] = {1.0, 2.0}; + QOCOFloat c_base[] = {1.0, 4.0}; + QOCOFloat c_item0[] = {2.0, -6.0}; + QOCOFloat Px_item1[] = {2.0, 4.0}; + + QOCOSettings settings; + set_batch_test_settings(&settings); + + QOCOCscMatrix P; + set_diag_qp_matrix(&P, Px_base); + + QOCOBatchSolver batch; + QOCOInt exit = qoco_batch_setup(&batch, 2, 2, 0, 0, &P, c_base, nullptr, + nullptr, nullptr, nullptr, 0, 0, nullptr, + &settings); + ASSERT_EQ(exit, QOCO_NO_ERROR); + + EXPECT_EQ(qoco_batch_update_vector_data(&batch, 0, c_item0, nullptr, nullptr), + QOCO_NO_ERROR); + EXPECT_EQ(qoco_batch_update_matrix_data(&batch, 1, Px_item1, nullptr, + nullptr), + QOCO_NO_ERROR); + + QOCOSolver* serial0 = setup_serial(Px_base, c_base, &settings); + ASSERT_NE(serial0, nullptr); + qoco_update_vector_data(serial0, c_item0, nullptr, nullptr); + QOCOInt serial0_status = qoco_solve(serial0); + + QOCOSolver* serial1 = setup_serial(Px_base, c_base, &settings); + ASSERT_NE(serial1, nullptr); + qoco_update_matrix_data(serial1, Px_item1, nullptr, nullptr); + QOCOInt serial1_status = qoco_solve(serial1); + + exit = qoco_batch_solve(&batch); + ASSERT_EQ(exit, QOCO_NO_ERROR); + ASSERT_EQ(batch.statuses[0], serial0_status); + ASSERT_EQ(batch.statuses[1], serial1_status); + + QOCOSolution* sol0 = qoco_batch_get_solution(&batch, 0); + QOCOSolution* sol1 = qoco_batch_get_solution(&batch, 1); + ASSERT_NE(sol0, nullptr); + ASSERT_NE(sol1, nullptr); + + expect_eq_vectorf(sol0->x, serial0->sol->x, 2, 1e-6); + expect_eq_vectorf(sol1->x, serial1->sol->x, 2, 1e-6); + EXPECT_NEAR((double)sol0->obj, (double)serial0->sol->obj, 1e-6); + EXPECT_NEAR((double)sol1->obj, (double)serial1->sol->obj, 1e-6); + + qoco_cleanup(serial0); + qoco_cleanup(serial1); + qoco_batch_cleanup(&batch); +} + +TEST(batch, batch_count_one_matches_serial_solve) +{ + QOCOFloat Px[] = {1.0, 2.0}; + QOCOFloat c[] = {-3.0, 4.0}; + + QOCOSettings settings; + set_batch_test_settings(&settings); + + QOCOCscMatrix P; + set_diag_qp_matrix(&P, Px); + + QOCOBatchSolver batch; + QOCOInt exit = qoco_batch_setup(&batch, 1, 2, 0, 0, &P, c, nullptr, nullptr, + nullptr, nullptr, 0, 0, nullptr, &settings); + ASSERT_EQ(exit, QOCO_NO_ERROR); + + QOCOSolver* serial = setup_serial(Px, c, &settings); + ASSERT_NE(serial, nullptr); + QOCOInt serial_status = qoco_solve(serial); + + exit = qoco_batch_solve(&batch); + ASSERT_EQ(exit, QOCO_NO_ERROR); + ASSERT_EQ(batch.statuses[0], serial_status); + + QOCOSolution* sol = qoco_batch_get_solution(&batch, 0); + ASSERT_NE(sol, nullptr); + expect_eq_vectorf(sol->x, serial->sol->x, 2, 1e-6); + EXPECT_EQ(qoco_batch_get_solution(&batch, 1), nullptr); + + qoco_cleanup(serial); + qoco_batch_cleanup(&batch); +} From 1da1bbea9c3c3d124276f58b37dc22199a08d00f Mon Sep 17 00:00:00 2001 From: govindchari Date: Tue, 26 May 2026 11:44:59 -0700 Subject: [PATCH 2/2] Clang format --- include/qoco_batch.h | 10 +++++----- src/qoco_batch.c | 18 ++++++++---------- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/include/qoco_batch.h b/include/qoco_batch.h index 1c9e5ed5..071294dd 100644 --- a/include/qoco_batch.h +++ b/include/qoco_batch.h @@ -43,11 +43,11 @@ extern "C" { * @param settings Settings struct. * @return 0 if no error or flag containing error code. */ -QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, - QOCOInt n, QOCOInt m, QOCOInt p, QOCOCscMatrix* P, - QOCOFloat* c, QOCOCscMatrix* A, QOCOFloat* b, - QOCOCscMatrix* G, QOCOFloat* h, QOCOInt l, - QOCOInt nsoc, QOCOInt* q, QOCOSettings* settings); +QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, QOCOInt n, + QOCOInt m, QOCOInt p, QOCOCscMatrix* P, QOCOFloat* c, + QOCOCscMatrix* A, QOCOFloat* b, QOCOCscMatrix* G, + QOCOFloat* h, QOCOInt l, QOCOInt nsoc, QOCOInt* q, + QOCOSettings* settings); /** * @brief Updates vector data for one item in a batch. diff --git a/src/qoco_batch.c b/src/qoco_batch.c index f1fcbb89..91ffa870 100644 --- a/src/qoco_batch.c +++ b/src/qoco_batch.c @@ -36,11 +36,11 @@ static void qoco_batch_zero(QOCOBatchSolver* batch) batch->batch_linsys_stale = 0; } -QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, - QOCOInt n, QOCOInt m, QOCOInt p, QOCOCscMatrix* P, - QOCOFloat* c, QOCOCscMatrix* A, QOCOFloat* b, - QOCOCscMatrix* G, QOCOFloat* h, QOCOInt l, - QOCOInt nsoc, QOCOInt* q, QOCOSettings* settings) +QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, QOCOInt n, + QOCOInt m, QOCOInt p, QOCOCscMatrix* P, QOCOFloat* c, + QOCOCscMatrix* A, QOCOFloat* b, QOCOCscMatrix* G, + QOCOFloat* h, QOCOInt l, QOCOInt nsoc, QOCOInt* q, + QOCOSettings* settings) { if (!batch || batch_count <= 0) { return qoco_error(QOCO_DATA_VALIDATION_ERROR); @@ -48,8 +48,7 @@ QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, qoco_batch_zero(batch); batch->batch_count = batch_count; - batch->solvers = - (QOCOSolver**)qoco_calloc(batch_count, sizeof(QOCOSolver*)); + batch->solvers = (QOCOSolver**)qoco_calloc(batch_count, sizeof(QOCOSolver*)); batch->statuses = (QOCOInt*)qoco_calloc(batch_count, sizeof(QOCOInt)); if (!batch->solvers || !batch->statuses) { qoco_batch_cleanup(batch); @@ -63,9 +62,8 @@ QOCOInt qoco_batch_setup(QOCOBatchSolver* batch, QOCOInt batch_count, return qoco_error(QOCO_MALLOC_ERROR); } - QOCOInt exit = - qoco_setup(batch->solvers[item], n, m, p, P, c, A, b, G, h, l, nsoc, - q, settings); + QOCOInt exit = qoco_setup(batch->solvers[item], n, m, p, P, c, A, b, G, h, + l, nsoc, q, settings); if (exit != QOCO_NO_ERROR) { qoco_free(batch->solvers[item]); batch->solvers[item] = NULL;