Actual source code: isutil.c

  1: #include <petsctao.h>
  2: #include <petsc/private/vecimpl.h>
  3: #include <petsc/private/taoimpl.h>
  4: #include <../src/tao/matrix/submatfree.h>

  6: /*@C
  7:   TaoVecGetSubVec - Gets a subvector using the IS

  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)

 15:   Output Parameter:
 16: . vreduced - the subvector

 18:   Notes:
 19:   maskvalue should usually be 0.0, unless a pointwise divide will be used.

 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;


 37:   VecGetSize(vfull, &nfull);
 38:   ISGetSize(is, &nreduced);

 40:   if (nreduced == nfull) {
 41:     VecDestroy(vreduced);
 42:     VecDuplicate(vfull, vreduced);
 43:     VecCopy(vfull, *vreduced);
 44:   } else {
 45:     switch (reduced_type) {
 46:     case TAO_SUBSET_SUBVEC:
 47:       VecGetType(vfull, &vtype);
 48:       VecGetOwnershipRange(vfull, &flow, &fhigh);
 49:       ISGetLocalSize(is, &nreduced_local);
 50:       PetscObjectGetComm((PetscObject)vfull, &comm);
 51:       if (*vreduced) VecDestroy(vreduced);
 52:       VecCreate(comm, vreduced);
 53:       VecSetType(*vreduced, vtype);

 55:       VecSetSizes(*vreduced, nreduced_local, nreduced);
 56:       VecGetOwnershipRange(*vreduced, &rlow, &rhigh);
 57:       ISCreateStride(comm, nreduced_local, rlow, 1, &ident);
 58:       VecScatterCreate(vfull, is, *vreduced, ident, &scatter);
 59:       VecScatterBegin(scatter, vfull, *vreduced, INSERT_VALUES, SCATTER_FORWARD);
 60:       VecScatterEnd(scatter, vfull, *vreduced, INSERT_VALUES, SCATTER_FORWARD);
 61:       VecScatterDestroy(&scatter);
 62:       ISDestroy(&ident);
 63:       break;

 65:     case TAO_SUBSET_MASK:
 66:     case TAO_SUBSET_MATRIXFREE:
 67:       /* vr[i] = vf[i]   if i in is
 68:        vr[i] = 0       otherwise */
 69:       if (!*vreduced) VecDuplicate(vfull, vreduced);

 71:       VecSet(*vreduced, maskvalue);
 72:       ISGetLocalSize(is, &nlocal);
 73:       VecGetOwnershipRange(vfull, &flow, &fhigh);
 74:       VecGetArray(vfull, &fv);
 75:       VecGetArray(*vreduced, &rv);
 76:       ISGetIndices(is, &s);
 78:       for (i = 0; i < nlocal; ++i) rv[s[i] - flow] = fv[s[i] - flow];
 79:       ISRestoreIndices(is, &s);
 80:       VecRestoreArray(vfull, &fv);
 81:       VecRestoreArray(*vreduced, &rv);
 82:       break;
 83:     }
 84:   }
 85:   return 0;
 86: }

 88: /*@C
 89:   TaoMatGetSubMat - Gets a submatrix using the IS

 91:   Input Parameters:
 92: + M - the full matrix (n x n)
 93: . is - the index set for the submatrix (both row and column index sets need to be the same)
 94: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
 95: - subset_type <TAO_SUBSET_SUBVEC,TAO_SUBSET_MASK,TAO_SUBSET_MATRIXFREE> - the method TAO is using for subsetting

 97:   Output Parameter:
 98: . Msub - the submatrix

100:   Level: developer
101: @*/
102: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
103: {
104:   IS        iscomp;
105:   PetscBool flg = PETSC_TRUE;

109:   MatDestroy(Msub);
110:   switch (subset_type) {
111:   case TAO_SUBSET_SUBVEC:
112:     MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
113:     break;

115:   case TAO_SUBSET_MASK:
116:     /* Get Reduced Hessian
117:      Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
118:      Msub[i,j] = 0      if i!=j and i or j not in Free_Local
119:      */
120:     PetscObjectOptionsBegin((PetscObject)M);
121:     PetscOptionsBool("-overwrite_hessian", "modify the existing hessian matrix when computing submatrices", "TaoSubsetType", flg, &flg, NULL);
122:     PetscOptionsEnd();
123:     if (flg) {
124:       MatDuplicate(M, MAT_COPY_VALUES, Msub);
125:     } else {
126:       /* Act on hessian directly (default) */
127:       PetscObjectReference((PetscObject)M);
128:       *Msub = M;
129:     }
130:     /* Save the diagonal to temporary vector */
131:     MatGetDiagonal(*Msub, v1);

133:     /* Zero out rows and columns */
134:     ISComplementVec(is, v1, &iscomp);

136:     /* Use v1 instead of 0 here because of PETSc bug */
137:     MatZeroRowsColumnsIS(*Msub, iscomp, 1.0, v1, v1);

139:     ISDestroy(&iscomp);
140:     break;
141:   case TAO_SUBSET_MATRIXFREE:
142:     ISComplementVec(is, v1, &iscomp);
143:     MatCreateSubMatrixFree(M, iscomp, iscomp, Msub);
144:     ISDestroy(&iscomp);
145:     break;
146:   }
147:   return 0;
148: }

