1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <petsc/private/vecimpl.h> 3 #include <petsc/private/taoimpl.h> 4 #include <../src/tao/matrix/submatfree.h> 5 6 /*@C 7 TaoVecGetSubVec - Gets a subvector using the IS 8 9 Input Parameters: 10 + vfull - the full matrix 11 . is - the index set for the subvector 12 . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, TAO_SUBSET_MATRIXFREE) 13 - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE) 14 15 Output Parameter: 16 . vreduced - the subvector 17 18 Notes: 19 maskvalue should usually be 0.0, unless a pointwise divide will be used. 20 21 Level: developer 22 @*/ 23 PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced) { 24 PetscInt nfull, nreduced, nreduced_local, rlow, rhigh, flow, fhigh; 25 PetscInt i, nlocal; 26 PetscReal *fv, *rv; 27 const PetscInt *s; 28 IS ident; 29 VecType vtype; 30 VecScatter scatter; 31 MPI_Comm comm; 32 33 PetscFunctionBegin; 34 PetscValidHeaderSpecific(vfull, VEC_CLASSID, 1); 35 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 36 37 PetscCall(VecGetSize(vfull, &nfull)); 38 PetscCall(ISGetSize(is, &nreduced)); 39 40 if (nreduced == nfull) { 41 PetscCall(VecDestroy(vreduced)); 42 PetscCall(VecDuplicate(vfull, vreduced)); 43 PetscCall(VecCopy(vfull, *vreduced)); 44 } else { 45 switch (reduced_type) { 46 case TAO_SUBSET_SUBVEC: 47 PetscCall(VecGetType(vfull, &vtype)); 48 PetscCall(VecGetOwnershipRange(vfull, &flow, &fhigh)); 49 PetscCall(ISGetLocalSize(is, &nreduced_local)); 50 PetscCall(PetscObjectGetComm((PetscObject)vfull, &comm)); 51 if (*vreduced) PetscCall(VecDestroy(vreduced)); 52 PetscCall(VecCreate(comm, vreduced)); 53 PetscCall(VecSetType(*vreduced, vtype)); 54 55 PetscCall(VecSetSizes(*vreduced, nreduced_local, nreduced)); 56 PetscCall(VecGetOwnershipRange(*vreduced, &rlow, &rhigh)); 57 PetscCall(ISCreateStride(comm, nreduced_local, rlow, 1, &ident)); 58 PetscCall(VecScatterCreate(vfull, is, *vreduced, ident, &scatter)); 59 PetscCall(VecScatterBegin(scatter, vfull, *vreduced, INSERT_VALUES, SCATTER_FORWARD)); 60 PetscCall(VecScatterEnd(scatter, vfull, *vreduced, INSERT_VALUES, SCATTER_FORWARD)); 61 PetscCall(VecScatterDestroy(&scatter)); 62 PetscCall(ISDestroy(&ident)); 63 break; 64 65 case TAO_SUBSET_MASK: 66 case TAO_SUBSET_MATRIXFREE: 67 /* vr[i] = vf[i] if i in is 68 vr[i] = 0 otherwise */ 69 if (!*vreduced) PetscCall(VecDuplicate(vfull, vreduced)); 70 71 PetscCall(VecSet(*vreduced, maskvalue)); 72 PetscCall(ISGetLocalSize(is, &nlocal)); 73 PetscCall(VecGetOwnershipRange(vfull, &flow, &fhigh)); 74 PetscCall(VecGetArray(vfull, &fv)); 75 PetscCall(VecGetArray(*vreduced, &rv)); 76 PetscCall(ISGetIndices(is, &s)); 77 PetscCheck(nlocal <= (fhigh - flow), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "IS local size %" PetscInt_FMT " > Vec local size %" PetscInt_FMT, nlocal, fhigh - flow); 78 for (i = 0; i < nlocal; ++i) rv[s[i] - flow] = fv[s[i] - flow]; 79 PetscCall(ISRestoreIndices(is, &s)); 80 PetscCall(VecRestoreArray(vfull, &fv)); 81 PetscCall(VecRestoreArray(*vreduced, &rv)); 82 break; 83 } 84 } 85 PetscFunctionReturn(0); 86 } 87 88 /*@C 89 TaoMatGetSubMat - Gets a submatrix using the IS 90 91 Input Parameters: 92 + M - the full matrix (n x n) 93 . is - the index set for the submatrix (both row and column index sets need to be the same) 94 . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option 95 - subset_type <TAO_SUBSET_SUBVEC,TAO_SUBSET_MASK,TAO_SUBSET_MATRIXFREE> - the method TAO is using for subsetting 96 97 Output Parameter: 98 . Msub - the submatrix 99 100 Level: developer 101 @*/ 102 PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub) { 103 IS iscomp; 104 PetscBool flg = PETSC_TRUE; 105 106 PetscFunctionBegin; 107 PetscValidHeaderSpecific(M, MAT_CLASSID, 1); 108 PetscValidHeaderSpecific(is, IS_CLASSID, 2); 109 PetscCall(MatDestroy(Msub)); 110 switch (subset_type) { 111 case TAO_SUBSET_SUBVEC: PetscCall(MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub)); break; 112 113 case TAO_SUBSET_MASK: 114 /* Get Reduced Hessian 115 Msub[i,j] = M[i,j] if i,j in Free_Local or i==j 116 Msub[i,j] = 0 if i!=j and i or j not in Free_Local 117 */ 118 PetscObjectOptionsBegin((PetscObject)M); 119 PetscCall(PetscOptionsBool("-overwrite_hessian", "modify the existing hessian matrix when computing submatrices", "TaoSubsetType", flg, &flg, NULL)); 120 PetscOptionsEnd(); 121 if (flg) { 122 PetscCall(MatDuplicate(M, MAT_COPY_VALUES, Msub)); 123 } else { 124 /* Act on hessian directly (default) */ 125 PetscCall(PetscObjectReference((PetscObject)M)); 126 *Msub = M; 127 } 128 /* Save the diagonal to temporary vector */ 129 PetscCall(MatGetDiagonal(*Msub, v1)); 130 131 /* Zero out rows and columns */ 132 PetscCall(ISComplementVec(is, v1, &iscomp)); 133 134 /* Use v1 instead of 0 here because of PETSc bug */ 135 PetscCall(MatZeroRowsColumnsIS(*Msub, iscomp, 1.0, v1, v1)); 136 137 PetscCall(ISDestroy(&iscomp)); 138 break; 139 case TAO_SUBSET_MATRIXFREE: 140 PetscCall(ISComplementVec(is, v1, &iscomp)); 141 PetscCall(MatCreateSubMatrixFree(M, iscomp, iscomp, Msub)); 142 PetscCall(ISDestroy(&iscomp)); 143 break; 144 } 145 PetscFunctionReturn(0); 146 } 147 148 /*@C 149 TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper 150 bounds, as well as fixed variables where lower and upper bounds equal each other. 151 152 Input Parameters: 153 + X - solution vector 154 . XL - lower bound vector 155 . XU - upper bound vector 156 . G - unprojected gradient 157 . S - step direction with which the active bounds will be estimated 158 . W - work vector of type and size of X 159 - steplen - the step length at which the active bounds will be estimated (needs to be conservative) 160 161 Output Parameters: 162 + bound_tol - tolerance for the bound estimation 163 . active_lower - index set for active variables at the lower bound 164 . active_upper - index set for active variables at the upper bound 165 . active_fixed - index set for fixed variables 166 . active - index set for all active variables 167 - inactive - complementary index set for inactive variables 168 169 Notes: 170 This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3. 171 172 Level: developer 173 @*/ 174 PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol, IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive) { 175 PetscReal wnorm; 176 PetscReal zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0); 177 PetscInt i, n_isl = 0, n_isu = 0, n_isf = 0, n_isa = 0, n_isi = 0; 178 PetscInt N_isl, N_isu, N_isf, N_isa, N_isi; 179 PetscInt n, low, high, nDiff; 180 PetscInt *isl = NULL, *isu = NULL, *isf = NULL, *isa = NULL, *isi = NULL; 181 const PetscScalar *xl, *xu, *x, *g; 182 MPI_Comm comm = PetscObjectComm((PetscObject)X); 183 184 PetscFunctionBegin; 185 PetscValidHeaderSpecific(X, VEC_CLASSID, 1); 186 if (XL) PetscValidHeaderSpecific(XL, VEC_CLASSID, 2); 187 if (XU) PetscValidHeaderSpecific(XU, VEC_CLASSID, 3); 188 PetscValidHeaderSpecific(G, VEC_CLASSID, 4); 189 PetscValidHeaderSpecific(S, VEC_CLASSID, 5); 190 PetscValidHeaderSpecific(W, VEC_CLASSID, 6); 191 192 if (XL) PetscCheckSameType(X, 1, XL, 2); 193 if (XU) PetscCheckSameType(X, 1, XU, 3); 194 PetscCheckSameType(X, 1, G, 4); 195 PetscCheckSameType(X, 1, S, 5); 196 PetscCheckSameType(X, 1, W, 6); 197 if (XL) PetscCheckSameComm(X, 1, XL, 2); 198 if (XU) PetscCheckSameComm(X, 1, XU, 3); 199 PetscCheckSameComm(X, 1, G, 4); 200 PetscCheckSameComm(X, 1, S, 5); 201 PetscCheckSameComm(X, 1, W, 6); 202 if (XL) VecCheckSameSize(X, 1, XL, 2); 203 if (XU) VecCheckSameSize(X, 1, XU, 3); 204 VecCheckSameSize(X, 1, G, 4); 205 VecCheckSameSize(X, 1, S, 5); 206 VecCheckSameSize(X, 1, W, 6); 207 208 /* Update the tolerance for bound detection (this is based on Bertsekas' method) */ 209 PetscCall(VecCopy(X, W)); 210 PetscCall(VecAXPBY(W, steplen, 1.0, S)); 211 PetscCall(TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W)); 212 PetscCall(VecAXPBY(W, 1.0, -1.0, X)); 213 PetscCall(VecNorm(W, NORM_2, &wnorm)); 214 *bound_tol = PetscMin(*bound_tol, wnorm); 215 216 /* Clear all index sets */ 217 PetscCall(ISDestroy(active_lower)); 218 PetscCall(ISDestroy(active_upper)); 219 PetscCall(ISDestroy(active_fixed)); 220 PetscCall(ISDestroy(active)); 221 PetscCall(ISDestroy(inactive)); 222 223 PetscCall(VecGetOwnershipRange(X, &low, &high)); 224 PetscCall(VecGetLocalSize(X, &n)); 225 if (!XL && !XU) { 226 PetscCall(ISCreateStride(comm, n, low, 1, inactive)); 227 PetscFunctionReturn(0); 228 } 229 if (n > 0) { 230 PetscCall(VecGetArrayRead(X, &x)); 231 PetscCall(VecGetArrayRead(XL, &xl)); 232 PetscCall(VecGetArrayRead(XU, &xu)); 233 PetscCall(VecGetArrayRead(G, &g)); 234 235 /* Loop over variables and categorize the indexes */ 236 PetscCall(PetscMalloc1(n, &isl)); 237 PetscCall(PetscMalloc1(n, &isu)); 238 PetscCall(PetscMalloc1(n, &isf)); 239 PetscCall(PetscMalloc1(n, &isa)); 240 PetscCall(PetscMalloc1(n, &isi)); 241 for (i = 0; i < n; ++i) { 242 if (xl[i] == xu[i]) { 243 /* Fixed variables */ 244 isf[n_isf] = low + i; 245 ++n_isf; 246 isa[n_isa] = low + i; 247 ++n_isa; 248 } else if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + *bound_tol && g[i] > zero) { 249 /* Lower bounded variables */ 250 isl[n_isl] = low + i; 251 ++n_isl; 252 isa[n_isa] = low + i; 253 ++n_isa; 254 } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - *bound_tol && g[i] < zero) { 255 /* Upper bounded variables */ 256 isu[n_isu] = low + i; 257 ++n_isu; 258 isa[n_isa] = low + i; 259 ++n_isa; 260 } else { 261 /* Inactive variables */ 262 isi[n_isi] = low + i; 263 ++n_isi; 264 } 265 } 266 267 PetscCall(VecRestoreArrayRead(X, &x)); 268 PetscCall(VecRestoreArrayRead(XL, &xl)); 269 PetscCall(VecRestoreArrayRead(XU, &xu)); 270 PetscCall(VecRestoreArrayRead(G, &g)); 271 } 272 273 /* Collect global sizes */ 274 PetscCall(MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm)); 275 PetscCall(MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm)); 276 PetscCall(MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm)); 277 PetscCall(MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm)); 278 PetscCall(MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm)); 279 280 /* Create index set for lower bounded variables */ 281 if (N_isl > 0) { 282 PetscCall(ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower)); 283 } else { 284 PetscCall(PetscFree(isl)); 285 } 286 /* Create index set for upper bounded variables */ 287 if (N_isu > 0) { 288 PetscCall(ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper)); 289 } else { 290 PetscCall(PetscFree(isu)); 291 } 292 /* Create index set for fixed variables */ 293 if (N_isf > 0) { 294 PetscCall(ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed)); 295 } else { 296 PetscCall(PetscFree(isf)); 297 } 298 /* Create index set for all actively bounded variables */ 299 if (N_isa > 0) { 300 PetscCall(ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active)); 301 } else { 302 PetscCall(PetscFree(isa)); 303 } 304 /* Create index set for all inactive variables */ 305 if (N_isi > 0) { 306 PetscCall(ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive)); 307 } else { 308 PetscCall(PetscFree(isi)); 309 } 310 PetscFunctionReturn(0); 311 } 312 313 /*@C 314 TaoBoundStep - Ensures the correct zero or adjusted step direction 315 values for active variables. 316 317 Input Parameters: 318 + X - solution vector 319 . XL - lower bound vector 320 . XU - upper bound vector 321 . active_lower - index set for lower bounded active variables 322 . active_upper - index set for lower bounded active variables 323 . active_fixed - index set for fixed active variables 324 - scale - amplification factor for the step that needs to be taken on actively bounded variables 325 326 Output Parameter: 327 . S - step direction to be modified 328 329 Level: developer 330 @*/ 331 PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S) { 332 Vec step_lower, step_upper, step_fixed; 333 Vec x_lower, x_upper; 334 Vec bound_lower, bound_upper; 335 336 PetscFunctionBegin; 337 /* Adjust step for variables at the estimated lower bound */ 338 if (active_lower) { 339 PetscCall(VecGetSubVector(S, active_lower, &step_lower)); 340 PetscCall(VecGetSubVector(X, active_lower, &x_lower)); 341 PetscCall(VecGetSubVector(XL, active_lower, &bound_lower)); 342 PetscCall(VecCopy(bound_lower, step_lower)); 343 PetscCall(VecAXPY(step_lower, -1.0, x_lower)); 344 PetscCall(VecScale(step_lower, scale)); 345 PetscCall(VecRestoreSubVector(S, active_lower, &step_lower)); 346 PetscCall(VecRestoreSubVector(X, active_lower, &x_lower)); 347 PetscCall(VecRestoreSubVector(XL, active_lower, &bound_lower)); 348 } 349 350 /* Adjust step for the variables at the estimated upper bound */ 351 if (active_upper) { 352 PetscCall(VecGetSubVector(S, active_upper, &step_upper)); 353 PetscCall(VecGetSubVector(X, active_upper, &x_upper)); 354 PetscCall(VecGetSubVector(XU, active_upper, &bound_upper)); 355 PetscCall(VecCopy(bound_upper, step_upper)); 356 PetscCall(VecAXPY(step_upper, -1.0, x_upper)); 357 PetscCall(VecScale(step_upper, scale)); 358 PetscCall(VecRestoreSubVector(S, active_upper, &step_upper)); 359 PetscCall(VecRestoreSubVector(X, active_upper, &x_upper)); 360 PetscCall(VecRestoreSubVector(XU, active_upper, &bound_upper)); 361 } 362 363 /* Zero out step for fixed variables */ 364 if (active_fixed) { 365 PetscCall(VecGetSubVector(S, active_fixed, &step_fixed)); 366 PetscCall(VecSet(step_fixed, 0.0)); 367 PetscCall(VecRestoreSubVector(S, active_fixed, &step_fixed)); 368 } 369 PetscFunctionReturn(0); 370 } 371 372 /*@C 373 TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance. 374 375 Collective on Vec 376 377 Input Parameters: 378 + X - solution vector 379 . XL - lower bound vector 380 . XU - upper bound vector 381 - bound_tol - absolute tolerance in enforcing the bound 382 383 Output Parameters: 384 + nDiff - total number of vector entries that have been bounded 385 - Xout - modified solution vector satisfying bounds to bound_tol 386 387 Level: developer 388 389 .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR` 390 @*/ 391 PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout) { 392 PetscInt i, n, low, high, nDiff_loc = 0; 393 PetscScalar *xout; 394 const PetscScalar *x, *xl, *xu; 395 396 PetscFunctionBegin; 397 PetscValidHeaderSpecific(X, VEC_CLASSID, 1); 398 if (XL) PetscValidHeaderSpecific(XL, VEC_CLASSID, 2); 399 if (XU) PetscValidHeaderSpecific(XU, VEC_CLASSID, 3); 400 PetscValidHeaderSpecific(Xout, VEC_CLASSID, 6); 401 if (!XL && !XU) { 402 PetscCall(VecCopy(X, Xout)); 403 *nDiff = 0.0; 404 PetscFunctionReturn(0); 405 } 406 PetscCheckSameType(X, 1, XL, 2); 407 PetscCheckSameType(X, 1, XU, 3); 408 PetscCheckSameType(X, 1, Xout, 6); 409 PetscCheckSameComm(X, 1, XL, 2); 410 PetscCheckSameComm(X, 1, XU, 3); 411 PetscCheckSameComm(X, 1, Xout, 6); 412 VecCheckSameSize(X, 1, XL, 2); 413 VecCheckSameSize(X, 1, XU, 3); 414 VecCheckSameSize(X, 1, Xout, 4); 415 416 PetscCall(VecGetOwnershipRange(X, &low, &high)); 417 PetscCall(VecGetLocalSize(X, &n)); 418 if (n > 0) { 419 PetscCall(VecGetArrayRead(X, &x)); 420 PetscCall(VecGetArrayRead(XL, &xl)); 421 PetscCall(VecGetArrayRead(XU, &xu)); 422 PetscCall(VecGetArray(Xout, &xout)); 423 424 for (i = 0; i < n; ++i) { 425 if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + bound_tol) { 426 xout[i] = xl[i]; 427 ++nDiff_loc; 428 } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - bound_tol) { 429 xout[i] = xu[i]; 430 ++nDiff_loc; 431 } 432 } 433 434 PetscCall(VecRestoreArrayRead(X, &x)); 435 PetscCall(VecRestoreArrayRead(XL, &xl)); 436 PetscCall(VecRestoreArrayRead(XU, &xu)); 437 PetscCall(VecRestoreArray(Xout, &xout)); 438 } 439 PetscCall(MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X))); 440 PetscFunctionReturn(0); 441 } 442