From 136c505259a8255aebffb9e7a8b66d0850d6e54a Mon Sep 17 00:00:00 2001 From: William A Perkins Date: Tue, 31 Oct 2023 07:51:35 -0700 Subject: [PATCH] Inequality Jacobian indexes are correct (for lines) --- .../model/power_bal_hiop/paramsrajahiop.cpp | 7 + .../model/power_bal_hiop/paramsrajahiop.h | 2 + .../power_bal_hiop/pbpolrajahiopsparse.cpp | 19 +- .../pbpolrajahiopsparsekernels.cpp | 243 ++++++++++++------ 4 files changed, 188 insertions(+), 83 deletions(-) diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp index 64a18b79..9259a49e 100644 --- a/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp +++ b/src/opflow/model/power_bal_hiop/paramsrajahiop.cpp @@ -235,6 +235,7 @@ int LINEParamsRajaHiop::copy(OPFLOW opflow) { resmgr.copy(gineqidx_dev_, gineqidx); resmgr.copy(gbineqidx_dev_, gbineqidx); resmgr.copy(linelimidx_dev_, linelimidx); + resmgr.copy(jac_ieq_idx_dev_, jac_ieq_idx); } #else Gff_dev_ = Gff; @@ -258,6 +259,7 @@ int LINEParamsRajaHiop::copy(OPFLOW opflow) { gineqidx_dev_ = gineqidx; gbineqidx_dev_ = gbineqidx; linelimidx_dev_ = linelimidx; + jac_ieq_idx_dev_ = jac_ieq_idx; } #endif return 0; @@ -289,6 +291,7 @@ int LINEParamsRajaHiop::destroy(OPFLOW opflow) { h_allocator_.deallocate(gineqidx); h_allocator_.deallocate(gbineqidx); h_allocator_.deallocate(linelimidx); + h_allocator_.deallocate(jac_ieq_idx); } #ifdef EXAGO_ENABLE_GPU @@ -317,6 +320,7 @@ int LINEParamsRajaHiop::destroy(OPFLOW opflow) { d_allocator_.deallocate(gineqidx_dev_); d_allocator_.deallocate(gbineqidx_dev_); d_allocator_.deallocate(linelimidx_dev_); + d_allocator_.deallocate(jac_ieq_idx_dev_); } #endif @@ -369,6 +373,7 @@ int LINEParamsRajaHiop::allocate(OPFLOW opflow) { linelimidx = paramAlloc(h_allocator_, nlinelim); gineqidx = paramAlloc(h_allocator_, nlinelim); gbineqidx = paramAlloc(h_allocator_, nlinelim); + jac_ieq_idx = paramAlloc(h_allocator_, nlinelim); } PetscInt j = 0; @@ -416,6 +421,7 @@ int LINEParamsRajaHiop::allocate(OPFLOW opflow) { gbineqidx[j] = opflow->nconeq + line->startineqloc; gineqidx[j] = line->startineqloc; linelimidx[j] = linei; + jac_ieq_idx[j] = 0; j++; } @@ -450,6 +456,7 @@ int LINEParamsRajaHiop::allocate(OPFLOW opflow) { gineqidx_dev_ = paramAlloc(d_allocator_, nlinelim); gbineqidx_dev_ = paramAlloc(d_allocator_, nlinelim); linelimidx_dev_ = paramAlloc(d_allocator_, nlinelim); + jac_ieq_idx_dev_ = paramAlloc(d_allocator_, nlinelim); } #endif return 0; diff --git a/src/opflow/model/power_bal_hiop/paramsrajahiop.h b/src/opflow/model/power_bal_hiop/paramsrajahiop.h index 73948f7b..da9f922c 100644 --- a/src/opflow/model/power_bal_hiop/paramsrajahiop.h +++ b/src/opflow/model/power_bal_hiop/paramsrajahiop.h @@ -200,6 +200,7 @@ struct LINEParamsRajaHiop { int *bust_idx; /* To bus index */ int *jacf_idx; /* Location number in the sparse Jacobian (from) */ int *jact_idx; /* Location number in the sparse Jacobian (to) */ + int *jac_ieq_idx;/* Location number in sparse inequality Jacobian */ // Device data double *Gff_dev_; /* From side self conductance */ @@ -228,6 +229,7 @@ struct LINEParamsRajaHiop { int *bust_idx_dev_; /* To bus index */ int *jacf_idx_dev_; /* Location number in the sparse Jacobian (from) */ int *jact_idx_dev_; /* Location number in the sparse Jacobian (to) */ + int *jac_ieq_idx_dev_;/* Location number in sparse inequality Jacobian */ int allocate(OPFLOW); int destroy(OPFLOW); diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp index cdaf6874..8e25aad9 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparse.cpp @@ -310,19 +310,25 @@ PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { nnz_eqjac += 4; } + // if there are lines, non-zeros were over counted + if (ps->nline > 0) { + nnz_eqjac -= 8; + } + std::cout << "Equality Jacobian nonzero count: " << nnz_eqjac << std::endl; if (opflow->has_gensetpoint) { - for (int ibus = 0; ibus < ps->nbus; ++ibus) { + for (int ibus = 0, igen = 0; ibus < ps->nbus; ++ibus) { PSBUS bus = &(ps->bus[ibus]); - for (int igen = 0; igen < bus->ngen; ++igen) { + for (int bgen = 0; bgen < bus->ngen; ++bgen) { PSGEN gen; - ierr = PSBUSGetGen(bus, igen, &gen); + ierr = PSBUSGetGen(bus, bgen, &gen); CHKERRQ(ierr); if (!gen->status) continue; - genparams->ineqjacspgen_idx[igen] = nnz_ineqjac; + genparams->ineqjacspgen_idx[igen] = nnz_eqjac + nnz_ineqjac; nnz_ineqjac += 6; + igen++; } } } @@ -331,9 +337,9 @@ PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { for (int ibus = 0; ibus < ps->nbus; ++ibus) { PSBUS bus = &(ps->bus[ibus]); if (bus->ide == PV_BUS || bus->ide == REF_BUS) { - for (int igen = 0; igen < bus->ngen; ++igen) { + for (int bgen = 0; bgen < bus->ngen; ++bgen) { PSGEN gen; - ierr = PSBUSGetGen(bus, igen, &gen); + ierr = PSBUSGetGen(bus, bgen, &gen); CHKERRQ(ierr); if (!gen->status) continue; @@ -346,6 +352,7 @@ PetscErrorCode OPFLOWModelSetUp_PBPOLRAJAHIOPSPARSE(OPFLOW opflow) { if (!opflow->ignore_lineflow_constraints) { for (int iline = 0; iline < opflow->nlinesmon; ++iline) { // PSLINE line = &ps->line[opflow->linesmon[iline]]; + lineparams->jac_ieq_idx[iline] = nnz_eqjac + nnz_ineqjac; nnz_ineqjac += 8; } } diff --git a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp index ef5cf779..487c4dd7 100644 --- a/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp +++ b/src/opflow/model/power_bal_hiop/pbpolrajahiopsparsekernels.cpp @@ -18,6 +18,9 @@ #include "pbpolrajahiopsparsekernels.hpp" #include "pbpolrajahiopsparse.hpp" +static const bool debugmsg(true); +static const bool oldhostway(true); + PetscErrorCode OPFLOWSetInitialGuessArray_PBPOLRAJAHIOPSPARSE(OPFLOW opflow, double *x0_dev) { PetscErrorCode ierr; @@ -471,6 +474,54 @@ PetscErrorCode OPFLOWComputeGradientArray_PBPOLRAJAHIOPSPARSE( PetscFunctionReturn(0); } +// A routine to get the triplet arrays from the device and print them out +static void +PrintTriplets(const std::string& title, const int& n, int *i, int *j, double *v) +{ + auto &resmgr = umpire::ResourceManager::getInstance(); + umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); + + int *itemp(NULL); + + if (i != NULL) { + itemp = (int *)(h_allocator_.allocate(n * sizeof(int))); + resmgr.copy(itemp, i, n*sizeof(int)); + } + + int *jtemp(NULL); + if (j != 0) { + jtemp = (int *)(h_allocator_.allocate(n * sizeof(int))); + resmgr.copy(jtemp, j, n*sizeof(int)); + } + + double *vtemp(NULL); + if (v != NULL) { + vtemp = (double *)(h_allocator_.allocate(n * sizeof(double))); + resmgr.copy(vtemp, v, n*sizeof(double)); + } + + std::cout << title << std::endl; + for (int idx = 0; idx < n; ++idx) { + std::cout << std::setw(5) << idx << " "; + if (itemp != NULL) { + std::cout << std::setw(5) << std::right << itemp[idx] << " "; + } + if (jtemp != NULL) { + std::cout << std::setw(5) << std::right << jtemp[idx]; + } + if (vtemp != NULL) { + std::cout << std::setw(12) << std::right + << std::scientific << std::setprecision(3) + << vtemp[idx]; + } + std::cout << std::endl; + } + h_allocator_.deallocate(itemp); + h_allocator_.deallocate(jtemp); + if (vtemp != NULL) h_allocator_.deallocate(vtemp); +} + + PetscErrorCode OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( OPFLOW opflow, const double *x_dev, int *iJacS_dev, int *jJacS_dev, @@ -480,7 +531,7 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( PetscErrorCode ierr; double *x, *values; PetscInt *iRowstart, *jColstart; - PetscInt roffset, coffset; + PetscInt roffset, coffset, idxoffset; PetscInt nrow, ncol; PetscInt nvals; const PetscInt *cols; @@ -491,11 +542,91 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( PetscFunctionBegin; if (MJacS_dev == NULL) { + + if (debugmsg) + std::cout << "Official Inequality Jacobian nonzero count: " + << opflow->nnz_ineqjacsp << std::endl; + /* Set locations only */ if (opflow->Nconineq) { ierr = PetscLogEventBegin(opflow->ineqconsjaclogger, 0, 0, 0, 0); + /* Inequality constraints start after equality constraints + Hence the offset + */ + roffset = opflow->nconeq; + coffset = 0; + idxoffset = opflow->nnz_eqjacsp; + + resmgr.memset(iJacS_dev + idxoffset, 0, + opflow->nnz_ineqjacsp*sizeof(int)); + resmgr.memset(jJacS_dev + idxoffset, 0, + opflow->nnz_ineqjacsp*sizeof(int)); + + if (!opflow->ignore_lineflow_constraints) { + LINEParamsRajaHiop *lineparams = &pbpolrajahiopsparse->lineparams; + int *jac_ieq_idx = lineparams->jac_ieq_idx_dev_; + int *linelimidx = lineparams->linelimidx_dev_; + int *xidxf = lineparams->xidxf_dev_; + int *xidxt = lineparams->xidxt_dev_; + int *gbineqidx = lineparams->gbineqidx_dev_; + + RAJA::forall( + RAJA::RangeSegment(0, lineparams->nlinelim), + RAJA_LAMBDA(RAJA::Index_type i) { + int iline(linelimidx[i]); + int offset(0); + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i]; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxf[iline]; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i]; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxf[iline] + 1; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i]; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxt[iline]; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i]; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxt[iline] + 1; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i] + 1; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxf[iline]; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i] + 1; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxf[iline] + 1; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i] + 1; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxt[iline]; + offset++; + + iJacS_dev[jac_ieq_idx[i] + offset] = gbineqidx[i] + 1; + jJacS_dev[jac_ieq_idx[i] + offset] = xidxt[iline] + 1; + }); + + } + + if (debugmsg) + PrintTriplets("Nonzero indexes for Inequality Constraint Jacobian (GPU):", + opflow->nnz_ineqjacsp, + iJacS_dev + opflow->nnz_eqjacsp, + jJacS_dev + opflow->nnz_eqjacsp, + NULL); + + if (oldhostway) { + + /* Inequality constraints start after equality constraints + Hence the offset + */ + roffset = opflow->nconeq; + coffset = 0; + // Create arrays on host to store i,j, and val arrays umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); @@ -509,12 +640,6 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( iRowstart = pbpolrajahiopsparse->i_jacineq; jColstart = pbpolrajahiopsparse->j_jacineq; - /* Inequality constraints start after equality constraints - Hence the offset - */ - roffset = opflow->nconeq; - coffset = 0; - ierr = (*opflow->modelops.computeinequalityconstraintjacobian)( opflow, opflow->X, opflow->Jac_Gi); CHKERRQ(ierr); @@ -536,30 +661,30 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( CHKERRQ(ierr); } - // Dump out the matrix indexes as a check - std::cout << "Nonzero indexes for Inequality Constraint Jacobian: " - << opflow->nnz_ineqjacsp - << std::endl; - for (int idx = 0; idx < opflow->nnz_ineqjacsp; ++idx) { - std::cout << std::setw(5) << idx << " " - << std::setw(5) << pbpolrajahiopsparse->i_jacineq[idx] << " " - << std::setw(5) << pbpolrajahiopsparse->j_jacineq[idx] << std::endl; - } - // 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); + if (debugmsg) + PrintTriplets("Nonzero indexes for Inequality Constraint Jacobian:", + opflow->nnz_ineqjacsp, + iJacS_dev + opflow->nnz_eqjacsp, + jJacS_dev + opflow->nnz_eqjacsp, + NULL); + ierr = PetscLogEventEnd(opflow->ineqconsjaclogger, 0, 0, 0, 0); CHKERRQ(ierr); + } } } else { if (opflow->Nconineq) { ierr = PetscLogEventBegin(opflow->ineqconsjaclogger, 0, 0, 0, 0); CHKERRQ(ierr); + if (oldhostway) { + ierr = VecGetArray(opflow->X, &x); CHKERRQ(ierr); @@ -595,8 +720,16 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( resmgr.copy(MJacS_dev + opflow->nnz_eqjacsp, pbpolrajahiopsparse->val_jacineq); + if (debugmsg) + PrintTriplets("Inequality Constraint Jacobian:", + opflow->nnz_ineqjacsp, + (iJacS_dev == NULL ? NULL : iJacS_dev + opflow->nnz_eqjacsp), + (jJacS_dev == NULL ? NULL : jJacS_dev + opflow->nnz_eqjacsp), + MJacS_dev + opflow->nnz_eqjacsp); + ierr = PetscLogEventEnd(opflow->ineqconsjaclogger, 0, 0, 0, 0); CHKERRQ(ierr); + } } } @@ -604,54 +737,6 @@ OPFLOWComputeSparseInequalityConstraintJacobian_PBPOLRAJAHIOPSPARSE( } -// A routine to get the triplet arrays from the device and print them out -static void -PrintTriplets(const std::string& title, const int& n, int *i, int *j, double *v) -{ - auto &resmgr = umpire::ResourceManager::getInstance(); - umpire::Allocator h_allocator_ = resmgr.getAllocator("HOST"); - - int *itemp(NULL); - - if (i != NULL) { - itemp = (int *)(h_allocator_.allocate(n * sizeof(int))); - resmgr.copy(itemp, i, n*sizeof(int)); - } - - int *jtemp(NULL); - if (j != 0) { - jtemp = (int *)(h_allocator_.allocate(n * sizeof(int))); - resmgr.copy(jtemp, j, n*sizeof(int)); - } - - double *vtemp(NULL); - if (v != NULL) { - vtemp = (double *)(h_allocator_.allocate(n * sizeof(double))); - resmgr.copy(vtemp, v, n*sizeof(double)); - } - - std::cout << title << std::endl; - for (int idx = 0; idx < n; ++idx) { - std::cout << std::setw(5) << idx << " "; - if (itemp != NULL) { - std::cout << std::setw(5) << std::right << itemp[idx] << " "; - } - if (jtemp != NULL) { - std::cout << std::setw(5) << std::right << jtemp[idx]; - } - if (vtemp != NULL) { - std::cout << std::setw(12) << std::right - << std::scientific << std::setprecision(3) - << vtemp[idx]; - } - std::cout << std::endl; - } - h_allocator_.deallocate(itemp); - h_allocator_.deallocate(jtemp); - if (vtemp != NULL) h_allocator_.deallocate(vtemp); -} - - PetscErrorCode OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( OPFLOW opflow, const double *x_dev, int *iJacS_dev, int *jJacS_dev, @@ -685,6 +770,10 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( if (MJacS_dev == NULL) { + if (debugmsg) + std::cout << "Official Equality Jacobian nonzero count: " + << opflow->nnz_eqjacsp << std::endl; + /* Set locations only */ resmgr.memset(iJacS_dev, 0, opflow->nnz_eqjacsp*sizeof(int)); @@ -698,7 +787,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( int *b_jacsq_idx = busparams->jacsq_idx_dev_; /* Bus */ - std::cout << "Begin with buses" << std::endl; + if (debugmsg) std::cout << "Begin with buses" << std::endl; RAJA::forall( RAJA::RangeSegment(0, busparams->nbus), RAJA_LAMBDA(RAJA::Index_type i) { @@ -714,7 +803,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( }); if (opflow->include_powerimbalance_variables) { - std::cout << "Bus power imbalance variables" << std::endl; + if (debugmsg) std::cout << "Bus power imbalance variables" << std::endl; RAJA::forall( RAJA::RangeSegment(0, busparams->nbus), RAJA_LAMBDA(RAJA::Index_type i) { @@ -732,7 +821,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( /* generation contributions */ - std::cout << "Generators " << std::endl; + if (debugmsg) std::cout << "Generators " << std::endl; int *g_gidxbus = genparams->gidxbus_dev_; int *g_xidx = genparams->xidx_dev_; @@ -752,7 +841,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( if (opflow->include_loadloss_variables) { - std::cout << "Load Loss" << std::endl; + if (debugmsg) std::cout << "Load Loss" << std::endl; int *l_gidx = loadparams->gidx_dev_; int *l_xidx = loadparams->xidx_dev_; int *l_jacsp_idx = loadparams->jacsp_idx_dev_; @@ -769,7 +858,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( /* Connected lines */ - std::cout << "Connected Lines" << std::endl; + if (debugmsg) std::cout << "Connected Lines" << std::endl; int *xidxf = lineparams->xidxf_dev_; int *xidxt = lineparams->xidxt_dev_; @@ -828,7 +917,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( if (opflow->has_gensetpoint) { - std::cout << "Generator set point" << std::endl; + if (debugmsg) std::cout << "Generator set point" << std::endl; int *eqjacspgen_idx = genparams->eqjacspgen_idx_dev_; int *g_geqidxgen = genparams->geqidxgen_dev_; int *g_xidx = genparams->xidx_dev_; @@ -853,11 +942,11 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( }); } - if (1) + if (debugmsg) PrintTriplets("Non-zero indexes for Equality Constraint Jacobian (GPU):", opflow->nnz_eqjacsp, iJacS_dev, jJacS_dev, NULL); - if (1) { + if (oldhostway) { roffset = 0; coffset = 0; @@ -897,7 +986,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( resmgr.copy(iJacS_dev, pbpolrajahiopsparse->i_jaceq); resmgr.copy(jJacS_dev, pbpolrajahiopsparse->j_jaceq); - if (1) + if (debugmsg) PrintTriplets("Non-zero indexes for Equality Constraint Jacobian:", opflow->nnz_eqjacsp, iJacS_dev, jJacS_dev, NULL); } @@ -1082,12 +1171,12 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( Vmt * (-Btf[i] * cos(thetatf) + Gtf[i] * sin(thetatf)); }); - if (1) + if (debugmsg) PrintTriplets("Equality Constraint Jacobian (GPU):", opflow->nnz_eqjacsp, iJacS_dev, jJacS_dev, MJacS_dev); - if (1) { + if (oldhostway) { ierr = VecGetArray(opflow->X, &x); CHKERRQ(ierr); @@ -1123,7 +1212,7 @@ OPFLOWComputeSparseEqualityConstraintJacobian_PBPOLRAJAHIOPSPARSE( // Copy over val_ineq to device resmgr.copy(MJacS_dev, pbpolrajahiopsparse->val_jaceq); - if (1) + if (debugmsg) PrintTriplets("Equality Constraint Jacobian:", opflow->nnz_eqjacsp, iJacS_dev, jJacS_dev, MJacS_dev); }