150: /*@C
151:   TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper
152:   bounds, as well as fixed variables where lower and upper bounds equal each other.

154:   Input Parameters:
155: + X - solution vector
156: . XL - lower bound vector
157: . XU - upper bound vector
158: . G - unprojected gradient
159: . S - step direction with which the active bounds will be estimated
160: . W - work vector of type and size of X
161: - steplen - the step length at which the active bounds will be estimated (needs to be conservative)

163:   Output Parameters:
164: + bound_tol - tolerance for the bound estimation
165: . active_lower - index set for active variables at the lower bound
166: . active_upper - index set for active variables at the upper bound
167: . active_fixed - index set for fixed variables
168: . active - index set for all active variables
169: - inactive - complementary index set for inactive variables

171:   Notes:
172:   This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3.

174:   Level: developer
175: @*/
176: PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol, IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive)
177: {
178:   PetscReal          wnorm;
179:   PetscReal          zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
180:   PetscInt           i, n_isl = 0, n_isu = 0, n_isf = 0, n_isa = 0, n_isi = 0;
181:   PetscInt           N_isl, N_isu, N_isf, N_isa, N_isi;
182:   PetscInt           n, low, high, nDiff;
183:   PetscInt          *isl = NULL, *isu = NULL, *isf = NULL, *isa = NULL, *isi = NULL;
184:   const PetscScalar *xl, *xu, *x, *g;
185:   MPI_Comm           comm = PetscObjectComm((PetscObject)X);


204:   if (XL) VecCheckSameSize(X, 1, XL, 2);
205:   if (XU) VecCheckSameSize(X, 1, XU, 3);
206:   VecCheckSameSize(X, 1, G, 4);
207:   VecCheckSameSize(X, 1, S, 5);
208:   VecCheckSameSize(X, 1, W, 6);

210:   /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
211:   VecCopy(X, W);
212:   VecAXPBY(W, steplen, 1.0, S);
213:   TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W);
214:   VecAXPBY(W, 1.0, -1.0, X);
215:   VecNorm(W, NORM_2, &wnorm);
216:   *bound_tol = PetscMin(*bound_tol, wnorm);

218:   /* Clear all index sets */
219:   ISDestroy(active_lower);
220:   ISDestroy(active_upper);
221:   ISDestroy(active_fixed);
222:   ISDestroy(active);
223:   ISDestroy(inactive);

225:   VecGetOwnershipRange(X, &low, &high);
226:   VecGetLocalSize(X, &n);
227:   if (!XL && !XU) {
228:     ISCreateStride(comm, n, low, 1, inactive);
229:     return 0;
230:   }
231:   if (n > 0) {
232:     VecGetArrayRead(X, &x);
233:     VecGetArrayRead(XL, &xl);
234:     VecGetArrayRead(XU, &xu);
235:     VecGetArrayRead(G, &g);

237:     /* Loop over variables and categorize the indexes */
238:     PetscMalloc1(n, &isl);
239:     PetscMalloc1(n, &isu);
240:     PetscMalloc1(n, &isf);
241:     PetscMalloc1(n, &isa);
242:     PetscMalloc1(n, &isi);
243:     for (i = 0; i < n; ++i) {
244:       if (xl[i] == xu[i]) {
245:         /* Fixed variables */
246:         isf[n_isf] = low + i;
247:         ++n_isf;
248:         isa[n_isa] = low + i;
249:         ++n_isa;
250:       } else if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + *bound_tol && g[i] > zero) {
251:         /* Lower bounded variables */
252:         isl[n_isl] = low + i;
253:         ++n_isl;
254:         isa[n_isa] = low + i;
255:         ++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;
259:         ++n_isu;
260:         isa[n_isa] = low + i;
261:         ++n_isa;
262:       } else {
263:         /* Inactive variables */
264:         isi[n_isi] = low + i;
265:         ++n_isi;
266:       }
267:     }

269:     VecRestoreArrayRead(X, &x);
270:     VecRestoreArrayRead(XL, &xl);
271:     VecRestoreArrayRead(XU, &xu);
272:     VecRestoreArrayRead(G, &g);
273:   }

275:   /* Collect global sizes */
276:   MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);
277:   MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);
278:   MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);
279:   MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);
280:   MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm);

