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

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

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

 75:       VecSet(*vreduced,maskvalue);
 76:       ISGetLocalSize(is,&nlocal);
 77:       VecGetOwnershipRange(vfull,&flow,&fhigh);
 78:       VecGetArray(vfull,&fv);
 79:       VecGetArray(*vreduced,&rv);
 80:       ISGetIndices(is,&s);
 82:       for (i=0;i<nlocal;++i) {
 83:         rv[s[i]-flow] = fv[s[i]-flow];
 84:       }
 85:       ISRestoreIndices(is,&s);
 86:       VecRestoreArray(vfull,&fv);
 87:       VecRestoreArray(*vreduced,&rv);
 88:       break;
 89:     }
 90:   }
 91:   return 0;
 92: }

 94: /*@C
 95:   TaoMatGetSubMat - Gets a submatrix using the IS

 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 <TAO_SUBSET_SUBVEC,TAO_SUBSET_MASK,TAO_SUBSET_MATRIXFREE> - the method TAO is using for subsetting

103:   Output Parameter:
104: . Msub - the submatrix

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

116:   MatDestroy(Msub);
117:   switch (subset_type) {
118:   case TAO_SUBSET_SUBVEC:
119:     MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
120:     break;

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

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

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

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

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.

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: . W - work vector of type and size of X
168: - steplen - the step length at which the active bounds will be estimated (needs to be conservative)

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

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

181:   Level: developer
182: @*/
183: PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol,
184:                                        IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive)
185: {
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);


218:   VecCheckSameSize(X,1,XL,2);
219:   VecCheckSameSize(X,1,XU,3);
220:   VecCheckSameSize(X,1,G,4);
221:   VecCheckSameSize(X,1,S,5);
222:   VecCheckSameSize(X,1,W,6);

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

232:   VecGetOwnershipRange(X, &low, &high);
233:   VecGetLocalSize(X, &n);
234:   if (n>0) {
235:     VecGetArrayRead(X, &x);
236:     VecGetArrayRead(XL, &xl);
237:     VecGetArrayRead(XU, &xu);
238:     VecGetArrayRead(G, &g);

240:     /* Loop over variables and categorize the indexes */
241:     PetscMalloc1(n, &isl);
242:     PetscMalloc1(n, &isu);
243:     PetscMalloc1(n, &isf);
244:     PetscMalloc1(n, &isa);
245:     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; ++n_isf;
250:         isa[n_isa]=low+i; ++n_isa;
251:       } else if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + *bound_tol) && (g[i] > zero)) {
252:         /* Lower bounded variables */
253:         isl[n_isl]=low+i; ++n_isl;
254:         isa[n_isa]=low+i; ++n_isa;
255:       } else if ((xu[i] < PETSC_INFINITY) && (x[i] >= xu[i] - *bound_tol) && (g[i] < zero)) {
256:         /* Upper bounded variables */
257:         isu[n_isu]=low+i; ++n_isu;
258:         isa[n_isa]=low+i; ++n_isa;
259:       } else {
260:         /* Inactive variables */
261:         isi[n_isi]=low+i; ++n_isi;
262:       }
263:     }

265:     VecRestoreArrayRead(X, &x);
266:     VecRestoreArrayRead(XL, &xl);
267:     VecRestoreArrayRead(XU, &xu);
268:     VecRestoreArrayRead(G, &g);
269:   }

271:   /* Clear all index sets */
272:   ISDestroy(active_lower);
273:   ISDestroy(active_upper);
274:   ISDestroy(active_fixed);
275:   ISDestroy(active);
276:   ISDestroy(inactive);

278:   /* Collect global sizes */
279:   MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);
280:   MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);
281:   MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);
282:   MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);
283:   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:     ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower);
288:   } else {
289:     PetscFree(isl);
290:   }
291:   /* Create index set for upper bounded variables */
292:   if (N_isu > 0) {
293:     ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper);
294:   } else {
295:     PetscFree(isu);
296:   }
297:   /* Create index set for fixed variables */
298:   if (N_isf > 0) {
299:     ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed);
300:   } else {
301:     PetscFree(isf);
302:   }
303:   /* Create index set for all actively bounded variables */
304:   if (N_isa > 0) {
305:     ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active);
306:   } else {
307:     PetscFree(isa);
308:   }
309:   /* Create index set for all inactive variables */
310:   if (N_isi > 0) {
311:     ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive);
312:   } else {
313:     PetscFree(isi);
314:   }

316:   /* Clean up and exit */
317:   return 0;
318: }

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

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

333:   Output Parameter:
334: . S - step direction to be modified

336:   Level: developer
337: @*/
338: PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S)
339: {

341:   Vec                          step_lower, step_upper, step_fixed;
342:   Vec                          x_lower, x_upper;
343:   Vec                          bound_lower, bound_upper;

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

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

371:   /* Zero out step for fixed variables */
372:   if (active_fixed) {
373:     VecGetSubVector(S, active_fixed, &step_fixed);
374:     VecSet(step_fixed, 0.0);
375:     VecRestoreSubVector(S, active_fixed, &step_fixed);
376:   }
377:   return 0;
378: }

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

383:   Collective on Vec

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

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

395:   Level: developer

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


420:   VecCheckSameSize(X,1,XL,2);
421:   VecCheckSameSize(X,1,XU,3);
422:   VecCheckSameSize(X,1,Xout,4);

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

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

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