Actual source code: isutil.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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 Parameters:
 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: {
 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;


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

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

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

 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:         VecDuplicate(vfull,vreduced);
 75:       }

 77:       VecSet(*vreduced,maskvalue);
 78:       ISGetLocalSize(is,&nlocal);
 79:       VecGetOwnershipRange(vfull,&flow,&fhigh);
 80:       VecGetArray(vfull,&fv);
 81:       VecGetArray(*vreduced,&rv);
 82:       ISGetIndices(is,&s);
 83:       if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_SELF,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:       ISRestoreIndices(is,&s);
 88:       VecRestoreArray(vfull,&fv);
 89:       VecRestoreArray(*vreduced,&rv);
 90:       break;
 91:     }
 92:   }
 93:   return(0);
 94: }

 96: /*@C
 97:   TaoMatGetSubMat - Gets a submatrix using the IS

 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

105:   Output Parameters:
106: . Msub - the submatrix

108:   Level: developer
109: @*/
110: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
111: {
113:   IS             iscomp;
114:   PetscBool      flg = PETSC_TRUE;

119:   MatDestroy(Msub);
120:   switch (subset_type) {
121:   case TAO_SUBSET_SUBVEC:
122:     MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
123:     break;

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:     PetscObjectOptionsBegin((PetscObject)M);
131:     PetscOptionsBool("-overwrite_hessian","modify the existing hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);
132:     PetscOptionsEnd();
133:     if (flg) {
134:       MatDuplicate(M, MAT_COPY_VALUES, Msub);
135:     } else {
136:       /* Act on hessian directly (default) */
137:       PetscObjectReference((PetscObject)M);
138:       *Msub = M;
139:     }
140:     /* Save the diagonal to temporary vector */
141:     MatGetDiagonal(*Msub,v1);

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

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

149:     ISDestroy(&iscomp);
150:     break;
151:   case TAO_SUBSET_MATRIXFREE:
152:     ISComplementVec(is,v1,&iscomp);
153:     MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
154:     ISDestroy(&iscomp);
155:     break;
156:   }
157:   return(0);
158: }

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.

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)

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

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

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);


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);

229:   /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
230:   VecCopy(X, W);
231:   VecAXPBY(W, steplen, 1.0, S);
232:   TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W);
233:   VecAXPBY(W, 1.0, -1.0, X);
234:   VecNorm(W, NORM_2, &wnorm);
235:   *bound_tol = PetscMin(*bound_tol, wnorm);

237:   VecGetOwnershipRange(X, &low, &high);
238:   VecGetLocalSize(X, &n);
239:   if (n>0){
240:     VecGetArrayRead(X, &x);
241:     VecGetArrayRead(XL, &xl);
242:     VecGetArrayRead(XU, &xu);
243:     VecGetArrayRead(G, &g);

245:     /* Loop over variables and categorize the indexes */
246:     PetscMalloc1(n, &isl);
247:     PetscMalloc1(n, &isu);
248:     PetscMalloc1(n, &isf);
249:     PetscMalloc1(n, &isa);
250:     PetscMalloc1(n, &isi);
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:     }

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

276:   /* Clear all index sets */
277:   ISDestroy(active_lower);
278:   ISDestroy(active_upper);
279:   ISDestroy(active_fixed);
280:   ISDestroy(active);
281:   ISDestroy(inactive);

283:   /* Collect global sizes */
284:   MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);
285:   MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);
286:   MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);
287:   MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);
288:   MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm);

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

321:   /* Clean up and exit */
322:   return(0);
323: }

325: /*@C
326:   TaoBoundStep - Ensures the correct zero or adjusted step direction
327:   values for active variables.

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

338:   Output Parameters:
339: . S - step direction to be modified

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;

347:   Vec                          step_lower, step_upper, step_fixed;
348:   Vec                          x_lower, x_upper;
349:   Vec                          bound_lower, bound_upper;

352:   /* Adjust step for variables at the estimated lower bound */
353:   if (active_lower) {
354:     VecGetSubVector(S, active_lower, &step_lower);
355:     VecGetSubVector(X, active_lower, &x_lower);
356:     VecGetSubVector(XL, active_lower, &bound_lower);
357:     VecCopy(bound_lower, step_lower);
358:     VecAXPY(step_lower, -1.0, x_lower);
359:     VecScale(step_lower, scale);
360:     VecRestoreSubVector(S, active_lower, &step_lower);
361:     VecRestoreSubVector(X, active_lower, &x_lower);
362:     VecRestoreSubVector(XL, active_lower, &bound_lower);
363:   }

365:   /* Adjust step for the variables at the estimated upper bound */
366:   if (active_upper) {
367:     VecGetSubVector(S, active_upper, &step_upper);
368:     VecGetSubVector(X, active_upper, &x_upper);
369:     VecGetSubVector(XU, active_upper, &bound_upper);
370:     VecCopy(bound_upper, step_upper);
371:     VecAXPY(step_upper, -1.0, x_upper);
372:     VecScale(step_upper, scale);
373:     VecRestoreSubVector(S, active_upper, &step_upper);
374:     VecRestoreSubVector(X, active_upper, &x_upper);
375:     VecRestoreSubVector(XU, active_upper, &bound_upper);
376:   }

378:   /* Zero out step for fixed variables */
379:   if (active_fixed) {
380:     VecGetSubVector(S, active_fixed, &step_fixed);
381:     VecSet(step_fixed, 0.0);
382:     VecRestoreSubVector(S, active_fixed, &step_fixed);
383:   }
384:   return(0);
385: }

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

390:   Collective on Vec

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

398:   Output Parameters:
399: + nDiff - total number of vector entries that have been bounded
400: - Xout - modified solution vector satisfying bounds to bound_tol

402:   Level: developer

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;


429:   VecCheckSameSize(X,1,XL,2);
430:   VecCheckSameSize(X,1,XU,3);
431:   VecCheckSameSize(X,1,Xout,4);

433:   VecGetOwnershipRange(X,&low,&high);
434:   VecGetLocalSize(X,&n);
435:   if (n>0){
436:     VecGetArrayRead(X, &x);
437:     VecGetArrayRead(XL, &xl);
438:     VecGetArrayRead(XU, &xu);
439:     VecGetArray(Xout, &xout);

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:     }

449:     VecRestoreArrayRead(X, &x);
450:     VecRestoreArrayRead(XL, &xl);
451:     VecRestoreArrayRead(XU, &xu);
452:     VecRestoreArray(Xout, &xout);
453:   }
454:   MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X));
455:   return(0);
456: }