282:   /* Create index set for lower bounded variables */
283:   if (N_isl > 0) {
284:     ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower);
285:   } else {
286:     PetscFree(isl);
287:   }
288:   /* Create index set for upper bounded variables */
289:   if (N_isu > 0) {
290:     ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper);
291:   } else {
292:     PetscFree(isu);
293:   }
294:   /* Create index set for fixed variables */
295:   if (N_isf > 0) {
296:     ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed);
297:   } else {
298:     PetscFree(isf);
299:   }
300:   /* Create index set for all actively bounded variables */
301:   if (N_isa > 0) {
302:     ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active);
303:   } else {
304:     PetscFree(isa);
305:   }
306:   /* Create index set for all inactive variables */
307:   if (N_isi > 0) {
308:     ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive);
309:   } else {
310:     PetscFree(isi);
311:   }
312:   return 0;
313: }

315: /*@C
316:   TaoBoundStep - Ensures the correct zero or adjusted step direction
317:   values for active variables.

319:   Input Parameters:
320: + X - solution vector
321: . XL - lower bound vector
322: . XU - upper bound vector
323: . active_lower - index set for lower bounded active variables
324: . active_upper - index set for lower bounded active variables
325: . active_fixed - index set for fixed active variables
326: - scale - amplification factor for the step that needs to be taken on actively bounded variables

328:   Output Parameter:
329: . S - step direction to be modified

331:   Level: developer
332: @*/
333: PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S)
334: {
335:   Vec step_lower, step_upper, step_fixed;
336:   Vec x_lower, x_upper;
337:   Vec bound_lower, bound_upper;

339:   /* Adjust step for variables at the estimated lower bound */
340:   if (active_lower) {
341:     VecGetSubVector(S, active_lower, &step_lower);
342:     VecGetSubVector(X, active_lower, &x_lower);
343:     VecGetSubVector(XL, active_lower, &bound_lower);
344:     VecCopy(bound_lower, step_lower);
345:     VecAXPY(step_lower, -1.0, x_lower);
346:     VecScale(step_lower, scale);
347:     VecRestoreSubVector(S, active_lower, &step_lower);
348:     VecRestoreSubVector(X, active_lower, &x_lower);
349:     VecRestoreSubVector(XL, active_lower, &bound_lower);
350:   }

352:   /* Adjust step for the variables at the estimated upper bound */
353:   if (active_upper) {
354:     VecGetSubVector(S, active_upper, &step_upper);
355:     VecGetSubVector(X, active_upper, &x_upper);
356:     VecGetSubVector(XU, active_upper, &bound_upper);
357:     VecCopy(bound_upper, step_upper);
358:     VecAXPY(step_upper, -1.0, x_upper);
359:     VecScale(step_upper, scale);
360:     VecRestoreSubVector(S, active_upper, &step_upper);
361:     VecRestoreSubVector(X, active_upper, &x_upper);
362:     VecRestoreSubVector(XU, active_upper, &bound_upper);
363:   }

365:   /* Zero out step for fixed variables */
366:   if (active_fixed) {
367:     VecGetSubVector(S, active_fixed, &step_fixed);
368:     VecSet(step_fixed, 0.0);
369:     VecRestoreSubVector(S, active_fixed, &step_fixed);
370:   }
371:   return 0;
372: }

374: /*@C
375:   TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance.

377:   Collective

379:   Input Parameters:
380: + X - solution vector
381: . XL - lower bound vector
382: . XU - upper bound vector
383: - bound_tol - absolute tolerance in enforcing the bound

385:   Output Parameters:
386: + nDiff - total number of vector entries that have been bounded
387: - Xout - modified solution vector satisfying bounds to bound_tol

389:   Level: developer

391: .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`
392: @*/
393: PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout)
394: {
395:   PetscInt           i, n, low, high, nDiff_loc = 0;
396:   PetscScalar       *xout;
397:   const PetscScalar *x, *xl, *xu;

403:   if (!XL && !XU) {
404:     VecCopy(X, Xout);
405:     *nDiff = 0.0;
406:     return 0;
407:   }
414:   VecCheckSameSize(X, 1, XL, 2);
415:   VecCheckSameSize(X, 1, XU, 3);
416:   VecCheckSameSize(X, 1, Xout, 4);

418:   VecGetOwnershipRange(X, &low, &high);
419:   VecGetLocalSize(X, &n);
420:   if (n > 0) {
421:     VecGetArrayRead(X, &x);
422:     VecGetArrayRead(XL, &xl);
423:     VecGetArrayRead(XU, &xu);
424:     VecGetArray(Xout, &xout);

426:     for (i = 0; i < n; ++i) {
427:       if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + bound_tol) {
428:         xout[i] = xl[i];
429:         ++nDiff_loc;
430:       } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - bound_tol) {
431:         xout[i] = xu[i];
432:         ++nDiff_loc;
433:       }
434:     }

436:     VecRestoreArrayRead(X, &x);
437:     VecRestoreArrayRead(XL, &xl);
438:     VecRestoreArrayRead(XU, &xu);
439:     VecRestoreArray(Xout, &xout);
440:   }
441:   MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X));
442:   return 0;
443: }