From 910b08f46757e5bfdab70d16efbb326e14306bdc Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sun, 4 Oct 2026 15:22:09 -0400 Subject: [PATCH 1/3] Build inequality Jacobian structure without PETSc. --- .../model/power_bal_hiop/paramsrajahiop.cpp | 27 ++-- .../model/power_bal_hiop/paramsrajahiop.h | 4 + .../pbpolrajahiopsparse_gpu.cpp | 45 +++--- .../pbpolrajahiopsparse_gpu.hpp | 9 +- .../pbpolrajahiopsparsekernels.cpp | 145 ++++++++++++++---- tests/unit/opflow/test_ineqjac_gpu.cpp | 26 +++- 6 files changed, 184 insertions(+), 72 deletions(-) diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp index 4ff8581f..d66787e6 100644 --- a/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp +++ b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp @@ -914,29 +914,32 @@ void PbpolModelRajaHiop::destroy(OPFLOW opflow) { genparams.destroy(opflow); #ifdef EXAGO_ENABLE_GPU - if (i_jaceq != NULL) { - // These arrays get allocated only with sparse GPU model. For other models, - // they are not allocated. Hence, this business of checking if the array is - // NULL - - auto &resmgr = umpire::ResourceManager::getInstance(); - umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); - umpire::Allocator d_allocator_ = resmgr.getAllocator("DEVICE"); + // These arrays get allocated only with sparse GPU model, each group when its + // sparsity pattern is first computed. Hence, this business of checking if + // the array is NULL + auto &resmgr = umpire::ResourceManager::getInstance(); + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + umpire::Allocator d_allocator_ = resmgr.getAllocator("DEVICE"); + if (i_jaceq != NULL) { h_allocator_.deallocate(i_jaceq); h_allocator_.deallocate(j_jaceq); h_allocator_.deallocate(perm_jaceq); d_allocator_.deallocate(perm_jaceq_dev); + } + if (i_hess != NULL) { h_allocator_.deallocate(i_hess); h_allocator_.deallocate(j_hess); h_allocator_.deallocate(perm_hess); d_allocator_.deallocate(perm_hess_dev); + } - if (opflow->nconineq) { - h_allocator_.deallocate(i_jacineq); - h_allocator_.deallocate(j_jacineq); - } + if (i_jacineq != NULL) { + h_allocator_.deallocate(i_jacineq); + h_allocator_.deallocate(j_jacineq); + h_allocator_.deallocate(perm_jacineq); + d_allocator_.deallocate(perm_jacineq_dev); } #endif } diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.h b/src/opflow/model/power_bal_hiop/paramsrajahiop.h index 087708fb..8f7a1ebb 100644 --- a/src/opflow/model/power_bal_hiop/paramsrajahiop.h +++ b/src/opflow/model/power_bal_hiop/paramsrajahiop.h @@ -282,6 +282,7 @@ struct PbpolModelRajaHiop : public _p_FormPBPOLRAJAHIOP { i_jacineq = j_jacineq = NULL; i_hess = j_hess = NULL; perm_jaceq = perm_jaceq_dev = NULL; + perm_jacineq = perm_jacineq_dev = NULL; perm_hess = perm_hess_dev = NULL; } @@ -307,6 +308,9 @@ struct PbpolModelRajaHiop : public _p_FormPBPOLRAJAHIOP { *j_hess; // Row and column indices for hessian int *perm_jaceq, *perm_jaceq_dev; // Permutation for equality constraints Jacobian indices + int *perm_jacineq, + *perm_jacineq_dev; // Permutation for inequality constraints Jacobian + // indices int *perm_hess, *perm_hess_dev; // Permutation for equality constraints Jacobian indices }; diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.cpp index eaeb04e7..1190af6c 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.cpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.cpp @@ -11,6 +11,7 @@ void ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, const double *x_dev, + const int *perm_dev, double *jacd_dev) { PbpolModelRajaHiop *pbpolrajahiopsparse = reinterpret_cast(opflow->model); @@ -55,14 +56,14 @@ void ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, double apf_delP_minus_delPg = apf * delP - delPg; /* Row 0: d/d{Pg, delPg, delP} of (apf*delP - delPg)*(Pg - pt) */ - jacd_dev[base + 0] = apf_delP_minus_delPg; - jacd_dev[base + 1] = -(Pg - pt); - jacd_dev[base + 2] = apf * (Pg - pt); + jacd_dev[perm_dev[base + 0]] = apf_delP_minus_delPg; + jacd_dev[perm_dev[base + 1]] = -(Pg - pt); + jacd_dev[perm_dev[base + 2]] = apf * (Pg - pt); /* Row 1: d/d{Pg, delPg, delP} of (delPg - apf*delP)*(pb - Pg) */ - jacd_dev[base + 3] = apf_delP_minus_delPg; - jacd_dev[base + 4] = pb - Pg; - jacd_dev[base + 5] = -apf * (pb - Pg); + jacd_dev[perm_dev[base + 3]] = apf_delP_minus_delPg; + jacd_dev[perm_dev[base + 4]] = pb - Pg; + jacd_dev[perm_dev[base + 5]] = -apf * (pb - Pg); }); } @@ -118,15 +119,15 @@ void ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, double dQg = Vset - V; /* Row 0 entry for this gen's Qg */ - jacd_dev[base + k] = dQg; + jacd_dev[perm_dev[base + k]] = dQg; /* Row 1 entry for this gen's Qg */ - jacd_dev[base + ngen + 1 + k] = dQg; + jacd_dev[perm_dev[base + ngen + 1 + k]] = dQg; } /* Row 0: V derivative entry (last in this row) */ - jacd_dev[base + ngen] = Qmax - Q; + jacd_dev[perm_dev[base + ngen]] = Qmax - Q; /* Row 1: V derivative entry (last in this row) */ - jacd_dev[base + ngen + 1 + ngen] = Qmin - Q; + jacd_dev[perm_dev[base + ngen + 1 + ngen]] = Qmin - Q; }); } @@ -216,25 +217,29 @@ void ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, double dQt_dVmf = Vmt * (-Btf * cos_tf + Gtf * sin_tf); /* Row 0 (Sf2): derivatives w.r.t. thetaf, Vmf, thetat, Vmt */ - jacd_dev[base + 0] = dSf2_dPf * dPf_dthetaf + dSf2_dQf * dQf_dthetaf; - jacd_dev[base + 1] = dSf2_dPf * dPf_dVmf + dSf2_dQf * dQf_dVmf; - jacd_dev[base + 2] = dSf2_dPf * dPf_dthetat + dSf2_dQf * dQf_dthetat; - jacd_dev[base + 3] = dSf2_dPf * dPf_dVmt + dSf2_dQf * dQf_dVmt; + jacd_dev[perm_dev[base + 0]] = + dSf2_dPf * dPf_dthetaf + dSf2_dQf * dQf_dthetaf; + jacd_dev[perm_dev[base + 1]] = + dSf2_dPf * dPf_dVmf + dSf2_dQf * dQf_dVmf; + jacd_dev[perm_dev[base + 2]] = + dSf2_dPf * dPf_dthetat + dSf2_dQf * dQf_dthetat; + jacd_dev[perm_dev[base + 3]] = + dSf2_dPf * dPf_dVmt + dSf2_dQf * dQf_dVmt; /* Row 1 (St2): derivatives w.r.t. thetaf, Vmf, thetat, Vmt */ - jacd_dev[base + row_stride + 0] = + jacd_dev[perm_dev[base + row_stride + 0]] = dSt2_dPt * dPt_dthetaf + dSt2_dQt * dQt_dthetaf; - jacd_dev[base + row_stride + 1] = + jacd_dev[perm_dev[base + row_stride + 1]] = dSt2_dPt * dPt_dVmf + dSt2_dQt * dQt_dVmf; - jacd_dev[base + row_stride + 2] = + jacd_dev[perm_dev[base + row_stride + 2]] = dSt2_dPt * dPt_dthetat + dSt2_dQt * dQt_dthetat; - jacd_dev[base + row_stride + 3] = + jacd_dev[perm_dev[base + row_stride + 3]] = dSt2_dPt * dPt_dVmt + dSt2_dQt * dQt_dVmt; /* Slack variable entries (Step 5) */ if (has_slack) { - jacd_dev[base + 4] = -1.0; - jacd_dev[base + row_stride + 4] = -1.0; + jacd_dev[perm_dev[base + 4]] = -1.0; + jacd_dev[perm_dev[base + row_stride + 4]] = -1.0; } }); } diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.hpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.hpp index c946aa0a..fcd31ff0 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.hpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse_gpu.hpp @@ -11,19 +11,20 @@ /** * GPU-only (PETSc-free) computation of inequality constraint Jacobian values. * - * Replaces the PETSc-based path that calls - * opflow->modelops.computeinequalityconstraintjacobian followed by - * MatGetRow extraction. Writes directly into device memory using RAJA - * kernels, without PETSc Mat/Vec operations; no H2D, D2H copies back and forth. + * Writes directly into device memory using RAJA kernels, without PETSc + * Mat/Vec operations; no H2D, D2H copies back and forth. * * @param opflow The OPFLOW problem context * @param x_dev Device array of variable values + * @param perm_dev Device array mapping setup-layout positions to sorted + * (row, column) positions in the inequality block * @param jacd_dev Device output array for inequality Jacobian values * (points to the ineq portion of the sparse Jacobian, * i.e. MJacS_dev + nnz_eqjacsp) */ void ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, const double *x_dev, + const int *perm_dev, double *jacd_dev); /** diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp index 2326a61c..be470d1c 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp @@ -578,13 +578,7 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( PbpolModelRajaHiop *pbpolrajahiopsparse = reinterpret_cast(opflow->model); PetscErrorCode ierr; - PetscInt *iRowstart, *jColstart; PetscInt roffset, coffset; - PetscInt nrow, ncol; - PetscInt nvals; - const PetscInt *cols; - const PetscScalar *vals; - PetscInt i, j; auto &resmgr = umpire::ResourceManager::getInstance(); PetscFunctionBegin; @@ -592,53 +586,135 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( // If sparsity pattern does not exist, create it! if (iJacS_dev != NULL && jJacS_dev != NULL) { - // Create arrays on host to store i,j, and val arrays if (opflow->Nconineq) { + // Compute sparsity pattern on host, matching the flat-array layout + // defined during setup in OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE. + /* Inequality constraints start after equality constraints + Hence the offset + */ + roffset = opflow->nconeq; + coffset = 0; + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); pbpolrajahiopsparse->i_jacineq = (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); pbpolrajahiopsparse->j_jacineq = (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); + pbpolrajahiopsparse->perm_jacineq = + (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); - iRowstart = pbpolrajahiopsparse->i_jacineq; - jColstart = pbpolrajahiopsparse->j_jacineq; + int *iRow_temp = + (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); + int *jCol_temp = + (int *)(h_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); - /* Inequality constraints start after equality constraints - Hence the offset - */ - roffset = opflow->nconeq; - coffset = 0; + BUSParamsRajaHiop *busparams = &pbpolrajahiopsparse->busparams; + GENParamsRajaHiop *genparams = &pbpolrajahiopsparse->genparams; + LINEParamsRajaHiop *lineparams = &pbpolrajahiopsparse->lineparams; + + /* Generator set-point entries */ + if (opflow->has_gensetpoint && opflow->use_agc) { + int delP_col = coffset + pbpolrajahiopsparse->agc_xidx; + for (int g = 0; g < genparams->ngenON; g++) { + if (genparams->isrenewable[g]) + continue; + int base = genparams->ineqjacspgen_idx[g]; + int row0 = roffset + genparams->gineqidxgen[g]; + int Pg_col = coffset + genparams->xidx[g]; + int delPg_col = coffset + genparams->xpdevidx[g]; + + for (int r = 0; r < 2; r++) { + iRow_temp[base + 3 * r + 0] = row0 + r; + jCol_temp[base + 3 * r + 0] = Pg_col; + iRow_temp[base + 3 * r + 1] = row0 + r; + jCol_temp[base + 3 * r + 1] = delPg_col; + iRow_temp[base + 3 * r + 2] = row0 + r; + jCol_temp[base + 3 * r + 2] = delP_col; + } + } + } - ierr = (*opflow->modelops.computeinequalityconstraintjacobian)( - opflow, opflow->X, opflow->Jac_Gi); - CHKERRQ(ierr); + /* FIXED_WITHIN_QBOUNDS voltage entries */ + if (opflow->genbusvoltagetype == FIXED_WITHIN_QBOUNDS) { + for (int ibus = 0; ibus < busparams->nbus; ibus++) { + if (!busparams->ispv[ibus] && !busparams->isref[ibus]) + continue; + int base = busparams->ineqjacsp_idx[ibus]; + int ngen = busparams->ngenONbus[ibus]; + int goff = busparams->genoffset[ibus]; + int row0 = roffset + busparams->gineqidx[ibus]; + int Vm_col = coffset + busparams->xidx[ibus] + 1; + + for (int r = 0; r < 2; r++) { + int rbase = base + r * (ngen + 1); + for (int k = 0; k < ngen; k++) { + iRow_temp[rbase + k] = row0 + r; + jCol_temp[rbase + k] = coffset + genparams->xidx[goff + k] + 1; + } + iRow_temp[rbase + ngen] = row0 + r; + jCol_temp[rbase + ngen] = Vm_col; + } + } + } - ierr = MatGetSize(opflow->Jac_Gi, &nrow, &ncol); - CHKERRQ(ierr); - /* Copy over locations to triplet format */ - for (i = 0; i < nrow; i++) { - ierr = MatGetRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); - CHKERRQ(ierr); - for (j = 0; j < nvals; j++) { - iRowstart[j] = roffset + i; - jColstart[j] = coffset + cols[j]; + /* Line flow entries */ + int has_slack = (int)opflow->allow_lineflow_violation; + int row_stride = 4 + has_slack; + for (int i = 0; i < lineparams->nlinelim; i++) { + int l = lineparams->linelimidx[i]; + int base = lineparams->ineqjacsp_idx[i]; + int row0 = roffset + lineparams->gineqidx[i]; + int thetaf_col = coffset + lineparams->xidxf[l]; + int thetat_col = coffset + lineparams->xidxt[l]; + + for (int r = 0; r < 2; r++) { + int rbase = base + r * row_stride; + iRow_temp[rbase + 0] = row0 + r; + jCol_temp[rbase + 0] = thetaf_col; + iRow_temp[rbase + 1] = row0 + r; + jCol_temp[rbase + 1] = thetaf_col + 1; + iRow_temp[rbase + 2] = row0 + r; + jCol_temp[rbase + 2] = thetat_col; + iRow_temp[rbase + 3] = row0 + r; + jCol_temp[rbase + 3] = thetat_col + 1; + if (has_slack) { + iRow_temp[rbase + 4] = row0 + r; + jCol_temp[rbase + 4] = coffset + lineparams->xslackidx[i] + r; + } } - /* Increment iRow,jCol pointers */ - iRowstart += nvals; - jColstart += nvals; - ierr = MatRestoreRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); - CHKERRQ(ierr); } + // Sort and permute indices + std::vector perm_temp(opflow->nnz_ineqjacsp); + std::iota(perm_temp.begin(), perm_temp.end(), 0); + std::sort(perm_temp.begin(), perm_temp.end(), [&](int i, int j) { + return (iRow_temp[i] != iRow_temp[j]) ? iRow_temp[i] < iRow_temp[j] + : jCol_temp[i] < jCol_temp[j]; + }); + + int *iRow = pbpolrajahiopsparse->i_jacineq; + int *jCol = pbpolrajahiopsparse->j_jacineq; + int *perm = pbpolrajahiopsparse->perm_jacineq; + for (int i = 0; i < opflow->nnz_ineqjacsp; i++) { + iRow[i] = iRow_temp[perm_temp[i]]; + jCol[i] = jCol_temp[perm_temp[i]]; + perm[perm_temp[i]] = i; + } + h_allocator_.deallocate(iRow_temp); + h_allocator_.deallocate(jCol_temp); + // Copy over i_jacineq and j_jacineq arrays to device resmgr.copy(iJacS_dev + opflow->nnz_eqjacsp, pbpolrajahiopsparse->i_jacineq); resmgr.copy(jJacS_dev + opflow->nnz_eqjacsp, pbpolrajahiopsparse->j_jacineq); - ierr = PetscLogEventEnd(opflow->ineqconsjaclogger, 0, 0, 0, 0); - CHKERRQ(ierr); + umpire::Allocator d_allocator_ = resmgr.getAllocator("DEVICE"); + pbpolrajahiopsparse->perm_jacineq_dev = + (int *)(d_allocator_.allocate(opflow->nnz_ineqjacsp * sizeof(int))); + resmgr.copy(pbpolrajahiopsparse->perm_jacineq_dev, + pbpolrajahiopsparse->perm_jacineq); } } @@ -653,7 +729,8 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( No H2D, D2H copies: x_dev is already on device, output goes straight into the ineq portion of MJacS_dev. */ ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE( - opflow, x_dev, MJacS_dev + opflow->nnz_eqjacsp); + opflow, x_dev, pbpolrajahiopsparse->perm_jacineq_dev, + MJacS_dev + opflow->nnz_eqjacsp); ierr = PetscLogEventEnd(opflow->ineqconsjaclogger, 0, 0, 0, 0); CHKERRQ(ierr); diff --git a/tests/unit/opflow/test_ineqjac_gpu.cpp b/tests/unit/opflow/test_ineqjac_gpu.cpp index 0d54a98d..f9b5918d 100644 --- a/tests/unit/opflow/test_ineqjac_gpu.cpp +++ b/tests/unit/opflow/test_ineqjac_gpu.cpp @@ -241,8 +241,25 @@ int main(int argc, char **argv) { resmgr.copy(x_dev, x_sd); #endif + /* Build the sparsity pattern and the permutation used by the kernel */ + int nnz_jac = opflow->nnz_eqjacsp + nnz; + int *iJacS_dev, *jJacS_dev; +#ifdef EXAGO_ENABLE_GPU + iJacS_dev = static_cast(d_allocator.allocate(nnz_jac * sizeof(int))); + jJacS_dev = static_cast(d_allocator.allocate(nnz_jac * sizeof(int))); +#else + iJacS_dev = static_cast(h_allocator.allocate(nnz_jac * sizeof(int))); + jJacS_dev = static_cast(h_allocator.allocate(nnz_jac * sizeof(int))); +#endif + ierr = (*opflow->modelops.computesparseinequalityconstraintjacobianhiop)( + opflow, x_dev, iJacS_dev, jJacS_dev, NULL); + CHKERRQ(ierr); + const int *perm_dev = + reinterpret_cast(opflow->model)->perm_jacineq_dev; + std::cout << " Running RAJA GPU inequality Jacobian kernel..." << std::endl; - ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(opflow, x_dev, gpu_vals_dev); + ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(opflow, x_dev, perm_dev, + gpu_vals_dev); gpu_vals = static_cast(h_allocator.allocate(nnz * sizeof(double))); #ifdef EXAGO_ENABLE_GPU @@ -362,7 +379,8 @@ int main(int argc, char **argv) { #endif auto t0 = std::chrono::high_resolution_clock::now(); for (int iter = 0; iter < niters; iter++) { - ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(opflow, x_dev, bench_dev); + ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(opflow, x_dev, perm_dev, + bench_dev); } #ifdef EXAGO_ENABLE_HIP (void)hipDeviceSynchronize(); @@ -393,8 +411,12 @@ int main(int argc, char **argv) { #ifdef EXAGO_ENABLE_GPU d_allocator.deallocate(x_dev); d_allocator.deallocate(gpu_vals_dev); + d_allocator.deallocate(iJacS_dev); + d_allocator.deallocate(jJacS_dev); #else h_allocator.deallocate(gpu_vals_dev); + h_allocator.deallocate(iJacS_dev); + h_allocator.deallocate(jJacS_dev); #endif ierr = OPFLOWDestroy(&opflow); From 30c34a10a9d5b55594770ed30197a293dda6017b Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Sun, 4 Oct 2026 16:30:20 -0400 Subject: [PATCH 2/3] Clean up GPU Jacobian and Hessian unit tests. --- tests/unit/opflow/CMakeLists.txt | 1 + tests/unit/opflow/test_eqjac_gpu.cpp | 30 +++++--- tests/unit/opflow/test_hess_gpu.cpp | 27 ++++++-- tests/unit/opflow/test_ineqjac_gpu.cpp | 95 ++++++++++---------------- 4 files changed, 78 insertions(+), 75 deletions(-) diff --git a/tests/unit/opflow/CMakeLists.txt b/tests/unit/opflow/CMakeLists.txt index 82957bc0..0a88bee5 100644 --- a/tests/unit/opflow/CMakeLists.txt +++ b/tests/unit/opflow/CMakeLists.txt @@ -77,6 +77,7 @@ if(EXAGO_RUN_TESTS) DEPENDS GPU HIOP_SPARSE + IPOPT RAJA COMMAND ${RUNCMD} diff --git a/tests/unit/opflow/test_eqjac_gpu.cpp b/tests/unit/opflow/test_eqjac_gpu.cpp index dbdf4484..db9b1922 100644 --- a/tests/unit/opflow/test_eqjac_gpu.cpp +++ b/tests/unit/opflow/test_eqjac_gpu.cpp @@ -23,41 +23,47 @@ struct TripletEntry { double val; }; -static void computeReferenceJacobian(OPFLOW opflow, Vec X, - std::vector &entries) { +static PetscErrorCode +computeReferenceJacobian(OPFLOW opflow, Vec X, + std::vector &entries) { PetscErrorCode ierr; ierr = (*opflow->modelops.computeequalityconstraintjacobian)(opflow, X, opflow->Jac_Ge); + CHKERRQ(ierr); PetscInt nrow, ncol; ierr = MatGetSize(opflow->Jac_Ge, &nrow, &ncol); + CHKERRQ(ierr); for (PetscInt i = 0; i < nrow; i++) { PetscInt nvals; const PetscInt *cols; const PetscScalar *vals; ierr = MatGetRow(opflow->Jac_Ge, i, &nvals, &cols, &vals); + CHKERRQ(ierr); for (PetscInt j = 0; j < nvals; j++) { entries.push_back({(int)i, (int)cols[j], vals[j]}); } ierr = MatRestoreRow(opflow->Jac_Ge, i, &nvals, &cols, &vals); + CHKERRQ(ierr); } + return 0; } -static double benchmarkPETSc(OPFLOW opflow, Vec X, int niters) { +static PetscErrorCode benchmarkPETSc(OPFLOW opflow, Vec X, int niters, + double *ms) { PetscErrorCode ierr; - PetscScalar *x_arr; - ierr = VecGetArray(X, &x_arr); auto t0 = Clock::now(); for (int iter = 0; iter < niters; iter++) { ierr = (*opflow->modelops.computeequalityconstraintjacobian)( opflow, X, opflow->Jac_Ge); + CHKERRQ(ierr); } auto t1 = Clock::now(); - ierr = VecRestoreArray(X, &x_arr); - return Ms(t1 - t0).count() / niters; + *ms = Ms(t1 - t0).count() / niters; + return 0; } int main(int argc, char **argv) { @@ -108,7 +114,8 @@ int main(int argc, char **argv) { CHKERRQ(ierr); std::vector ref_entries; - computeReferenceJacobian(opflow_ref, X_ref, ref_entries); + ierr = computeReferenceJacobian(opflow_ref, X_ref, ref_entries); + CHKERRQ(ierr); /* ---------------------------------------------------------------- * Set up OPFLOW with HIOPSPARSE to exercise the GPU path. @@ -301,8 +308,11 @@ int main(int argc, char **argv) { * Benchmark performance * ---------------------------------------------------------------- */ /* Warmup and benchmark PETSc*/ - benchmarkPETSc(opflow_ref, X_ref, 5); - double petsc_ms = benchmarkPETSc(opflow_ref, X_ref, niters); + double petsc_ms; + ierr = benchmarkPETSc(opflow_ref, X_ref, 5, &petsc_ms); + CHKERRQ(ierr); + ierr = benchmarkPETSc(opflow_ref, X_ref, niters, &petsc_ms); + CHKERRQ(ierr); /* Warmup the GPU values kernel */ for (int i = 0; i < 5; i++) { diff --git a/tests/unit/opflow/test_hess_gpu.cpp b/tests/unit/opflow/test_hess_gpu.cpp index 5d28870f..02d3ae6f 100644 --- a/tests/unit/opflow/test_hess_gpu.cpp +++ b/tests/unit/opflow/test_hess_gpu.cpp @@ -23,40 +23,47 @@ struct TripletEntry { double val; }; -static void computeReferenceHessian(OPFLOW opflow, - std::vector &entries) { +static PetscErrorCode +computeReferenceHessian(OPFLOW opflow, std::vector &entries) { PetscErrorCode ierr; ierr = (*opflow->modelops.computehessian)(opflow, opflow->X, opflow->Lambdae, opflow->Lambdai, opflow->Hes); + CHKERRQ(ierr); PetscInt nrow, ncol; ierr = MatGetSize(opflow->Hes, &nrow, &ncol); + CHKERRQ(ierr); for (PetscInt i = 0; i < nrow; i++) { PetscInt nvals; const PetscInt *cols; const PetscScalar *vals; ierr = MatGetRow(opflow->Hes, i, &nvals, &cols, &vals); + CHKERRQ(ierr); for (PetscInt j = 0; j < nvals; j++) { if (cols[j] >= i) { // upper triangle entries.push_back({(int)i, (int)cols[j], vals[j]}); } } ierr = MatRestoreRow(opflow->Hes, i, &nvals, &cols, &vals); + CHKERRQ(ierr); } + return 0; } -static double benchmarkPETSc(OPFLOW opflow, int niters) { +static PetscErrorCode benchmarkPETSc(OPFLOW opflow, int niters, double *ms) { PetscErrorCode ierr; auto t0 = Clock::now(); for (int iter = 0; iter < niters; iter++) { ierr = (*opflow->modelops.computehessian)( opflow, opflow->X, opflow->Lambdae, opflow->Lambdai, opflow->Hes); + CHKERRQ(ierr); } auto t1 = Clock::now(); - return Ms(t1 - t0).count() / niters; + *ms = Ms(t1 - t0).count() / niters; + return 0; } int main(int argc, char **argv) { @@ -111,7 +118,8 @@ int main(int argc, char **argv) { CHKERRQ(ierr); std::vector ref_entries; - computeReferenceHessian(opflow_ref, ref_entries); + ierr = computeReferenceHessian(opflow_ref, ref_entries); + CHKERRQ(ierr); /* ---------------------------------------------------------------- * Set up OPFLOW with HIOPSPARSE to exercise the GPU path. @@ -332,8 +340,11 @@ int main(int argc, char **argv) { * Benchmark performance * ---------------------------------------------------------------- */ /* Warmup and benchmark PETSc*/ - benchmarkPETSc(opflow_ref, 5); - double petsc_ms = benchmarkPETSc(opflow_ref, niters); + double petsc_ms; + ierr = benchmarkPETSc(opflow_ref, 5, &petsc_ms); + CHKERRQ(ierr); + ierr = benchmarkPETSc(opflow_ref, niters, &petsc_ms); + CHKERRQ(ierr); /* Warmup the GPU values kernel */ for (int i = 0; i < 5; i++) { @@ -376,6 +387,8 @@ int main(int argc, char **argv) { ierr = VecRestoreArray(X_gpu, &x_host); CHKERRQ(ierr); + ierr = VecRestoreArray(Lambda_gpu, &lambda_host); + CHKERRQ(ierr); /* ---------------------------------------------------------------- * Print benchmark results diff --git a/tests/unit/opflow/test_ineqjac_gpu.cpp b/tests/unit/opflow/test_ineqjac_gpu.cpp index f9b5918d..2d4db339 100644 --- a/tests/unit/opflow/test_ineqjac_gpu.cpp +++ b/tests/unit/opflow/test_ineqjac_gpu.cpp @@ -174,26 +174,29 @@ int main(int argc, char **argv) { double *ref_vals = static_cast(h_allocator.allocate(nnz * sizeof(double))); + int *ref_rows = static_cast(h_allocator.allocate(nnz * sizeof(int))); + int *ref_cols = static_cast(h_allocator.allocate(nnz * sizeof(int))); PetscInt nrow, ncol; ierr = MatGetSize(opflow->Jac_Gi, &nrow, &ncol); CHKERRQ(ierr); - double *vptr = ref_vals; + int ref_count = 0; for (int i = 0; i < nrow; i++) { PetscInt nvals; const PetscInt *cols; const PetscScalar *vals; ierr = MatGetRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); CHKERRQ(ierr); - for (int j = 0; j < nvals; j++) - vptr[j] = vals[j]; - vptr += nvals; + for (int j = 0; j < nvals; j++, ref_count++) { + ref_rows[ref_count] = opflow->nconeq + i; + ref_cols[ref_count] = cols[j]; + ref_vals[ref_count] = vals[j]; + } ierr = MatRestoreRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); CHKERRQ(ierr); } - int ref_count = (int)(vptr - ref_vals); std::cout << " PETSc extracted " << ref_count << " ineq Jacobian values" << " (expected " << nnz << ")" << std::endl; @@ -257,6 +260,31 @@ int main(int argc, char **argv) { const int *perm_dev = reinterpret_cast(opflow->model)->perm_jacineq_dev; + int *gpu_rows = static_cast(h_allocator.allocate(nnz * sizeof(int))); + int *gpu_cols = static_cast(h_allocator.allocate(nnz * sizeof(int))); + resmgr.copy(gpu_rows, iJacS_dev + opflow->nnz_eqjacsp, nnz * sizeof(int)); + resmgr.copy(gpu_cols, jJacS_dev + opflow->nnz_eqjacsp, nnz * sizeof(int)); + + std::cout << " Comparing " << nnz + << " inequality Jacobian (row, col) pairs..." << std::endl; + int struct_fail = 0; + for (int i = 0; i < nnz; i++) { + if (gpu_rows[i] != ref_rows[i] || gpu_cols[i] != ref_cols[i]) { + std::cout << " MISMATCH ineqjac structure[" << i << "]: PETSc=(" + << ref_rows[i] << ", " << ref_cols[i] << ") GPU=(" + << gpu_rows[i] << ", " << gpu_cols[i] << ")" << std::endl; + struct_fail++; + } + } + fail += struct_fail; + + if (struct_fail == 0) + std::cout << " PASS: All " << nnz + << " inequality Jacobian (row, col) pairs match" << std::endl; + else + std::cout << " FAIL: " << struct_fail << " of " << nnz + << " (row, col) pairs differ" << std::endl; + std::cout << " Running RAJA GPU inequality Jacobian kernel..." << std::endl; ComputeIneqJacValuesGPU_PBPOLRAJAHIOPSPARSE(opflow, x_dev, perm_dev, gpu_vals_dev); @@ -291,61 +319,9 @@ int main(int argc, char **argv) { * ------------------------------------------------------------------ */ { int niters = 1000; - PetscInt bench_nrow, bench_ncol; - ierr = MatGetSize(opflow->Jac_Gi, &bench_nrow, &bench_ncol); - CHKERRQ(ierr); - - double *bench_vals = - static_cast(h_allocator.allocate(nnz * sizeof(double))); - std::cout << "\n=== Performance benchmark (" << niters << " iterations) ===" << std::endl; - /* --- PETSc path: compute + MatGetRow extraction + copy to device --- */ - { - double *bench_dev; - size_t nnz_bytes = nnz * sizeof(double); -#ifdef EXAGO_ENABLE_GPU - bench_dev = static_cast(d_allocator.allocate(nnz_bytes)); -#else - bench_dev = static_cast(h_allocator.allocate(nnz_bytes)); -#endif - - auto t0 = std::chrono::high_resolution_clock::now(); - for (int iter = 0; iter < niters; iter++) { - ierr = (*opflow->modelops.computeinequalityconstraintjacobian)( - opflow, opflow->X, opflow->Jac_Gi); - - double *vp = bench_vals; - for (int i = 0; i < bench_nrow; i++) { - PetscInt nv; - const PetscInt *c; - const PetscScalar *v; - MatGetRow(opflow->Jac_Gi, i, &nv, &c, &v); - for (int j = 0; j < nv; j++) - vp[j] = v[j]; - vp += nv; - MatRestoreRow(opflow->Jac_Gi, i, &nv, &c, &v); - } -#ifdef EXAGO_ENABLE_GPU - resmgr.copy(bench_dev, bench_vals); -#else - memcpy(bench_dev, bench_vals, nnz_bytes); -#endif - } - auto t1 = std::chrono::high_resolution_clock::now(); - double petsc_us = - std::chrono::duration(t1 - t0).count() / niters; - std::cout << " PETSc path (compute + MatGetRow + copy): " << petsc_us - << " us/iter" << std::endl; - -#ifdef EXAGO_ENABLE_GPU - d_allocator.deallocate(bench_dev); -#else - h_allocator.deallocate(bench_dev); -#endif - } - /* --- PETSc path: CPU compute only --- */ { auto t0 = std::chrono::high_resolution_clock::now(); @@ -398,7 +374,6 @@ int main(int argc, char **argv) { #endif } - h_allocator.deallocate(bench_vals); std::cout << "=== End benchmark ===" << std::endl; } @@ -406,6 +381,10 @@ int main(int argc, char **argv) { * Cleanup * ------------------------------------------------------------------ */ h_allocator.deallocate(ref_vals); + h_allocator.deallocate(ref_rows); + h_allocator.deallocate(ref_cols); + h_allocator.deallocate(gpu_rows); + h_allocator.deallocate(gpu_cols); h_allocator.deallocate(gpu_vals); h_allocator.deallocate(x_sd); #ifdef EXAGO_ENABLE_GPU From eb60d1f47e4d8e6342cc0ebdb4a3777974152338 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Wed, 7 Oct 2026 12:23:11 -0400 Subject: [PATCH 3/3] Avoid overflow when extracting the reference inequality Jacobian. --- tests/unit/opflow/test_ineqjac_gpu.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/unit/opflow/test_ineqjac_gpu.cpp b/tests/unit/opflow/test_ineqjac_gpu.cpp index 2d4db339..f0a7225e 100644 --- a/tests/unit/opflow/test_ineqjac_gpu.cpp +++ b/tests/unit/opflow/test_ineqjac_gpu.cpp @@ -189,9 +189,11 @@ int main(int argc, char **argv) { ierr = MatGetRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); CHKERRQ(ierr); for (int j = 0; j < nvals; j++, ref_count++) { - ref_rows[ref_count] = opflow->nconeq + i; - ref_cols[ref_count] = cols[j]; - ref_vals[ref_count] = vals[j]; + if (ref_count < nnz) { // Excess entries counted, reported below + ref_rows[ref_count] = opflow->nconeq + i; + ref_cols[ref_count] = cols[j]; + ref_vals[ref_count] = vals[j]; + } } ierr = MatRestoreRow(opflow->Jac_Gi, i, &nvals, &cols, &vals); CHKERRQ(ierr);