1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <petsc/private/taoimpl.h> 3 #include <../src/tao/matrix/submatfree.h> 4 5 /*@C 6 TaoVecGetSubVec - Gets a subvector using the IS 7 8 Input Parameters: 9 + vfull - the full matrix 10 . is - the index set for the subvector 11 . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, TAO_SUBSET_MATRIXFREE) 12 - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE) 13 14 Output Parameters: 15 . vreduced - the subvector 16 17 Notes: 18 maskvalue should usually be 0.0, unless a pointwise divide will be used. 19 20 @*/ 21 PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced) 22 { 23 PetscErrorCode ierr; 24 PetscInt nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh; 25 PetscInt i,nlocal; 26 PetscReal *fv,*rv; 27 const PetscInt *s; 28 IS ident; 29 VecType vtype; 30 VecScatter scatter; 31 MPI_Comm comm; 32 33 PetscFunctionBegin; 34 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1); 35 PetscValidHeaderSpecific(is,IS_CLASSID,2); 36 37 ierr = VecGetSize(vfull, &nfull);CHKERRQ(ierr); 38 ierr = ISGetSize(is, &nreduced);CHKERRQ(ierr); 39 40 if (nreduced == nfull) { 41 ierr = VecDestroy(vreduced);CHKERRQ(ierr); 42 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr); 43 ierr = VecCopy(vfull,*vreduced);CHKERRQ(ierr); 44 } else { 45 switch (reduced_type) { 46 case TAO_SUBSET_SUBVEC: 47 ierr = VecGetType(vfull,&vtype);CHKERRQ(ierr); 48 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr); 49 ierr = ISGetLocalSize(is,&nreduced_local);CHKERRQ(ierr); 50 ierr = PetscObjectGetComm((PetscObject)vfull,&comm);CHKERRQ(ierr); 51 if (*vreduced) { 52 ierr = VecDestroy(vreduced);CHKERRQ(ierr); 53 } 54 ierr = VecCreate(comm,vreduced);CHKERRQ(ierr); 55 ierr = VecSetType(*vreduced,vtype);CHKERRQ(ierr); 56 57 ierr = VecSetSizes(*vreduced,nreduced_local,nreduced);CHKERRQ(ierr); 58 ierr = VecGetOwnershipRange(*vreduced,&rlow,&rhigh);CHKERRQ(ierr); 59 ierr = ISCreateStride(comm,nreduced_local,rlow,1,&ident);CHKERRQ(ierr); 60 ierr = VecScatterCreate(vfull,is,*vreduced,ident,&scatter);CHKERRQ(ierr); 61 ierr = VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 62 ierr = VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 63 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 64 ierr = ISDestroy(&ident);CHKERRQ(ierr); 65 break; 66 67 case TAO_SUBSET_MASK: 68 case TAO_SUBSET_MATRIXFREE: 69 /* vr[i] = vf[i] if i in is 70 vr[i] = 0 otherwise */ 71 if (!*vreduced) { 72 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr); 73 } 74 75 ierr = VecSet(*vreduced,maskvalue);CHKERRQ(ierr); 76 ierr = ISGetLocalSize(is,&nlocal);CHKERRQ(ierr); 77 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr); 78 ierr = VecGetArray(vfull,&fv);CHKERRQ(ierr); 79 ierr = VecGetArray(*vreduced,&rv);CHKERRQ(ierr); 80 ierr = ISGetIndices(is,&s);CHKERRQ(ierr); 81 if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow); 82 for (i=0;i<nlocal;i++) { 83 rv[s[i]-flow] = fv[s[i]-flow]; 84 } 85 ierr = ISRestoreIndices(is,&s);CHKERRQ(ierr); 86 ierr = VecRestoreArray(vfull,&fv);CHKERRQ(ierr); 87 ierr = VecRestoreArray(*vreduced,&rv);CHKERRQ(ierr); 88 break; 89 } 90 } 91 PetscFunctionReturn(0); 92 } 93 94 /*@C 95 TaoMatGetSubMat - Gets a submatrix using the IS 96 97 Input Parameters: 98 + M - the full matrix (n x n) 99 . is - the index set for the submatrix (both row and column index sets need to be the same) 100 . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option 101 - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, 102 TAO_SUBSET_MATRIXFREE) 103 104 Output Parameters: 105 . Msub - the submatrix 106 @*/ 107 PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub) 108 { 109 PetscErrorCode ierr; 110 IS iscomp; 111 PetscBool flg = PETSC_TRUE; 112 113 PetscFunctionBegin; 114 PetscValidHeaderSpecific(M,MAT_CLASSID,1); 115 PetscValidHeaderSpecific(is,IS_CLASSID,2); 116 ierr = MatDestroy(Msub);CHKERRQ(ierr); 117 switch (subset_type) { 118 case TAO_SUBSET_SUBVEC: 119 ierr = MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);CHKERRQ(ierr); 120 break; 121 122 case TAO_SUBSET_MASK: 123 /* Get Reduced Hessian 124 Msub[i,j] = M[i,j] if i,j in Free_Local or i==j 125 Msub[i,j] = 0 if i!=j and i or j not in Free_Local 126 */ 127 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)M),NULL,NULL,NULL);CHKERRQ(ierr); 128 ierr = PetscOptionsBool("-overwrite_hessian","modify the existing hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);CHKERRQ(ierr); 129 ierr = PetscOptionsEnd();CHKERRQ(ierr); 130 if (flg) { 131 ierr = MatDuplicate(M, MAT_COPY_VALUES, Msub);CHKERRQ(ierr); 132 } else { 133 /* Act on hessian directly (default) */ 134 ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr); 135 *Msub = M; 136 } 137 /* Save the diagonal to temporary vector */ 138 ierr = MatGetDiagonal(*Msub,v1);CHKERRQ(ierr); 139 140 /* Zero out rows and columns */ 141 ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr); 142 143 /* Use v1 instead of 0 here because of PETSc bug */ 144 ierr = MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);CHKERRQ(ierr); 145 146 ierr = ISDestroy(&iscomp);CHKERRQ(ierr); 147 break; 148 case TAO_SUBSET_MATRIXFREE: 149 ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr); 150 ierr = MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);CHKERRQ(ierr); 151 ierr = ISDestroy(&iscomp);CHKERRQ(ierr); 152 break; 153 } 154 PetscFunctionReturn(0); 155 } 156 157 /*@C 158 TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper 159 bounds, as well as fixed variables where lower and upper bounds equal each other. 160 161 Input Parameters: 162 + X - solution vector 163 . XL - lower bound vector 164 . XU - upper bound vector 165 . G - unprojected gradient 166 . S - step direction with which the active bounds will be estimated 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 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 @*/ 181 PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, PetscReal steplen, PetscReal *bound_tol, 182 IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive) 183 { 184 PetscErrorCode ierr; 185 186 Vec W; 187 PetscReal wnorm; 188 PetscInt i, n_isl=0, n_isu=0, n_isf=0; 189 PetscInt n, low, high; 190 PetscInt *isl=NULL, *isu=NULL, *isf=NULL; 191 const PetscScalar *xl, *xu, *x, *g; 192 193 PetscFunctionBegin; 194 /* Update the tolerance for bound detection (this is based on Bertsekas' method) */ 195 ierr = VecDuplicate(S, &W);CHKERRQ(ierr); 196 ierr = VecCopy(S, W);CHKERRQ(ierr); 197 ierr = VecAXPBY(W, 1.0, steplen, X);CHKERRQ(ierr); 198 ierr = VecMedian(XL, W, XU, W);CHKERRQ(ierr); 199 ierr = VecAXPBY(W, 1.0, -1.0, X);CHKERRQ(ierr); 200 ierr = VecNorm(W, NORM_2, &wnorm);CHKERRQ(ierr); 201 *bound_tol = PetscMin(*bound_tol, wnorm); 202 203 ierr = VecGetOwnershipRange(X, &low, &high);CHKERRQ(ierr); 204 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr); 205 if (n>0){ 206 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr); 207 ierr = VecGetArrayRead(XL, &xl);CHKERRQ(ierr); 208 ierr = VecGetArrayRead(XU, &xu);CHKERRQ(ierr); 209 ierr = VecGetArrayRead(G, &g);CHKERRQ(ierr); 210 211 /* Loop over variables and categorize the indexes */ 212 ierr = PetscMalloc1(n, &isl);CHKERRQ(ierr); 213 ierr = PetscMalloc1(n, &isu);CHKERRQ(ierr); 214 ierr = PetscMalloc1(n, &isf);CHKERRQ(ierr); 215 for (i=0; i<n; i++) { 216 if (xl[i] == xu[i]) { 217 /* Fixed variables here */ 218 isf[n_isf]=low+i; ++n_isf; 219 } else if ((x[i] <= xl[i] + *bound_tol) && (g[i] > 0.0)) { 220 /* Lower bounded variables here */ 221 isl[n_isl]=low+i; ++n_isl; 222 } else if ((x[i] >= xu[i] - *bound_tol) && (g[i] < 0.0)) { 223 /* Upper bounded variables here */ 224 isu[n_isu]=low+i; ++n_isu; 225 } 226 } 227 228 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr); 229 ierr = VecRestoreArrayRead(XL, &xl);CHKERRQ(ierr); 230 ierr = VecRestoreArrayRead(XU, &xu);CHKERRQ(ierr); 231 ierr = VecRestoreArrayRead(G, &g);CHKERRQ(ierr); 232 } 233 234 /* Create index set for lower bounded variables */ 235 ierr = ISDestroy(active_lower);CHKERRQ(ierr); 236 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)X), n_isl, isl, PETSC_OWN_POINTER, active_lower);CHKERRQ(ierr); 237 /* Create index set for upper bounded variables */ 238 ierr = ISDestroy(active_upper);CHKERRQ(ierr); 239 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)X), n_isu, isu, PETSC_OWN_POINTER, active_upper);CHKERRQ(ierr); 240 /* Create index set for fixed variables */ 241 ierr = ISDestroy(active_fixed);CHKERRQ(ierr); 242 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)X), n_isf, isf, PETSC_OWN_POINTER, active_fixed);CHKERRQ(ierr); 243 244 /* Create the combined active set */ 245 ierr = ISDestroy(active);CHKERRQ(ierr); 246 if (*active_lower && *active_upper && *active_fixed) { 247 /* All three types of active variables exist */ 248 const IS islist[3] = {*active_lower, *active_upper, *active_fixed}; 249 ierr = ISConcatenate(PetscObjectComm((PetscObject)X), 3, islist, active);CHKERRQ(ierr); 250 ierr = ISSort(*active);CHKERRQ(ierr); 251 } else if (*active_lower && *active_upper) { 252 /* Only lower and upper bounded active variables exist */ 253 ierr = ISSum(*active_lower, *active_upper, active);CHKERRQ(ierr); 254 } else if (*active_lower && *active_fixed) { 255 /* Only lower bounded and fixed active variables exist */ 256 ierr = ISSum(*active_lower, *active_fixed, active);CHKERRQ(ierr); 257 } else if (*active_upper && *active_fixed) { 258 /* Only upper bounded and fixed active variables exist */ 259 ierr = ISSum(*active_upper, *active_fixed, active);CHKERRQ(ierr); 260 } else if (*active_lower) { 261 /* Only lower bounded active variables exist */ 262 *active = *active_lower; 263 } else if (*active_upper) { 264 /* Only upper bounded active variables exist */ 265 *active = *active_upper; 266 } else if (*active_fixed) { 267 /* Only fixed active variables exist */ 268 *active = *active_fixed; 269 } 270 /* Create the inactive set */ 271 ierr = ISDestroy(inactive);CHKERRQ(ierr); 272 if (*active) { ierr = ISComplementVec(*active, X, inactive);CHKERRQ(ierr); } 273 274 /* Clean up and exit */ 275 ierr = VecDestroy(&W);CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 /*@C 280 TaoBoundStep - Ensures the correct zero or adjusted step direction 281 values for active variables. 282 283 Input Parameters: 284 + X - solution vector 285 . XL - lower bound vector 286 . XU - upper bound vector 287 . active_lower - index set for lower bounded active variables 288 . active_upper - index set for lower bounded active variables 289 - active_fixed - index set for fixed active variables 290 291 Output Parameters: 292 . S - step direction to be modified 293 @*/ 294 PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, Vec S) 295 { 296 PetscErrorCode ierr; 297 298 Vec step_lower, step_upper, step_fixed; 299 Vec x_lower, x_upper; 300 Vec bound_lower, bound_upper; 301 302 PetscFunctionBegin; 303 /* Adjust step for variables at the estimated lower bound */ 304 if (active_lower) { 305 ierr = VecGetSubVector(S, active_lower, &step_lower);CHKERRQ(ierr); 306 ierr = VecGetSubVector(X, active_lower, &x_lower);CHKERRQ(ierr); 307 ierr = VecGetSubVector(XL, active_lower, &bound_lower);CHKERRQ(ierr); 308 ierr = VecCopy(bound_lower, step_lower);CHKERRQ(ierr); 309 ierr = VecAXPY(step_lower, -1.0, x_lower);CHKERRQ(ierr); 310 ierr = VecRestoreSubVector(S, active_lower, &step_lower);CHKERRQ(ierr); 311 ierr = VecRestoreSubVector(X, active_lower, &x_lower);CHKERRQ(ierr); 312 ierr = VecRestoreSubVector(XL, active_lower, &bound_lower);CHKERRQ(ierr); 313 } 314 315 /* Adjust step for the variables at the estimated upper bound */ 316 if (active_upper) { 317 ierr = VecGetSubVector(S, active_upper, &step_upper);CHKERRQ(ierr); 318 ierr = VecGetSubVector(X, active_upper, &x_upper);CHKERRQ(ierr); 319 ierr = VecGetSubVector(XU, active_upper, &bound_upper);CHKERRQ(ierr); 320 ierr = VecCopy(bound_upper, step_upper);CHKERRQ(ierr); 321 ierr = VecAXPY(step_upper, -1.0, x_upper);CHKERRQ(ierr); 322 ierr = VecRestoreSubVector(S, active_upper, &step_upper);CHKERRQ(ierr); 323 ierr = VecRestoreSubVector(X, active_upper, &x_upper);CHKERRQ(ierr); 324 ierr = VecRestoreSubVector(XU, active_upper, &bound_upper);CHKERRQ(ierr); 325 } 326 327 /* Zero out step for fixed variables */ 328 if (active_fixed) { 329 ierr = VecGetSubVector(S, active_fixed, &step_fixed);CHKERRQ(ierr); 330 ierr = VecSet(step_fixed, 0.0);CHKERRQ(ierr); 331 ierr = VecRestoreSubVector(S, active_fixed, &step_fixed);CHKERRQ(ierr); 332 } 333 PetscFunctionReturn(0); 334 } 335