Actual source code: isutil.c
petsc-3.12.5 2020-03-29
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: @*/
22: PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced)
23: {
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;
38: VecGetSize(vfull, &nfull);
39: ISGetSize(is, &nreduced);
41: if (nreduced == nfull) {
42: VecDestroy(vreduced);
43: VecDuplicate(vfull,vreduced);
44: VecCopy(vfull,*vreduced);
45: } else {
46: switch (reduced_type) {
47: case TAO_SUBSET_SUBVEC:
48: VecGetType(vfull,&vtype);
49: VecGetOwnershipRange(vfull,&flow,&fhigh);
50: ISGetLocalSize(is,&nreduced_local);
51: PetscObjectGetComm((PetscObject)vfull,&comm);
52: if (*vreduced) {
53: VecDestroy(vreduced);
54: }
55: VecCreate(comm,vreduced);
56: VecSetType(*vreduced,vtype);
58: VecSetSizes(*vreduced,nreduced_local,nreduced);
59: VecGetOwnershipRange(*vreduced,&rlow,&rhigh);
60: ISCreateStride(comm,nreduced_local,rlow,1,&ident);
61: VecScatterCreate(vfull,is,*vreduced,ident,&scatter);
62: VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
63: VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
64: VecScatterDestroy(&scatter);
65: 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) {
73: VecDuplicate(vfull,vreduced);
74: }
76: VecSet(*vreduced,maskvalue);
77: ISGetLocalSize(is,&nlocal);
78: VecGetOwnershipRange(vfull,&flow,&fhigh);
79: VecGetArray(vfull,&fv);
80: VecGetArray(*vreduced,&rv);
81: ISGetIndices(is,&s);
82: if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow);
83: for (i=0;i<nlocal;++i) {
84: rv[s[i]-flow] = fv[s[i]-flow];
85: }
86: ISRestoreIndices(is,&s);
87: VecRestoreArray(vfull,&fv);
88: VecRestoreArray(*vreduced,&rv);
89: break;
90: }
91: }
92: return(0);
93: }
95: /*@C
96: TaoMatGetSubMat - Gets a submatrix using the IS
98: Input Parameters:
99: + M - the full matrix (n x n)
100: . is - the index set for the submatrix (both row and column index sets need to be the same)
101: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
102: - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,
103: TAO_SUBSET_MATRIXFREE)
105: Output Parameters:
106: . Msub - the submatrix
107: @*/
108: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
109: {
111: IS iscomp;
112: PetscBool flg = PETSC_TRUE;
117: MatDestroy(Msub);
118: switch (subset_type) {
119: case TAO_SUBSET_SUBVEC:
120: MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
121: break;
123: case TAO_SUBSET_MASK:
124: /* Get Reduced Hessian
125: Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
126: Msub[i,j] = 0 if i!=j and i or j not in Free_Local
127: */
128: PetscOptionsBegin(PetscObjectComm((PetscObject)M),NULL,NULL,NULL);
129: PetscOptionsBool("-overwrite_hessian","modify the existing hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);
130: PetscOptionsEnd();
131: if (flg) {
132: MatDuplicate(M, MAT_COPY_VALUES, Msub);
133: } else {
134: /* Act on hessian directly (default) */
135: PetscObjectReference((PetscObject)M);
136: *Msub = M;
137: }
138: /* Save the diagonal to temporary vector */
139: MatGetDiagonal(*Msub,v1);
141: /* Zero out rows and columns */
142: ISComplementVec(is,v1,&iscomp);
144: /* Use v1 instead of 0 here because of PETSc bug */
145: MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);
147: ISDestroy(&iscomp);
148: break;
149: case TAO_SUBSET_MATRIXFREE:
150: ISComplementVec(is,v1,&iscomp);
151: MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
152: ISDestroy(&iscomp);
153: break;
154: }
155: return(0);
156: }
158: /*@C
159: TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper
160: bounds, as well as fixed variables where lower and upper bounds equal each other.
162: Input Parameters:
163: + X - solution vector
164: . XL - lower bound vector
165: . XU - upper bound vector
166: . G - unprojected gradient
167: . S - step direction with which the active bounds will be estimated
168: - steplen - the step length at which the active bounds will be estimated (needs to be conservative)
170: Output Parameters:
171: . bound_tol - tolerance for for the bound estimation
172: . active_lower - index set for active variables at the lower bound
173: . active_upper - index set for active variables at the upper bound
174: . active_fixed - index set for fixed variables
175: . active - index set for all active variables
176: . inactive - complementary index set for inactive variables
178: Notes:
179: This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3.
180:
181: @*/
182: PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol,
183: IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive)
184: {
185: PetscErrorCode ierr;
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);
202:
219: VecCheckSameSize(X,1,XL,2);
220: VecCheckSameSize(X,1,XU,3);
221: VecCheckSameSize(X,1,G,4);
222: VecCheckSameSize(X,1,S,5);
223: VecCheckSameSize(X,1,W,6);
224:
225: /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
226: VecCopy(X, W);
227: VecAXPBY(W, steplen, 1.0, S);
228: TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W);
229: VecAXPBY(W, 1.0, -1.0, X);
230: VecNorm(W, NORM_2, &wnorm);
231: *bound_tol = PetscMin(*bound_tol, wnorm);
232:
233: VecGetOwnershipRange(X, &low, &high);
234: VecGetLocalSize(X, &n);
235: if (n>0){
236: VecGetArrayRead(X, &x);
237: VecGetArrayRead(XL, &xl);
238: VecGetArrayRead(XU, &xu);
239: VecGetArrayRead(G, &g);
240:
241: /* Loop over variables and categorize the indexes */
242: PetscMalloc1(n, &isl);
243: PetscMalloc1(n, &isu);
244: PetscMalloc1(n, &isf);
245: PetscMalloc1(n, &isa);
246: PetscMalloc1(n, &isi);
247: for (i=0; i<n; ++i) {
248: if (xl[i] == xu[i]) {
249: /* Fixed variables */
250: isf[n_isf]=low+i; ++n_isf;
251: isa[n_isa]=low+i; ++n_isa;
252: } else if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + *bound_tol) && (g[i] > zero)) {
253: /* Lower bounded variables */
254: isl[n_isl]=low+i; ++n_isl;
255: isa[n_isa]=low+i; ++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; ++n_isu;
259: isa[n_isa]=low+i; ++n_isa;
260: } else {
261: /* Inactive variables */
262: isi[n_isi]=low+i; ++n_isi;
263: }
264: }
265:
266: VecRestoreArrayRead(X, &x);
267: VecRestoreArrayRead(XL, &xl);
268: VecRestoreArrayRead(XU, &xu);
269: VecRestoreArrayRead(G, &g);
270: }
271:
272: /* Clear all index sets */
273: ISDestroy(active_lower);
274: ISDestroy(active_upper);
275: ISDestroy(active_fixed);
276: ISDestroy(active);
277: ISDestroy(inactive);
278:
279: /* Collect global sizes */
280: MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);
281: MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);
282: MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);
283: MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);
284: MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm);
285:
286: /* Create index set for lower bounded variables */
287: if (N_isl > 0) {
288: ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower);
289: } else {
290: PetscFree(isl);
291: }
292: /* Create index set for upper bounded variables */
293: if (N_isu > 0) {
294: ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper);
295: } else {
296: PetscFree(isu);
297: }
298: /* Create index set for fixed variables */
299: if (N_isf > 0) {
300: ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed);
301: } else {
302: PetscFree(isf);
303: }
304: /* Create index set for all actively bounded variables */
305: if (N_isa > 0) {
306: ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active);
307: } else {
308: PetscFree(isa);
309: }
310: /* Create index set for all inactive variables */
311: if (N_isi > 0) {
312: ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive);
313: } else {
314: PetscFree(isi);
315: }
317: /* Clean up and exit */
318: return(0);
319: }
321: /*@C
322: TaoBoundStep - Ensures the correct zero or adjusted step direction
323: values for active variables.
325: Input Parameters:
326: + X - solution vector
327: . XL - lower bound vector
328: . XU - upper bound vector
329: . active_lower - index set for lower bounded active variables
330: . active_upper - index set for lower bounded active variables
331: - active_fixed - index set for fixed active variables
333: Output Parameters:
334: . S - step direction to be modified
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: PetscErrorCode ierr;
339:
340: Vec step_lower, step_upper, step_fixed;
341: Vec x_lower, x_upper;
342: Vec bound_lower, bound_upper;
343:
345: /* Adjust step for variables at the estimated lower bound */
346: if (active_lower) {
347: VecGetSubVector(S, active_lower, &step_lower);
348: VecGetSubVector(X, active_lower, &x_lower);
349: VecGetSubVector(XL, active_lower, &bound_lower);
350: VecCopy(bound_lower, step_lower);
351: VecAXPY(step_lower, -1.0, x_lower);
352: VecScale(step_lower, scale);
353: VecRestoreSubVector(S, active_lower, &step_lower);
354: VecRestoreSubVector(X, active_lower, &x_lower);
355: VecRestoreSubVector(XL, active_lower, &bound_lower);
356: }
357:
358: /* Adjust step for the variables at the estimated upper bound */
359: if (active_upper) {
360: VecGetSubVector(S, active_upper, &step_upper);
361: VecGetSubVector(X, active_upper, &x_upper);
362: VecGetSubVector(XU, active_upper, &bound_upper);
363: VecCopy(bound_upper, step_upper);
364: VecAXPY(step_upper, -1.0, x_upper);
365: VecScale(step_upper, scale);
366: VecRestoreSubVector(S, active_upper, &step_upper);
367: VecRestoreSubVector(X, active_upper, &x_upper);
368: VecRestoreSubVector(XU, active_upper, &bound_upper);
369: }
370:
371: /* Zero out step for fixed variables */
372: if (active_fixed) {
373: VecGetSubVector(S, active_fixed, &step_fixed);
374: VecSet(step_fixed, 0.0);
375: VecRestoreSubVector(S, active_fixed, &step_fixed);
376: }
377: return(0);
378: }
380: /*@C
381: TaoBoundSolution - Ensures that the solution vector is snapped into the bounds.
383: Input Parameters:
384: + XL - lower bound vector
385: . XU - upper bound vector
386: . X - solution vector
387: .
388: -
390: Output Parameters:
391: . X - modified solution vector
392: @*/
393: PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout)
394: {
395: PetscErrorCode ierr;
396: PetscInt i,n,low,high,nDiff_loc=0;
397: PetscScalar *xout;
398: const PetscScalar *x,*xl,*xu;
416: VecCheckSameSize(X,1,XL,2);
417: VecCheckSameSize(X,1,XU,3);
418: VecCheckSameSize(X,1,Xout,4);
420: VecGetOwnershipRange(X,&low,&high);
421: VecGetLocalSize(X,&n);
422: if (n>0){
423: VecGetArrayRead(X, &x);
424: VecGetArrayRead(XL, &xl);
425: VecGetArrayRead(XU, &xu);
426: VecGetArray(Xout, &xout);
428: for (i=0;i<n;++i){
429: if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + bound_tol)) {
430: xout[i] = xl[i]; ++nDiff_loc;
431: } else if ((xu[i] < PETSC_INFINITY) && (x[i] >= xu[i] - bound_tol)) {
432: xout[i] = xu[i]; ++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: }