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) { 71 PetscCall(VecDuplicate(vfull,vreduced)); 72 } 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) { 82 rv[s[i]-flow] = fv[s[i]-flow]; 83 } 84 PetscCall(ISRestoreIndices(is,&s)); 85 PetscCall(VecRestoreArray(vfull,&fv)); 86 PetscCall(VecRestoreArray(*vreduced,&rv)); 87 break; 88 } 89 } 90 PetscFunctionReturn(0); 91 } 92 93 /*@C 94 TaoMatGetSubMat - Gets a submatrix using the IS 95 96 Input Parameters: 97 + M - the full matrix (n x n) 98 . is - the index set for the submatrix (both row and column index sets need to be the same) 99 . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option 100 - subset_type <TAO_SUBSET_SUBVEC,TAO_SUBSET_MASK,TAO_SUBSET_MATRIXFREE> - the method TAO is using for subsetting 101 102 Output Parameter: 103 . Msub - the submatrix 104 105 Level: developer 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(0); 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 Notes: 178 This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3. 179 180 Level: developer 181 @*/ 182 PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol, 183 IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive) 184 { 185 PetscReal wnorm; 186 PetscReal zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0); 187 PetscInt i, n_isl=0, n_isu=0, n_isf=0, n_isa=0, n_isi=0; 188 PetscInt N_isl, N_isu, N_isf, N_isa, N_isi; 189 PetscInt n, low, high, nDiff; 190 PetscInt *isl=NULL, *isu=NULL, *isf=NULL, *isa=NULL, *isi=NULL; 191 const PetscScalar *xl, *xu, *x, *g; 192 MPI_Comm comm = PetscObjectComm((PetscObject)X); 193 194 PetscFunctionBegin; 195 PetscValidHeaderSpecific(X,VEC_CLASSID,1); 196 if (XL) PetscValidHeaderSpecific(XL,VEC_CLASSID,2); 197 if (XU) PetscValidHeaderSpecific(XU,VEC_CLASSID,3); 198 PetscValidHeaderSpecific(G,VEC_CLASSID,4); 199 PetscValidHeaderSpecific(S,VEC_CLASSID,5); 200 PetscValidHeaderSpecific(W,VEC_CLASSID,6); 201 202 if (XL) PetscCheckSameType(X,1,XL,2); 203 if (XU) PetscCheckSameType(X,1,XU,3); 204 PetscCheckSameType(X,1,G,4); 205 PetscCheckSameType(X,1,S,5); 206 PetscCheckSameType(X,1,W,6); 207 if (XL) PetscCheckSameComm(X,1,XL,2); 208 if (XU) PetscCheckSameComm(X,1,XU,3); 209 PetscCheckSameComm(X,1,G,4); 210 PetscCheckSameComm(X,1,S,5); 211 PetscCheckSameComm(X,1,W,6); 212 if (XL) VecCheckSameSize(X,1,XL,2); 213 if (XU) VecCheckSameSize(X,1,XU,3); 214 VecCheckSameSize(X,1,G,4); 215 VecCheckSameSize(X,1,S,5); 216 VecCheckSameSize(X,1,W,6); 217 218 /* Update the tolerance for bound detection (this is based on Bertsekas' method) */ 219 PetscCall(VecCopy(X, W)); 220 PetscCall(VecAXPBY(W, steplen, 1.0, S)); 221 PetscCall(TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W)); 222 PetscCall(VecAXPBY(W, 1.0, -1.0, X)); 223 PetscCall(VecNorm(W, NORM_2, &wnorm)); 224 *bound_tol = PetscMin(*bound_tol, wnorm); 225 226 /* Clear all index sets */ 227 PetscCall(ISDestroy(active_lower)); 228 PetscCall(ISDestroy(active_upper)); 229 PetscCall(ISDestroy(active_fixed)); 230 PetscCall(ISDestroy(active)); 231 PetscCall(ISDestroy(inactive)); 232 233 PetscCall(VecGetOwnershipRange(X, &low, &high)); 234 PetscCall(VecGetLocalSize(X, &n)); 235 if (!XL && !XU) { 236 PetscCall(ISCreateStride(comm,n,low,1,inactive)); 237 PetscFunctionReturn(0); 238 } 239 if (n>0) { 240 PetscCall(VecGetArrayRead(X, &x)); 241 PetscCall(VecGetArrayRead(XL, &xl)); 242 PetscCall(VecGetArrayRead(XU, &xu)); 243 PetscCall(VecGetArrayRead(G, &g)); 244 245 /* Loop over variables and categorize the indexes */ 246 PetscCall(PetscMalloc1(n, &isl)); 247 PetscCall(PetscMalloc1(n, &isu)); 248 PetscCall(PetscMalloc1(n, &isf)); 249 PetscCall(PetscMalloc1(n, &isa)); 250 PetscCall(PetscMalloc1(n, &isi)); 251 for (i=0; i<n; ++i) { 252 if (xl[i] == xu[i]) { 253 /* Fixed variables */ 254 isf[n_isf]=low+i; ++n_isf; 255 isa[n_isa]=low+i; ++n_isa; 256 } else if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + *bound_tol && g[i] > zero) { 257 /* Lower bounded variables */ 258 isl[n_isl]=low+i; ++n_isl; 259 isa[n_isa]=low+i; ++n_isa; 260 } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - *bound_tol && g[i] < zero) { 261 /* Upper bounded variables */ 262 isu[n_isu]=low+i; ++n_isu; 263 isa[n_isa]=low+i; ++n_isa; 264 } else { 265 /* Inactive variables */ 266 isi[n_isi]=low+i; ++n_isi; 267 } 268 } 269 270 PetscCall(VecRestoreArrayRead(X, &x)); 271 PetscCall(VecRestoreArrayRead(XL, &xl)); 272 PetscCall(VecRestoreArrayRead(XU, &xu)); 273 PetscCall(VecRestoreArrayRead(G, &g)); 274 } 275 276 /* Collect global sizes */ 277 PetscCall(MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm)); 278 PetscCall(MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm)); 279 PetscCall(MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm)); 280 PetscCall(MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm)); 281 PetscCall(MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm)); 282 283 /* Create index set for lower bounded variables */ 284 if (N_isl > 0) { 285 PetscCall(ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower)); 286 } else { 287 PetscCall(PetscFree(isl)); 288 } 289 /* Create index set for upper bounded variables */ 290 if (N_isu > 0) { 291 PetscCall(ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper)); 292 } else { 293 PetscCall(PetscFree(isu)); 294 } 295 /* Create index set for fixed variables */ 296 if (N_isf > 0) { 297 PetscCall(ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed)); 298 } else { 299 PetscCall(PetscFree(isf)); 300 } 301 /* Create index set for all actively bounded variables */ 302 if (N_isa > 0) { 303 PetscCall(ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active)); 304 } else { 305 PetscCall(PetscFree(isa)); 306 } 307 /* Create index set for all inactive variables */ 308 if (N_isi > 0) { 309 PetscCall(ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive)); 310 } else { 311 PetscCall(PetscFree(isi)); 312 } 313 PetscFunctionReturn(0); 314 } 315 316 /*@C 317 TaoBoundStep - Ensures the correct zero or adjusted step direction 318 values for active variables. 319 320 Input Parameters: 321 + X - solution vector 322 . XL - lower bound vector 323 . XU - upper bound vector 324 . active_lower - index set for lower bounded active variables 325 . active_upper - index set for lower bounded active variables 326 . active_fixed - index set for fixed active variables 327 - scale - amplification factor for the step that needs to be taken on actively bounded variables 328 329 Output Parameter: 330 . S - step direction to be modified 331 332 Level: developer 333 @*/ 334 PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S) 335 { 336 337 Vec step_lower, step_upper, step_fixed; 338 Vec x_lower, x_upper; 339 Vec bound_lower, bound_upper; 340 341 PetscFunctionBegin; 342 /* Adjust step for variables at the estimated lower bound */ 343 if (active_lower) { 344 PetscCall(VecGetSubVector(S, active_lower, &step_lower)); 345 PetscCall(VecGetSubVector(X, active_lower, &x_lower)); 346 PetscCall(VecGetSubVector(XL, active_lower, &bound_lower)); 347 PetscCall(VecCopy(bound_lower, step_lower)); 348 PetscCall(VecAXPY(step_lower, -1.0, x_lower)); 349 PetscCall(VecScale(step_lower, scale)); 350 PetscCall(VecRestoreSubVector(S, active_lower, &step_lower)); 351 PetscCall(VecRestoreSubVector(X, active_lower, &x_lower)); 352 PetscCall(VecRestoreSubVector(XL, active_lower, &bound_lower)); 353 } 354 355 /* Adjust step for the variables at the estimated upper bound */ 356 if (active_upper) { 357 PetscCall(VecGetSubVector(S, active_upper, &step_upper)); 358 PetscCall(VecGetSubVector(X, active_upper, &x_upper)); 359 PetscCall(VecGetSubVector(XU, active_upper, &bound_upper)); 360 PetscCall(VecCopy(bound_upper, step_upper)); 361 PetscCall(VecAXPY(step_upper, -1.0, x_upper)); 362 PetscCall(VecScale(step_upper, scale)); 363 PetscCall(VecRestoreSubVector(S, active_upper, &step_upper)); 364 PetscCall(VecRestoreSubVector(X, active_upper, &x_upper)); 365 PetscCall(VecRestoreSubVector(XU, active_upper, &bound_upper)); 366 } 367 368 /* Zero out step for fixed variables */ 369 if (active_fixed) { 370 PetscCall(VecGetSubVector(S, active_fixed, &step_fixed)); 371 PetscCall(VecSet(step_fixed, 0.0)); 372 PetscCall(VecRestoreSubVector(S, active_fixed, &step_fixed)); 373 } 374 PetscFunctionReturn(0); 375 } 376 377 /*@C 378 TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance. 379 380 Collective on Vec 381 382 Input Parameters: 383 + X - solution vector 384 . XL - lower bound vector 385 . XU - upper bound vector 386 - bound_tol - absolute tolerance in enforcing the bound 387 388 Output Parameters: 389 + nDiff - total number of vector entries that have been bounded 390 - Xout - modified solution vector satisfying bounds to bound_tol 391 392 Level: developer 393 394 .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR` 395 @*/ 396 PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout) 397 { 398 PetscInt i,n,low,high,nDiff_loc=0; 399 PetscScalar *xout; 400 const PetscScalar *x,*xl,*xu; 401 402 PetscFunctionBegin; 403 PetscValidHeaderSpecific(X,VEC_CLASSID,1); 404 if (XL) PetscValidHeaderSpecific(XL,VEC_CLASSID,2); 405 if (XU) PetscValidHeaderSpecific(XU,VEC_CLASSID,3); 406 PetscValidHeaderSpecific(Xout,VEC_CLASSID,6); 407 if (!XL && !XU) { 408 PetscCall(VecCopy(X,Xout)); 409 *nDiff = 0.0; 410 PetscFunctionReturn(0); 411 } 412 PetscCheckSameType(X,1,XL,2); 413 PetscCheckSameType(X,1,XU,3); 414 PetscCheckSameType(X,1,Xout,6); 415 PetscCheckSameComm(X,1,XL,2); 416 PetscCheckSameComm(X,1,XU,3); 417 PetscCheckSameComm(X,1,Xout,6); 418 VecCheckSameSize(X,1,XL,2); 419 VecCheckSameSize(X,1,XU,3); 420 VecCheckSameSize(X,1,Xout,4); 421 422 PetscCall(VecGetOwnershipRange(X,&low,&high)); 423 PetscCall(VecGetLocalSize(X,&n)); 424 if (n>0) { 425 PetscCall(VecGetArrayRead(X, &x)); 426 PetscCall(VecGetArrayRead(XL, &xl)); 427 PetscCall(VecGetArrayRead(XU, &xu)); 428 PetscCall(VecGetArray(Xout, &xout)); 429 430 for (i=0;i<n;++i) { 431 if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + bound_tol) { 432 xout[i] = xl[i]; ++nDiff_loc; 433 } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - bound_tol) { 434 xout[i] = xu[i]; ++nDiff_loc; 435 } 436 } 437 438 PetscCall(VecRestoreArrayRead(X, &x)); 439 PetscCall(VecRestoreArrayRead(XL, &xl)); 440 PetscCall(VecRestoreArrayRead(XU, &xu)); 441 PetscCall(VecRestoreArray(Xout, &xout)); 442 } 443 PetscCall(MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X))); 444 PetscFunctionReturn(0); 445 } 446