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