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;

 34:   PetscFunctionBegin;

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

 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) PetscCall(VecDestroy(vreduced));
 53:       PetscCall(VecCreate(comm, vreduced));
 54:       PetscCall(VecSetType(*vreduced, vtype));

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

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

 72:       PetscCall(VecSet(*vreduced, maskvalue));
 73:       PetscCall(ISGetLocalSize(is, &nlocal));
 74:       PetscCall(VecGetOwnershipRange(vfull, &flow, &fhigh));
 75:       PetscCall(VecGetArray(vfull, &fv));
 76:       PetscCall(VecGetArray(*vreduced, &rv));
 77:       PetscCall(ISGetIndices(is, &s));
 78:       PetscCheck(nlocal <= (fhigh - flow), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "IS local size %" PetscInt_FMT " > Vec local size %" PetscInt_FMT, nlocal, fhigh - flow);
 79:       for (i = 0; i < nlocal; ++i) rv[s[i] - flow] = fv[s[i] - flow];
 80:       PetscCall(ISRestoreIndices(is, &s));
 81:       PetscCall(VecRestoreArray(vfull, &fv));
 82:       PetscCall(VecRestoreArray(*vreduced, &rv));
 83:       break;
 84:     }
 85:   }
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

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

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

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

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

108:   PetscFunctionBegin;
111:   PetscCall(MatDestroy(Msub));
112:   switch (subset_type) {
113:   case TAO_SUBSET_SUBVEC:
114:     PetscCall(MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub));
115:     break;

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

135:     /* Zero out rows and columns */
136:     PetscCall(ISComplementVec(is, v1, &iscomp));

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

141:     PetscCall(ISDestroy(&iscomp));
142:     break;
143:   case TAO_SUBSET_MATRIXFREE:
144:     PetscCall(ISComplementVec(is, v1, &iscomp));
145:     PetscCall(MatCreateSubMatrixFree(M, iscomp, iscomp, Msub));
146:     PetscCall(ISDestroy(&iscomp));
147:     break;
148:   }
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

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

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

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

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

