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: /*@
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
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: Level: developer
20: Notes:
21: `maskvalue` should usually be `0.0`, unless a pointwise divide will be used.
23: .seealso: `TaoMatGetSubMat()`, `TaoSubsetType`
24: @*/
25: PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced)
26: {
27: PetscInt nfull, nreduced, nreduced_local, rlow, rhigh, flow, fhigh;
28: PetscInt i, nlocal;
29: PetscReal *fv, *rv;
30: const PetscInt *s;
31: IS ident;
32: VecType vtype;
33: VecScatter scatter;
34: MPI_Comm comm;
36: PetscFunctionBegin;
40: PetscCall(VecGetSize(vfull, &nfull));
41: PetscCall(ISGetSize(is, &nreduced));
43: if (nreduced == nfull) {
44: PetscCall(VecDestroy(vreduced));
45: PetscCall(VecDuplicate(vfull, vreduced));
46: PetscCall(VecCopy(vfull, *vreduced));
47: } else {
48: switch (reduced_type) {
49: case TAO_SUBSET_SUBVEC:
50: PetscCall(VecGetType(vfull, &vtype));
51: PetscCall(VecGetOwnershipRange(vfull, &flow, &fhigh));
52: PetscCall(ISGetLocalSize(is, &nreduced_local));
53: PetscCall(PetscObjectGetComm((PetscObject)vfull, &comm));
54: if (*vreduced) PetscCall(VecDestroy(vreduced));
55: PetscCall(VecCreate(comm, vreduced));
56: PetscCall(VecSetType(*vreduced, vtype));
58: PetscCall(VecSetSizes(*vreduced, nreduced_local, nreduced));
59: PetscCall(VecGetOwnershipRange(*vreduced, &rlow, &rhigh));
60: PetscCall(ISCreateStride(comm, nreduced_local, rlow, 1, &ident));
61: PetscCall(VecScatterCreate(vfull, is, *vreduced, ident, &scatter));
62: PetscCall(VecScatterBegin(scatter, vfull, *vreduced, INSERT_VALUES, SCATTER_FORWARD));
63: PetscCall(VecScatterEnd(scatter, vfull, *vreduced, INSERT_VALUES, SCATTER_FORWARD));
64: PetscCall(VecScatterDestroy(&scatter));
65: PetscCall(ISDestroy(&ident));
66: break;
68: case TAO_SUBSET_MASK:
69: case TAO_SUBSET_MATRIXFREE:
70: /* vr[i] = vf[i] if i in is
71: vr[i] = 0 otherwise */
72: if (!*vreduced) PetscCall(VecDuplicate(vfull, vreduced));
74: PetscCall(VecSet(*vreduced, maskvalue));
75: PetscCall(ISGetLocalSize(is, &nlocal));
76: PetscCall(VecGetOwnershipRange(vfull, &flow, &fhigh));
77: PetscCall(VecGetArray(vfull, &fv));
78: PetscCall(VecGetArray(*vreduced, &rv));
79: PetscCall(ISGetIndices(is, &s));
80: PetscCheck(nlocal <= (fhigh - flow), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "IS local size %" PetscInt_FMT " > Vec local size %" PetscInt_FMT, nlocal, fhigh - flow);
81: for (i = 0; i < nlocal; ++i) rv[s[i] - flow] = fv[s[i] - flow];
82: PetscCall(ISRestoreIndices(is, &s));
83: PetscCall(VecRestoreArray(vfull, &fv));
84: PetscCall(VecRestoreArray(*vreduced, &rv));
85: break;
86: }
87: }
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: /*@
92: TaoMatGetSubMat - Gets a submatrix using the `IS`
94: Input Parameters:
95: + M - the full matrix (`n x n`)
96: . is - the index set for the submatrix (both row and column index sets need to be the same)
97: . v1 - work vector of dimension n, needed for `TAO_SUBSET_MASK` option
98: - subset_type - the method `Tao` is using for subsetting
100: Output Parameter:
101: . Msub - the submatrix
103: Level: developer
105: .seealso: `TaoVecGetSubVec()`, `TaoSubsetType`
106: @*/
107: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
108: {
109: IS iscomp;
110: PetscBool flg = PETSC_TRUE;
112: PetscFunctionBegin;
115: PetscCall(MatDestroy(Msub));
116: switch (subset_type) {
117: case TAO_SUBSET_SUBVEC:
118: PetscCall(MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub));
119: break;
121: case TAO_SUBSET_MASK:
122: /* Get Reduced Hessian
123: Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
124: Msub[i,j] = 0 if i!=j and i or j not in Free_Local
125: */
126: PetscObjectOptionsBegin((PetscObject)M);
127: PetscCall(PetscOptionsBool("-overwrite_hessian", "modify the existing hessian matrix when computing submatrices", "TaoSubsetType", flg, &flg, NULL));
128: PetscOptionsEnd();
129: if (flg) PetscCall(MatDuplicate(M, MAT_COPY_VALUES, Msub));
130: else {
131: /* Act on hessian directly (default) */
132: PetscCall(PetscObjectReference((PetscObject)M));
133: *Msub = M;
134: }
135: /* Save the diagonal to temporary vector */
136: PetscCall(MatGetDiagonal(*Msub, v1));
138: /* Zero out rows and columns */
139: PetscCall(ISComplementVec(is, v1, &iscomp));
141: /* Use v1 instead of 0 here because of PETSc bug */
142: PetscCall(MatZeroRowsColumnsIS(*Msub, iscomp, 1.0, v1, v1));
144: PetscCall(ISDestroy(&iscomp));
145: break;
146: case TAO_SUBSET_MATRIXFREE:
147: PetscCall(ISComplementVec(is, v1, &iscomp));
148: PetscCall(MatCreateSubMatrixFree(M, iscomp, iscomp, Msub));
149: PetscCall(ISDestroy(&iscomp));
150: break;
151: }
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*@C
156: TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper
157: bounds, as well as fixed variables where lower and upper bounds equal each other.
159: Input Parameters:
160: + X - solution vector
161: . XL - lower bound vector
162: . XU - upper bound vector
163: . G - unprojected gradient
164: . S - step direction with which the active bounds will be estimated
165: . W - work vector of type and size of `X`
166: - steplen - the step length at which the active bounds will be estimated (needs to be conservative)
168: Output Parameters:
169: + bound_tol - tolerance for the bound estimation
170: . active_lower - index set for active variables at the lower bound
171: . active_upper - index set for active variables at the upper bound
172: . active_fixed - index set for fixed variables
173: . active - index set for all active variables
174: - inactive - complementary index set for inactive variables
176: Level: developer
178: Notes:
179: This estimation is based on Bertsekas' method, with a built in diagonal scaling value of `1.0e-3`.
181: .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundSolution()`
182: @*/
183: 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)
184: {
185: PetscReal wnorm;
186: PetscReal zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0);
187: PetscInt i, n_isl = 0, n_isu = 0, n_isf = 0, n_isa = 0, n_isi = 0;
188: PetscInt N_isl, N_isu, N_isf, N_isa, N_isi;
189: PetscInt n, low, high, nDiff;
190: PetscInt *isl = NULL, *isu = NULL, *isf = NULL, *isa = NULL, *isi = NULL;
191: const PetscScalar *xl, *xu, *x, *g;
192: MPI_Comm comm = PetscObjectComm((PetscObject)X);
194: PetscFunctionBegin;
202: if (XL) PetscCheckSameType(X, 1, XL, 2);
203: if (XU) PetscCheckSameType(X, 1, XU, 3);
204: PetscCheckSameType(X, 1, G, 4);
205: PetscCheckSameType(X, 1, S, 5);
206: PetscCheckSameType(X, 1, W, 6);
207: if (XL) PetscCheckSameComm(X, 1, XL, 2);
208: if (XU) PetscCheckSameComm(X, 1, XU, 3);
209: PetscCheckSameComm(X, 1, G, 4);
210: PetscCheckSameComm(X, 1, S, 5);
211: PetscCheckSameComm(X, 1, W, 6);
212: if (XL) VecCheckSameSize(X, 1, XL, 2);
213: if (XU) VecCheckSameSize(X, 1, XU, 3);
214: VecCheckSameSize(X, 1, G, 4);
215: VecCheckSameSize(X, 1, S, 5);
216: VecCheckSameSize(X, 1, W, 6);
218: /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
219: PetscCall(VecCopy(X, W));
220: PetscCall(VecAXPBY(W, steplen, 1.0, S));
221: PetscCall(TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W));
222: PetscCall(VecAXPBY(W, 1.0, -1.0, X));
223: PetscCall(VecNorm(W, NORM_2, &wnorm));
224: *bound_tol = PetscMin(*bound_tol, wnorm);
226: /* Clear all index sets */
227: PetscCall(ISDestroy(active_lower));
228: PetscCall(ISDestroy(active_upper));
229: PetscCall(ISDestroy(active_fixed));
230: PetscCall(ISDestroy(active));
231: PetscCall(ISDestroy(inactive));
233: PetscCall(VecGetOwnershipRange(X, &low, &high));
234: PetscCall(VecGetLocalSize(X, &n));
235: if (!XL && !XU) {
236: PetscCall(ISCreateStride(comm, n, low, 1, inactive));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
239: if (n > 0) {
240: PetscCall(VecGetArrayRead(X, &x));
241: PetscCall(VecGetArrayRead(XL, &xl));
242: PetscCall(VecGetArrayRead(XU, &xu));
243: PetscCall(VecGetArrayRead(G, &g));
245: /* Loop over variables and categorize the indexes */
246: PetscCall(PetscMalloc1(n, &isl));
247: PetscCall(PetscMalloc1(n, &isu));
248: PetscCall(PetscMalloc1(n, &isf));
249: PetscCall(PetscMalloc1(n, &isa));
250: PetscCall(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;
255: ++n_isf;
256: isa[n_isa] = low + i;
257: ++n_isa;
258: } else if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + *bound_tol && g[i] > zero) {
259: /* Lower bounded variables */
260: isl[n_isl] = low + i;
261: ++n_isl;
262: isa[n_isa] = low + i;
263: ++n_isa;
264: } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - *bound_tol && g[i] < zero) {
265: /* Upper bounded variables */
266: isu[n_isu] = low + i;
267: ++n_isu;
268: isa[n_isa] = low + i;
269: ++n_isa;
270: } else {
271: /* Inactive variables */
272: isi[n_isi] = low + i;
273: ++n_isi;
274: }
275: }
277: PetscCall(VecRestoreArrayRead(X, &x));
278: PetscCall(VecRestoreArrayRead(XL, &xl));
279: PetscCall(VecRestoreArrayRead(XU, &xu));
280: PetscCall(VecRestoreArrayRead(G, &g));
281: }
283: /* Collect global sizes */
284: PetscCallMPI(MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm));
285: PetscCallMPI(MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm));
286: PetscCallMPI(MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm));
287: PetscCallMPI(MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm));
288: PetscCallMPI(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: PetscCall(ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower));
293: } else {
294: PetscCall(PetscFree(isl));
295: }
296: /* Create index set for upper bounded variables */
297: if (N_isu > 0) {
298: PetscCall(ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper));
299: } else {
300: PetscCall(PetscFree(isu));
301: }
302: /* Create index set for fixed variables */
303: if (N_isf > 0) {
304: PetscCall(ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed));
305: } else {
306: PetscCall(PetscFree(isf));
307: }
308: /* Create index set for all actively bounded variables */
309: if (N_isa > 0) {
310: PetscCall(ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active));
311: } else {
312: PetscCall(PetscFree(isa));
313: }
314: /* Create index set for all inactive variables */
315: if (N_isi > 0) {
316: PetscCall(ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive));
317: } else {
318: PetscCall(PetscFree(isi));
319: }
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@
324: TaoBoundStep - Ensures the correct zero or adjusted step direction values for active
325: variables.
327: Input Parameters:
328: + X - solution vector
329: . XL - lower bound vector
330: . XU - upper bound vector
331: . active_lower - index set for lower bounded active variables
332: . active_upper - index set for lower bounded active variables
333: . active_fixed - index set for fixed active variables
334: - scale - amplification factor for the step that needs to be taken on actively bounded variables
336: Output Parameter:
337: . S - step direction to be modified
339: Level: developer
341: .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundSolution()`
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: Vec step_lower, step_upper, step_fixed;
346: Vec x_lower, x_upper;
347: Vec bound_lower, bound_upper;
349: PetscFunctionBegin;
350: /* Adjust step for variables at the estimated lower bound */
351: if (active_lower) {
352: PetscCall(VecGetSubVector(S, active_lower, &step_lower));
353: PetscCall(VecGetSubVector(X, active_lower, &x_lower));
354: PetscCall(VecGetSubVector(XL, active_lower, &bound_lower));
355: PetscCall(VecCopy(bound_lower, step_lower));
356: PetscCall(VecAXPY(step_lower, -1.0, x_lower));
357: PetscCall(VecScale(step_lower, scale));
358: PetscCall(VecRestoreSubVector(S, active_lower, &step_lower));
359: PetscCall(VecRestoreSubVector(X, active_lower, &x_lower));
360: PetscCall(VecRestoreSubVector(XL, active_lower, &bound_lower));
361: }
363: /* Adjust step for the variables at the estimated upper bound */
364: if (active_upper) {
365: PetscCall(VecGetSubVector(S, active_upper, &step_upper));
366: PetscCall(VecGetSubVector(X, active_upper, &x_upper));
367: PetscCall(VecGetSubVector(XU, active_upper, &bound_upper));
368: PetscCall(VecCopy(bound_upper, step_upper));
369: PetscCall(VecAXPY(step_upper, -1.0, x_upper));
370: PetscCall(VecScale(step_upper, scale));
371: PetscCall(VecRestoreSubVector(S, active_upper, &step_upper));
372: PetscCall(VecRestoreSubVector(X, active_upper, &x_upper));
373: PetscCall(VecRestoreSubVector(XU, active_upper, &bound_upper));
374: }
376: /* Zero out step for fixed variables */
377: if (active_fixed) {
378: PetscCall(VecGetSubVector(S, active_fixed, &step_fixed));
379: PetscCall(VecSet(step_fixed, 0.0));
380: PetscCall(VecRestoreSubVector(S, active_fixed, &step_fixed));
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance.
388: Collective
390: Input Parameters:
391: + X - solution vector
392: . XL - lower bound vector
393: . XU - upper bound vector
394: - bound_tol - absolute tolerance in enforcing the bound
396: Output Parameters:
397: + nDiff - total number of vector entries that have been bounded
398: - Xout - modified solution vector satisfying bounds to `bound_tol`
400: Level: developer
402: .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundStep()`
403: @*/
404: PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout)
405: {
406: PetscInt i, n, low, high, nDiff_loc = 0;
407: PetscScalar *xout;
408: const PetscScalar *x, *xl, *xu;
410: PetscFunctionBegin;
415: if (!XL && !XU) {
416: PetscCall(VecCopy(X, Xout));
417: *nDiff = 0.0;
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
420: PetscCheckSameType(X, 1, XL, 2);
421: PetscCheckSameType(X, 1, XU, 3);
422: PetscCheckSameType(X, 1, Xout, 6);
423: PetscCheckSameComm(X, 1, XL, 2);
424: PetscCheckSameComm(X, 1, XU, 3);
425: PetscCheckSameComm(X, 1, Xout, 6);
426: VecCheckSameSize(X, 1, XL, 2);
427: VecCheckSameSize(X, 1, XU, 3);
428: VecCheckSameSize(X, 1, Xout, 4);
430: PetscCall(VecGetOwnershipRange(X, &low, &high));
431: PetscCall(VecGetLocalSize(X, &n));
432: if (n > 0) {
433: PetscCall(VecGetArrayRead(X, &x));
434: PetscCall(VecGetArrayRead(XL, &xl));
435: PetscCall(VecGetArrayRead(XU, &xu));
436: PetscCall(VecGetArray(Xout, &xout));
438: for (i = 0; i < n; ++i) {
439: if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + bound_tol) {
440: xout[i] = xl[i];
441: ++nDiff_loc;
442: } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - bound_tol) {
443: xout[i] = xu[i];
444: ++nDiff_loc;
445: }
446: }
448: PetscCall(VecRestoreArrayRead(X, &x));
449: PetscCall(VecRestoreArrayRead(XL, &xl));
450: PetscCall(VecRestoreArrayRead(XU, &xu));
451: PetscCall(VecRestoreArray(Xout, &xout));
452: }
453: PetscCallMPI(MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X)));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }