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