1 #include <petsctao.h> /*I "petsctao.h" I*/
2 #include <petsc/private/vecimpl.h>
3 #include <petsc/private/taoimpl.h>
4 #include <../src/tao/matrix/submatfree.h>
5
6 /*@
7 TaoVecGetSubVec - Gets a subvector using the `IS`
8
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`)
14
15 Output Parameter:
16 . vreduced - the subvector
17
18 Level: developer
19
20 Notes:
21 `maskvalue` should usually be `0.0`, unless a pointwise divide will be used.
22
23 .seealso: `TaoMatGetSubMat()`, `TaoSubsetType`
24 @*/
TaoVecGetSubVec(Vec vfull,IS is,TaoSubsetType reduced_type,PetscReal maskvalue,Vec * vreduced)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;
35
36 PetscFunctionBegin;
37 PetscValidHeaderSpecific(vfull, VEC_CLASSID, 1);
38 PetscValidHeaderSpecific(is, IS_CLASSID, 2);
39
40 PetscCall(VecGetSize(vfull, &nfull));
41 PetscCall(ISGetSize(is, &nreduced));
42
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));
57
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;
67
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));
73
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 }
90
91 /*@
92 TaoMatGetSubMat - Gets a submatrix using the `IS`
93
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
99
100 Output Parameter:
101 . Msub - the submatrix
102
103 Level: developer
104
105 .seealso: `TaoVecGetSubVec()`, `TaoSubsetType`
106 @*/
TaoMatGetSubMat(Mat M,IS is,Vec v1,TaoSubsetType subset_type,Mat * Msub)107 PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
108 {
109 IS iscomp;
110 PetscBool flg = PETSC_TRUE;
111
112 PetscFunctionBegin;
113 PetscValidHeaderSpecific(M, MAT_CLASSID, 1);
114 PetscValidHeaderSpecific(is, IS_CLASSID, 2);
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;
120
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) {
130 PetscCall(MatDuplicate(M, MAT_COPY_VALUES, Msub));
131 } else {
132 /* Act on hessian directly (default) */
133 PetscCall(PetscObjectReference((PetscObject)M));
134 *Msub = M;
135 }
136 /* Save the diagonal to temporary vector */
137 PetscCall(MatGetDiagonal(*Msub, v1));
138
139 /* Zero out rows and columns */
140 PetscCall(ISComplementVec(is, v1, &iscomp));
141
142 /* Use v1 instead of 0 here because of PETSc bug */
143 PetscCall(MatZeroRowsColumnsIS(*Msub, iscomp, 1.0, v1, v1));
144
145 PetscCall(ISDestroy(&iscomp));
146 break;
147 case TAO_SUBSET_MATRIXFREE:
148 PetscCall(ISComplementVec(is, v1, &iscomp));
149 PetscCall(MatCreateSubMatrixFree(M, iscomp, iscomp, Msub));
150 PetscCall(ISDestroy(&iscomp));
151 break;
152 }
153 PetscFunctionReturn(PETSC_SUCCESS);
154 }
155
156 /*@C
157 TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper
158 bounds, as well as fixed variables where lower and upper bounds equal each other.
159
160 Input Parameters:
161 + X - solution vector
162 . XL - lower bound vector
163 . XU - upper bound vector
164 . G - unprojected gradient
165 . S - step direction with which the active bounds will be estimated
166 . W - work vector of type and size of `X`
167 - steplen - the step length at which the active bounds will be estimated (needs to be conservative)
168
169 Output Parameters:
170 + bound_tol - tolerance for the bound estimation
171 . active_lower - index set for active variables at the lower bound
172 . active_upper - index set for active variables at the upper bound
173 . active_fixed - index set for fixed variables
174 . active - index set for all active variables
175 - inactive - complementary index set for inactive variables
176
177 Level: developer
178
179 Notes:
180 This estimation is based on Bertsekas' method, with a built in diagonal scaling value of `1.0e-3`.
181
182 .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundSolution()`
183 @*/
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 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)
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);
194
195 PetscFunctionBegin;
196 PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
197 if (XL) PetscValidHeaderSpecific(XL, VEC_CLASSID, 2);
198 if (XU) PetscValidHeaderSpecific(XU, VEC_CLASSID, 3);
199 PetscValidHeaderSpecific(G, VEC_CLASSID, 4);
200 PetscValidHeaderSpecific(S, VEC_CLASSID, 5);
201 PetscValidHeaderSpecific(W, VEC_CLASSID, 6);
202
203 if (XL) PetscCheckSameType(X, 1, XL, 2);
204 if (XU) PetscCheckSameType(X, 1, XU, 3);
205 PetscCheckSameType(X, 1, G, 4);
206 PetscCheckSameType(X, 1, S, 5);
207 PetscCheckSameType(X, 1, W, 6);
208 if (XL) PetscCheckSameComm(X, 1, XL, 2);
209 if (XU) PetscCheckSameComm(X, 1, XU, 3);
210 PetscCheckSameComm(X, 1, G, 4);
211 PetscCheckSameComm(X, 1, S, 5);
212 PetscCheckSameComm(X, 1, W, 6);
213 if (XL) VecCheckSameSize(X, 1, XL, 2);
214 if (XU) VecCheckSameSize(X, 1, XU, 3);
215 VecCheckSameSize(X, 1, G, 4);
216 VecCheckSameSize(X, 1, S, 5);
217 VecCheckSameSize(X, 1, W, 6);
218
219 /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
220 PetscCall(VecCopy(X, W));
221 PetscCall(VecAXPBY(W, steplen, 1.0, S));
222 PetscCall(TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W));
223 PetscCall(VecAXPBY(W, 1.0, -1.0, X));
224 PetscCall(VecNorm(W, NORM_2, &wnorm));
225 *bound_tol = PetscMin(*bound_tol, wnorm);
226
227 /* Clear all index sets */
228 PetscCall(ISDestroy(active_lower));
229 PetscCall(ISDestroy(active_upper));
230 PetscCall(ISDestroy(active_fixed));
231 PetscCall(ISDestroy(active));
232 PetscCall(ISDestroy(inactive));
233
234 PetscCall(VecGetOwnershipRange(X, &low, &high));
235 PetscCall(VecGetLocalSize(X, &n));
236 if (!XL && !XU) {
237 PetscCall(ISCreateStride(comm, n, low, 1, inactive));
238 PetscFunctionReturn(PETSC_SUCCESS);
239 }
240 if (n > 0) {
241 PetscCall(VecGetArrayRead(X, &x));
242 PetscCall(VecGetArrayRead(XL, &xl));
243 PetscCall(VecGetArrayRead(XU, &xu));
244 PetscCall(VecGetArrayRead(G, &g));
245
246 /* Loop over variables and categorize the indexes */
247 PetscCall(PetscMalloc1(n, &isl));
248 PetscCall(PetscMalloc1(n, &isu));
249 PetscCall(PetscMalloc1(n, &isf));
250 PetscCall(PetscMalloc1(n, &isa));
251 PetscCall(PetscMalloc1(n, &isi));
252 for (i = 0; i < n; ++i) {
253 if (xl[i] == xu[i]) {
254 /* Fixed variables */
255 isf[n_isf] = low + i;
256 ++n_isf;
257 isa[n_isa] = low + i;
258 ++n_isa;
259 } else if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + *bound_tol && g[i] > zero) {
260 /* Lower bounded variables */
261 isl[n_isl] = low + i;
262 ++n_isl;
263 isa[n_isa] = low + i;
264 ++n_isa;
265 } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - *bound_tol && g[i] < zero) {
266 /* Upper bounded variables */
267 isu[n_isu] = low + i;
268 ++n_isu;
269 isa[n_isa] = low + i;
270 ++n_isa;
271 } else {
272 /* Inactive variables */
273 isi[n_isi] = low + i;
274 ++n_isi;
275 }
276 }
277
278 PetscCall(VecRestoreArrayRead(X, &x));
279 PetscCall(VecRestoreArrayRead(XL, &xl));
280 PetscCall(VecRestoreArrayRead(XU, &xu));
281 PetscCall(VecRestoreArrayRead(G, &g));
282 }
283
284 /* Collect global sizes */
285 PetscCallMPI(MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm));
286 PetscCallMPI(MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm));
287 PetscCallMPI(MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm));
288 PetscCallMPI(MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm));
289 PetscCallMPI(MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm));
290
291 /* Create index set for lower bounded variables */
292 if (N_isl > 0) {
293 PetscCall(ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower));
294 } else {
295 PetscCall(PetscFree(isl));
296 }
297 /* Create index set for upper bounded variables */
298 if (N_isu > 0) {
299 PetscCall(ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper));
300 } else {
301 PetscCall(PetscFree(isu));
302 }
303 /* Create index set for fixed variables */
304 if (N_isf > 0) {
305 PetscCall(ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed));
306 } else {
307 PetscCall(PetscFree(isf));
308 }
309 /* Create index set for all actively bounded variables */
310 if (N_isa > 0) {
311 PetscCall(ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active));
312 } else {
313 PetscCall(PetscFree(isa));
314 }
315 /* Create index set for all inactive variables */
316 if (N_isi > 0) {
317 PetscCall(ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive));
318 } else {
319 PetscCall(PetscFree(isi));
320 }
321 PetscFunctionReturn(PETSC_SUCCESS);
322 }
323
324 /*@
325 TaoBoundStep - Ensures the correct zero or adjusted step direction values for active
326 variables.
327
328 Input Parameters:
329 + X - solution vector
330 . XL - lower bound vector
331 . XU - upper bound vector
332 . active_lower - index set for lower bounded active variables
333 . active_upper - index set for lower bounded active variables
334 . active_fixed - index set for fixed active variables
335 - scale - amplification factor for the step that needs to be taken on actively bounded variables
336
337 Output Parameter:
338 . S - step direction to be modified
339
340 Level: developer
341
342 .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundSolution()`
343 @*/
TaoBoundStep(Vec X,Vec XL,Vec XU,IS active_lower,IS active_upper,IS active_fixed,PetscReal scale,Vec S)344 PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S)
345 {
346 Vec step_lower, step_upper, step_fixed;
347 Vec x_lower, x_upper;
348 Vec bound_lower, bound_upper;
349
350 PetscFunctionBegin;
351 /* Adjust step for variables at the estimated lower bound */
352 if (active_lower) {
353 PetscCall(VecGetSubVector(S, active_lower, &step_lower));
354 PetscCall(VecGetSubVector(X, active_lower, &x_lower));
355 PetscCall(VecGetSubVector(XL, active_lower, &bound_lower));
356 PetscCall(VecCopy(bound_lower, step_lower));
357 PetscCall(VecAXPY(step_lower, -1.0, x_lower));
358 PetscCall(VecScale(step_lower, scale));
359 PetscCall(VecRestoreSubVector(S, active_lower, &step_lower));
360 PetscCall(VecRestoreSubVector(X, active_lower, &x_lower));
361 PetscCall(VecRestoreSubVector(XL, active_lower, &bound_lower));
362 }
363
364 /* Adjust step for the variables at the estimated upper bound */
365 if (active_upper) {
366 PetscCall(VecGetSubVector(S, active_upper, &step_upper));
367 PetscCall(VecGetSubVector(X, active_upper, &x_upper));
368 PetscCall(VecGetSubVector(XU, active_upper, &bound_upper));
369 PetscCall(VecCopy(bound_upper, step_upper));
370 PetscCall(VecAXPY(step_upper, -1.0, x_upper));
371 PetscCall(VecScale(step_upper, scale));
372 PetscCall(VecRestoreSubVector(S, active_upper, &step_upper));
373 PetscCall(VecRestoreSubVector(X, active_upper, &x_upper));
374 PetscCall(VecRestoreSubVector(XU, active_upper, &bound_upper));
375 }
376
377 /* Zero out step for fixed variables */
378 if (active_fixed) {
379 PetscCall(VecGetSubVector(S, active_fixed, &step_fixed));
380 PetscCall(VecSet(step_fixed, 0.0));
381 PetscCall(VecRestoreSubVector(S, active_fixed, &step_fixed));
382 }
383 PetscFunctionReturn(PETSC_SUCCESS);
384 }
385
386 /*@
387 TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance.
388
389 Collective
390
391 Input Parameters:
392 + X - solution vector
393 . XL - lower bound vector
394 . XU - upper bound vector
395 - bound_tol - absolute tolerance in enforcing the bound
396
397 Output Parameters:
398 + nDiff - total number of vector entries that have been bounded
399 - Xout - modified solution vector satisfying bounds to `bound_tol`
400
401 Level: developer
402
403 .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundStep()`
404 @*/
TaoBoundSolution(Vec X,Vec XL,Vec XU,PetscReal bound_tol,PetscInt * nDiff,Vec Xout)405 PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout)
406 {
407 PetscInt i, n, low, high, nDiff_loc = 0;
408 PetscScalar *xout;
409 const PetscScalar *x, *xl, *xu;
410
411 PetscFunctionBegin;
412 PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
413 if (XL) PetscValidHeaderSpecific(XL, VEC_CLASSID, 2);
414 if (XU) PetscValidHeaderSpecific(XU, VEC_CLASSID, 3);
415 PetscValidHeaderSpecific(Xout, VEC_CLASSID, 6);
416 if (!XL && !XU) {
417 PetscCall(VecCopy(X, Xout));
418 *nDiff = 0.0;
419 PetscFunctionReturn(PETSC_SUCCESS);
420 }
421 PetscCheckSameType(X, 1, XL, 2);
422 PetscCheckSameType(X, 1, XU, 3);
423 PetscCheckSameType(X, 1, Xout, 6);
424 PetscCheckSameComm(X, 1, XL, 2);
425 PetscCheckSameComm(X, 1, XU, 3);
426 PetscCheckSameComm(X, 1, Xout, 6);
427 VecCheckSameSize(X, 1, XL, 2);
428 VecCheckSameSize(X, 1, XU, 3);
429 VecCheckSameSize(X, 1, Xout, 4);
430
431 PetscCall(VecGetOwnershipRange(X, &low, &high));
432 PetscCall(VecGetLocalSize(X, &n));
433 if (n > 0) {
434 PetscCall(VecGetArrayRead(X, &x));
435 PetscCall(VecGetArrayRead(XL, &xl));
436 PetscCall(VecGetArrayRead(XU, &xu));
437 PetscCall(VecGetArray(Xout, &xout));
438
439 for (i = 0; i < n; ++i) {
440 if (xl[i] > PETSC_NINFINITY && x[i] <= xl[i] + bound_tol) {
441 xout[i] = xl[i];
442 ++nDiff_loc;
443 } else if (xu[i] < PETSC_INFINITY && x[i] >= xu[i] - bound_tol) {
444 xout[i] = xu[i];
445 ++nDiff_loc;
446 }
447 }
448
449 PetscCall(VecRestoreArrayRead(X, &x));
450 PetscCall(VecRestoreArrayRead(XL, &xl));
451 PetscCall(VecRestoreArrayRead(XU, &xu));
452 PetscCall(VecRestoreArray(Xout, &xout));
453 }
454 PetscCallMPI(MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X)));
455 PetscFunctionReturn(PETSC_SUCCESS);
456 }
457