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