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