Actual source code: isutil.c
petsc-3.14.6 2021-03-30
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 Parameters:
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: {
26: PetscInt nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh;
27: PetscInt i,nlocal;
28: PetscReal *fv,*rv;
29: const PetscInt *s;
30: IS ident;
31: VecType vtype;
32: VecScatter scatter;
33: MPI_Comm comm;
39: VecGetSize(vfull, &nfull);
40: ISGetSize(is, &nreduced);
42: if (nreduced == nfull) {
43: VecDestroy(vreduced);
44: VecDuplicate(vfull,vreduced);
45: VecCopy(vfull,*vreduced);
46: } else {
47: switch (reduced_type) {
48: case TAO_SUBSET_SUBVEC:
49: VecGetType(vfull,&vtype);
50: VecGetOwnershipRange(vfull,&flow,&fhigh);
51: ISGetLocalSize(is,&nreduced_local);
52: PetscObjectGetComm((PetscObject)vfull,&comm);
53: if (*vreduced) {
54: VecDestroy(vreduced);
55: }
56: VecCreate(comm,vreduced);
57: VecSetType(*vreduced,vtype);
59: VecSetSizes(*vreduced,nreduced_local,nreduced);
60: VecGetOwnershipRange(*vreduced,&rlow,&rhigh);
61: ISCreateStride(comm,nreduced_local,rlow,1,&ident);
62: VecScatterCreate(vfull,is,*vreduced,ident,&scatter);
63: VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
64: VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
65: VecScatterDestroy(&scatter);
66: ISDestroy(&ident);
67: break;
69: case TAO_SUBSET_MASK:
70: case TAO_SUBSET_MATRIXFREE:
71: /* vr[i] = vf[i] if i in is
72: vr[i] = 0 otherwise */
73: if (!*vreduced) {
74: VecDuplicate(vfull,vreduced);
75: }
77: VecSet(*vreduced,maskvalue);
78: ISGetLocalSize(is,&nlocal);
79: VecGetOwnershipRange(vfull,&flow,&fhigh);
80: VecGetArray(vfull,&fv);
81: VecGetArray(*vreduced,&rv);
82: ISGetIndices(is,&s);
83: if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_SELF,1,"IS local size %D > Vec local size %D",nlocal,fhigh-flow);
84: for (i=0;i<nlocal;++i) {
85: rv[s[i]-flow] = fv[s[i]-flow];
86: }
87: ISRestoreIndices(is,&s);
88: VecRestoreArray(vfull,&fv);
89: VecRestoreArray(*vreduced,&rv);
90: break;
91: }
92: }
93: return(0);
94: }
96: /*@C
97: TaoMatGetSubMat - Gets a submatrix using the IS
99: Input Parameters:
100: + M - the full matrix (n x n)
101: . is - the index set for the submatrix (both row and column index sets need to be the same)
102: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
103: - subset_type <TAO_SUBSET_SUBVEC,TAO_SUBSET_MASK,TAO_SUBSET_MATRIXFREE> - the method TAO is using for subsetting
105: Output Parameters:
106: . Msub - the submatrix
108: Level: developer
109: @*/
110: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
111: {
113: IS iscomp;
114: PetscBool flg = PETSC_TRUE;
119: MatDestroy(Msub);
120: switch (subset_type) {
121: case TAO_SUBSET_SUBVEC:
122: MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
123: break;
125: case TAO_SUBSET_MASK:
126: /* Get Reduced Hessian
127: Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
128: Msub[i,j] = 0 if i!=j and i or j not in Free_Local
129: */
130: PetscObjectOptionsBegin((PetscObject)M);
131: PetscOptionsBool("-overwrite_hessian","modify the existing hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);
132: PetscOptionsEnd();
133: if (flg) {
134: MatDuplicate(M, MAT_COPY_VALUES, Msub);
135: } else {
136: /* Act on hessian directly (default) */
137: PetscObjectReference((PetscObject)M);
138: *Msub = M;
139: }
140: /* Save the diagonal to temporary vector */
141: MatGetDiagonal(*Msub,v1);
143: /* Zero out rows and columns */
144: ISComplementVec(is,v1,&iscomp);
146: /* Use v1 instead of 0 here because of PETSc bug */
147: MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);
149: ISDestroy(&iscomp);
150: break;
151: case TAO_SUBSET_MATRIXFREE:
152: ISComplementVec(is,v1,&iscomp);
153: MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
154: ISDestroy(&iscomp);
155: break;
156: }
157: return(0);
158: }
160: /*@C
161: TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper
162: bounds, as well as fixed variables where lower and upper bounds equal each other.
164: Input Parameters:
165: + X - solution vector
166: . XL - lower bound vector
167: . XU - upper bound vector
168: . G - unprojected gradient
169: . S - step direction with which the active bounds will be estimated
170: . W - work vector of type and size of X
171: - steplen - the step length at which the active bounds will be estimated (needs to be conservative)
173: Output Parameters:
174: + bound_tol - tolerance for for the bound estimation
175: . active_lower - index set for active variables at the lower bound
176: . active_upper - index set for active variables at the upper bound
177: . active_fixed - index set for fixed variables
178: . active - index set for all active variables
179: - inactive - complementary index set for inactive variables
181: Notes:
182: This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3.
184: Level: developer
185: @*/
186: PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol,
187: IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive)
188: {
189: PetscErrorCode ierr;
190: PetscReal wnorm;
191: PetscReal zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
192: PetscInt i, n_isl=0, n_isu=0, n_isf=0, n_isa=0, n_isi=0;
193: PetscInt N_isl, N_isu, N_isf, N_isa, N_isi;
194: PetscInt n, low, high, nDiff;
195: PetscInt *isl=NULL, *isu=NULL, *isf=NULL, *isa=NULL, *isi=NULL;
196: const PetscScalar *xl, *xu, *x, *g;
197: MPI_Comm comm = PetscObjectComm((PetscObject)X);
223: VecCheckSameSize(X,1,XL,2);
224: VecCheckSameSize(X,1,XU,3);
225: VecCheckSameSize(X,1,G,4);
226: VecCheckSameSize(X,1,S,5);
227: VecCheckSameSize(X,1,W,6);
229: /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
230: VecCopy(X, W);
231: VecAXPBY(W, steplen, 1.0, S);
232: TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W);
233: VecAXPBY(W, 1.0, -1.0, X);
234: VecNorm(W, NORM_2, &wnorm);
235: *bound_tol = PetscMin(*bound_tol, wnorm);
237: VecGetOwnershipRange(X, &low, &high);
238: VecGetLocalSize(X, &n);
239: if (n>0){
240: VecGetArrayRead(X, &x);
241: VecGetArrayRead(XL, &xl);
242: VecGetArrayRead(XU, &xu);
243: VecGetArrayRead(G, &g);
245: /* Loop over variables and categorize the indexes */
246: PetscMalloc1(n, &isl);
247: PetscMalloc1(n, &isu);
248: PetscMalloc1(n, &isf);
249: PetscMalloc1(n, &isa);
250: 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; ++n_isf;
255: isa[n_isa]=low+i; ++n_isa;
256: } else if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + *bound_tol) && (g[i] > zero)) {
257: /* Lower bounded variables */
258: isl[n_isl]=low+i; ++n_isl;
259: isa[n_isa]=low+i; ++n_isa;
260: } else if ((xu[i] < PETSC_INFINITY) && (x[i] >= xu[i] - *bound_tol) && (g[i] < zero)) {
261: /* Upper bounded variables */
262: isu[n_isu]=low+i; ++n_isu;
263: isa[n_isa]=low+i; ++n_isa;
264: } else {
265: /* Inactive variables */
266: isi[n_isi]=low+i; ++n_isi;
267: }
268: }
270: VecRestoreArrayRead(X, &x);
271: VecRestoreArrayRead(XL, &xl);
272: VecRestoreArrayRead(XU, &xu);
273: VecRestoreArrayRead(G, &g);
274: }
276: /* Clear all index sets */
277: ISDestroy(active_lower);
278: ISDestroy(active_upper);
279: ISDestroy(active_fixed);
280: ISDestroy(active);
281: ISDestroy(inactive);
283: /* Collect global sizes */
284: MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);
285: MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);
286: MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);
287: MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);
288: 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: ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower);
293: } else {
294: PetscFree(isl);
295: }
296: /* Create index set for upper bounded variables */
297: if (N_isu > 0) {
298: ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper);
299: } else {
300: PetscFree(isu);
301: }
302: /* Create index set for fixed variables */
303: if (N_isf > 0) {
304: ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed);
305: } else {
306: PetscFree(isf);
307: }
308: /* Create index set for all actively bounded variables */
309: if (N_isa > 0) {
310: ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active);
311: } else {
312: PetscFree(isa);
313: }
314: /* Create index set for all inactive variables */
315: if (N_isi > 0) {
316: ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive);
317: } else {
318: PetscFree(isi);
319: }
321: /* Clean up and exit */
322: return(0);
323: }
325: /*@C
326: TaoBoundStep - Ensures the correct zero or adjusted step direction
327: values for active variables.
329: Input Parameters:
330: + X - solution vector
331: . XL - lower bound vector
332: . XU - upper bound vector
333: . active_lower - index set for lower bounded active variables
334: . active_upper - index set for lower bounded active variables
335: . active_fixed - index set for fixed active variables
336: - scale - amplification factor for the step that needs to be taken on actively bounded variables
338: Output Parameters:
339: . S - step direction to be modified
341: Level: developer
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: PetscErrorCode ierr;
347: Vec step_lower, step_upper, step_fixed;
348: Vec x_lower, x_upper;
349: Vec bound_lower, bound_upper;
352: /* Adjust step for variables at the estimated lower bound */
353: if (active_lower) {
354: VecGetSubVector(S, active_lower, &step_lower);
355: VecGetSubVector(X, active_lower, &x_lower);
356: VecGetSubVector(XL, active_lower, &bound_lower);
357: VecCopy(bound_lower, step_lower);
358: VecAXPY(step_lower, -1.0, x_lower);
359: VecScale(step_lower, scale);
360: VecRestoreSubVector(S, active_lower, &step_lower);
361: VecRestoreSubVector(X, active_lower, &x_lower);
362: VecRestoreSubVector(XL, active_lower, &bound_lower);
363: }
365: /* Adjust step for the variables at the estimated upper bound */
366: if (active_upper) {
367: VecGetSubVector(S, active_upper, &step_upper);
368: VecGetSubVector(X, active_upper, &x_upper);
369: VecGetSubVector(XU, active_upper, &bound_upper);
370: VecCopy(bound_upper, step_upper);
371: VecAXPY(step_upper, -1.0, x_upper);
372: VecScale(step_upper, scale);
373: VecRestoreSubVector(S, active_upper, &step_upper);
374: VecRestoreSubVector(X, active_upper, &x_upper);
375: VecRestoreSubVector(XU, active_upper, &bound_upper);
376: }
378: /* Zero out step for fixed variables */
379: if (active_fixed) {
380: VecGetSubVector(S, active_fixed, &step_fixed);
381: VecSet(step_fixed, 0.0);
382: VecRestoreSubVector(S, active_fixed, &step_fixed);
383: }
384: return(0);
385: }
387: /*@C
388: TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance.
390: Collective on Vec
392: Input Parameters:
393: + X - solution vector
394: . XL - lower bound vector
395: . XU - upper bound vector
396: - bound_tol - absolute tolerance in enforcing the bound
398: Output Parameters:
399: + nDiff - total number of vector entries that have been bounded
400: - Xout - modified solution vector satisfying bounds to bound_tol
402: Level: developer
404: .seealso: TAOBNCG, TAOBNTL, TAOBNTR
405: @*/
406: PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout)
407: {
408: PetscErrorCode ierr;
409: PetscInt i,n,low,high,nDiff_loc=0;
410: PetscScalar *xout;
411: const PetscScalar *x,*xl,*xu;
429: VecCheckSameSize(X,1,XL,2);
430: VecCheckSameSize(X,1,XU,3);
431: VecCheckSameSize(X,1,Xout,4);
433: VecGetOwnershipRange(X,&low,&high);
434: VecGetLocalSize(X,&n);
435: if (n>0){
436: VecGetArrayRead(X, &x);
437: VecGetArrayRead(XL, &xl);
438: VecGetArrayRead(XU, &xu);
439: VecGetArray(Xout, &xout);
441: for (i=0;i<n;++i){
442: if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + bound_tol)) {
443: xout[i] = xl[i]; ++nDiff_loc;
444: } else if ((xu[i] < PETSC_INFINITY) && (x[i] >= xu[i] - bound_tol)) {
445: xout[i] = xu[i]; ++nDiff_loc;
446: }
447: }
449: VecRestoreArrayRead(X, &x);
450: VecRestoreArrayRead(XL, &xl);
451: VecRestoreArrayRead(XU, &xu);
452: VecRestoreArray(Xout, &xout);
453: }
454: MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X));
455: return(0);
456: }