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