Compute the HIOPSPARSEGPU inequality constraints Jacobian indices without PETSc

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
 }

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
 };

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<PbpolModelRajaHiop *>(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;
           }
         });
   }

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);
 
 /**

src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp

@@ -578,67 +578,143 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE(
   PbpolModelRajaHiop *pbpolrajahiopsparse =
       reinterpret_cast<PbpolModelRajaHiop *>(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;
 
   // 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<int> 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);

tests/unit/opflow/CMakeLists.txt

@@ -77,6 +77,7 @@ if(EXAGO_RUN_TESTS)
     DEPENDS
     GPU
     HIOP_SPARSE
+    IPOPT
     RAJA
     COMMAND
     ${RUNCMD}