176:   Level: developer
177: @*/
178: 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)
179: {
180:   PetscReal          wnorm;
181:   PetscReal          zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
182:   PetscInt           i, n_isl = 0, n_isu = 0, n_isf = 0, n_isa = 0, n_isi = 0;
183:   PetscInt           N_isl, N_isu, N_isf, N_isa, N_isi;
184:   PetscInt           n, low, high, nDiff;
185:   PetscInt          *isl = NULL, *isu = NULL, *isf = NULL, *isa = NULL, *isi = NULL;
186:   const PetscScalar *xl, *xu, *x, *g;
187:   MPI_Comm           comm = PetscObjectComm((PetscObject)X);

189:   PetscFunctionBegin;

197:   if (XL) PetscCheckSameType(X, 1, XL, 2);
198:   if (XU) PetscCheckSameType(X, 1, XU, 3);
199:   PetscCheckSameType(X, 1, G, 4);
200:   PetscCheckSameType(X, 1, S, 5);
201:   PetscCheckSameType(X, 1, W, 6);
202:   if (XL) PetscCheckSameComm(X, 1, XL, 2);
203:   if (XU) PetscCheckSameComm(X, 1, XU, 3);
204:   PetscCheckSameComm(X, 1, G, 4);
205:   PetscCheckSameComm(X, 1, S, 5);
206:   PetscCheckSameComm(X, 1, W, 6);
207:   if (XL) VecCheckSameSize(X, 1, XL, 2);
208:   if (XU) VecCheckSameSize(X, 1, XU, 3);
209:   VecCheckSameSize(X, 1, G, 4);
210:   VecCheckSameSize(X, 1, S, 5);
211:   VecCheckSameSize(X, 1, W, 6);

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

221:   /* Clear all index sets */
222:   PetscCall(ISDestroy(active_lower));
223:   PetscCall(ISDestroy(active_upper));
224:   PetscCall(ISDestroy(active_fixed));
225:   PetscCall(ISDestroy(active));
226:   PetscCall(ISDestroy(inactive));

228:   PetscCall(VecGetOwnershipRange(X, &low, &high));
229:   PetscCall(VecGetLocalSize(X, &n));
230:   if (!XL && !XU) {
231:     PetscCall(ISCreateStride(comm, n, low, 1, inactive));
232:     PetscFunctionReturn(PETSC_SUCCESS);
233:   }
234:   if (n > 0) {
235:     PetscCall(VecGetArrayRead(X, &x));
236:     PetscCall(VecGetArrayRead(XL, &xl));
237:     PetscCall(VecGetArrayRead(XU, &xu));
238:     PetscCall(VecGetArrayRead(G, &g));

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

272:     PetscCall(VecRestoreArrayRead(X, &x));
273:     PetscCall(VecRestoreArrayRead(XL, &xl));
274:     PetscCall(VecRestoreArrayRead(XU, &xu));
275:     PetscCall(VecRestoreArrayRead(G, &g));
276:   }

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

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

318: /*@C
319:   TaoBoundStep - Ensures the correct zero or adjusted step direction
320:   values for active variables.

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

331:   Output Parameter:
332: . S - step direction to be modified

334:   Level: developer
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:   Vec step_lower, step_upper, step_fixed;
339:   Vec x_lower, x_upper;
340:   Vec bound_lower, bound_upper;

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

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

369:   /* Zero out step for fixed variables */
370:   if (active_fixed) {
371:     PetscCall(VecGetSubVector(S, active_fixed, &step_fixed));
372:     PetscCall(VecSet(step_fixed, 0.0));
373:     PetscCall(VecRestoreSubVector(S, active_fixed, &step_fixed));
374:   }
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

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

381:   Collective

383:   Input Parameters:
384: + X - solution vector
385: . XL - lower bound vector
386: . XU - upper bound vector
387: - bound_tol - absolute tolerance in enforcing the bound

389:   Output Parameters:
390: + nDiff - total number of vector entries that have been bounded
391: - Xout - modified solution vector satisfying bounds to bound_tol

393:   Level: developer

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

403:   PetscFunctionBegin;
408:   if (!XL && !XU) {
409:     PetscCall(VecCopy(X, Xout));
410:     *nDiff = 0.0;
411:     PetscFunctionReturn(PETSC_SUCCESS);
412:   }
413:   PetscCheckSameType(X, 1, XL, 2);
414:   PetscCheckSameType(X, 1, XU, 3);
415:   PetscCheckSameType(X, 1, Xout, 6);
416:   PetscCheckSameComm(X, 1, XL, 2);
417:   PetscCheckSameComm(X, 1, XU, 3);
418:   PetscCheckSameComm(X, 1, Xout, 6);
419:   VecCheckSameSize(X, 1, XL, 2);
420:   VecCheckSameSize(X, 1, XU, 3);
421:   VecCheckSameSize(X, 1, Xout, 4);

423:   PetscCall(VecGetOwnershipRange(X, &low, &high));
424:   PetscCall(VecGetLocalSize(X, &n));
425:   if (n > 0) {
426:     PetscCall(VecGetArrayRead(X, &x));
427:     PetscCall(VecGetArrayRead(XL, &xl));
428:     PetscCall(VecGetArrayRead(XU, &xu));
429:     PetscCall(VecGetArray(Xout, &xout));

431:     for (i = 0; i < n; ++i) {
432:       if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + bound_tol) {
433:         xout[i] = xl[i];
434:         ++nDiff_loc;
435:       } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - bound_tol) {
436:         xout[i] = xu[i];
437:         ++nDiff_loc;
438:       }
439:     }

441:     PetscCall(VecRestoreArrayRead(X, &x));
442:     PetscCall(VecRestoreArrayRead(XL, &xl));
443:     PetscCall(VecRestoreArrayRead(XU, &xu));
444:     PetscCall(VecRestoreArray(Xout, &xout));
445:   }
446:   PetscCall(MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X)));
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }