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