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