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