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 /*@C
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 (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)
14
15 Output Parameters:
16 . vreduced - the subvector
17
18 Notes:
19 maskvalue should usually be 0.0, unless a pointwise divide will be used.
20
21 Level: developer
22 @*/
TaoVecGetSubVec(Vec vfull,IS is,TaoSubsetType reduced_type,PetscReal maskvalue,Vec * vreduced)23 PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced)
24 {
25 PetscErrorCode ierr;
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;
34
35 PetscFunctionBegin;
36 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1);
37 PetscValidHeaderSpecific(is,IS_CLASSID,2);
38
39 ierr = VecGetSize(vfull, &nfull);CHKERRQ(ierr);
40 ierr = ISGetSize(is, &nreduced);CHKERRQ(ierr);
41
42 if (nreduced == nfull) {
43 ierr = VecDestroy(vreduced);CHKERRQ(ierr);
44 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr);
45 ierr = VecCopy(vfull,*vreduced);CHKERRQ(ierr);
46 } else {
47 switch (reduced_type) {
48 case TAO_SUBSET_SUBVEC:
49 ierr = VecGetType(vfull,&vtype);CHKERRQ(ierr);
50 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr);
51 ierr = ISGetLocalSize(is,&nreduced_local);CHKERRQ(ierr);
52 ierr = PetscObjectGetComm((PetscObject)vfull,&comm);CHKERRQ(ierr);
53 if (*vreduced) {
54 ierr = VecDestroy(vreduced);CHKERRQ(ierr);
55 }
56 ierr = VecCreate(comm,vreduced);CHKERRQ(ierr);
57 ierr = VecSetType(*vreduced,vtype);CHKERRQ(ierr);
58
59 ierr = VecSetSizes(*vreduced,nreduced_local,nreduced);CHKERRQ(ierr);
60 ierr = VecGetOwnershipRange(*vreduced,&rlow,&rhigh);CHKERRQ(ierr);
61 ierr = ISCreateStride(comm,nreduced_local,rlow,1,&ident);CHKERRQ(ierr);
62 ierr = VecScatterCreate(vfull,is,*vreduced,ident,&scatter);CHKERRQ(ierr);
63 ierr = VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
64 ierr = VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
65 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr);
66 ierr = ISDestroy(&ident);CHKERRQ(ierr);
67 break;
68
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 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr);
75 }
76
77 ierr = VecSet(*vreduced,maskvalue);CHKERRQ(ierr);
78 ierr = ISGetLocalSize(is,&nlocal);CHKERRQ(ierr);
79 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr);
80 ierr = VecGetArray(vfull,&fv);CHKERRQ(ierr);
81 ierr = VecGetArray(*vreduced,&rv);CHKERRQ(ierr);
82 ierr = ISGetIndices(is,&s);CHKERRQ(ierr);
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 ierr = ISRestoreIndices(is,&s);CHKERRQ(ierr);
88 ierr = VecRestoreArray(vfull,&fv);CHKERRQ(ierr);
89 ierr = VecRestoreArray(*vreduced,&rv);CHKERRQ(ierr);
90 break;
91 }
92 }
93 PetscFunctionReturn(0);
94 }
95
96 /*@C
97 TaoMatGetSubMat - Gets a submatrix using the IS
98
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
104
105 Output Parameters:
106 . Msub - the submatrix
107
108 Level: developer
109 @*/
TaoMatGetSubMat(Mat M,IS is,Vec v1,TaoSubsetType subset_type,Mat * Msub)110 PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
111 {
112 PetscErrorCode ierr;
113 IS iscomp;
114 PetscBool flg = PETSC_TRUE;
115
116 PetscFunctionBegin;
117 PetscValidHeaderSpecific(M,MAT_CLASSID,1);
118 PetscValidHeaderSpecific(is,IS_CLASSID,2);
119 ierr = MatDestroy(Msub);CHKERRQ(ierr);
120 switch (subset_type) {
121 case TAO_SUBSET_SUBVEC:
122 ierr = MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);CHKERRQ(ierr);
123 break;
124
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 ierr = PetscObjectOptionsBegin((PetscObject)M);CHKERRQ(ierr);
131 ierr = PetscOptionsBool("-overwrite_hessian","modify the existing hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);CHKERRQ(ierr);
132 ierr = PetscOptionsEnd();CHKERRQ(ierr);
133 if (flg) {
134 ierr = MatDuplicate(M, MAT_COPY_VALUES, Msub);CHKERRQ(ierr);
135 } else {
136 /* Act on hessian directly (default) */
137 ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr);
138 *Msub = M;
139 }
140 /* Save the diagonal to temporary vector */
141 ierr = MatGetDiagonal(*Msub,v1);CHKERRQ(ierr);
142
143 /* Zero out rows and columns */
144 ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr);
145
146 /* Use v1 instead of 0 here because of PETSc bug */
147 ierr = MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);CHKERRQ(ierr);
148
149 ierr = ISDestroy(&iscomp);CHKERRQ(ierr);
150 break;
151 case TAO_SUBSET_MATRIXFREE:
152 ierr = ISComplementVec(is,v1,&iscomp);CHKERRQ(ierr);
153 ierr = MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);CHKERRQ(ierr);
154 ierr = ISDestroy(&iscomp);CHKERRQ(ierr);
155 break;
156 }
157 PetscFunctionReturn(0);
158 }
159
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.
163
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)
172
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
180
181 Notes:
182 This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3.
183
184 Level: developer
185 @*/
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)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);
198
199 PetscFunctionBegin;
200 PetscValidHeaderSpecific(X,VEC_CLASSID,1);
201 PetscValidHeaderSpecific(XL,VEC_CLASSID,2);
202 PetscValidHeaderSpecific(XU,VEC_CLASSID,3);
203 PetscValidHeaderSpecific(G,VEC_CLASSID,4);
204 PetscValidHeaderSpecific(S,VEC_CLASSID,5);
205 PetscValidHeaderSpecific(W,VEC_CLASSID,6);
206
207 PetscValidType(X,1);
208 PetscValidType(XL,2);
209 PetscValidType(XU,3);
210 PetscValidType(G,4);
211 PetscValidType(S,5);
212 PetscValidType(W,6);
213 PetscCheckSameType(X,1,XL,2);
214 PetscCheckSameType(X,1,XU,3);
215 PetscCheckSameType(X,1,G,4);
216 PetscCheckSameType(X,1,S,5);
217 PetscCheckSameType(X,1,W,6);
218 PetscCheckSameComm(X,1,XL,2);
219 PetscCheckSameComm(X,1,XU,3);
220 PetscCheckSameComm(X,1,G,4);
221 PetscCheckSameComm(X,1,S,5);
222 PetscCheckSameComm(X,1,W,6);
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);
228
229 /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
230 ierr = VecCopy(X, W);CHKERRQ(ierr);
231 ierr = VecAXPBY(W, steplen, 1.0, S);CHKERRQ(ierr);
232 ierr = TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W);CHKERRQ(ierr);
233 ierr = VecAXPBY(W, 1.0, -1.0, X);CHKERRQ(ierr);
234 ierr = VecNorm(W, NORM_2, &wnorm);CHKERRQ(ierr);
235 *bound_tol = PetscMin(*bound_tol, wnorm);
236
237 ierr = VecGetOwnershipRange(X, &low, &high);CHKERRQ(ierr);
238 ierr = VecGetLocalSize(X, &n);CHKERRQ(ierr);
239 if (n>0){
240 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr);
241 ierr = VecGetArrayRead(XL, &xl);CHKERRQ(ierr);
242 ierr = VecGetArrayRead(XU, &xu);CHKERRQ(ierr);
243 ierr = VecGetArrayRead(G, &g);CHKERRQ(ierr);
244
245 /* Loop over variables and categorize the indexes */
246 ierr = PetscMalloc1(n, &isl);CHKERRQ(ierr);
247 ierr = PetscMalloc1(n, &isu);CHKERRQ(ierr);
248 ierr = PetscMalloc1(n, &isf);CHKERRQ(ierr);
249 ierr = PetscMalloc1(n, &isa);CHKERRQ(ierr);
250 ierr = PetscMalloc1(n, &isi);CHKERRQ(ierr);
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 }
269
270 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr);
271 ierr = VecRestoreArrayRead(XL, &xl);CHKERRQ(ierr);
272 ierr = VecRestoreArrayRead(XU, &xu);CHKERRQ(ierr);
273 ierr = VecRestoreArrayRead(G, &g);CHKERRQ(ierr);
274 }
275
276 /* Clear all index sets */
277 ierr = ISDestroy(active_lower);CHKERRQ(ierr);
278 ierr = ISDestroy(active_upper);CHKERRQ(ierr);
279 ierr = ISDestroy(active_fixed);CHKERRQ(ierr);
280 ierr = ISDestroy(active);CHKERRQ(ierr);
281 ierr = ISDestroy(inactive);CHKERRQ(ierr);
282
283 /* Collect global sizes */
284 ierr = MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);CHKERRQ(ierr);
285 ierr = MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);CHKERRQ(ierr);
286 ierr = MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);CHKERRQ(ierr);
287 ierr = MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);CHKERRQ(ierr);
288 ierr = MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm);CHKERRQ(ierr);
289
290 /* Create index set for lower bounded variables */
291 if (N_isl > 0) {
292 ierr = ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower);CHKERRQ(ierr);
293 } else {
294 ierr = PetscFree(isl);CHKERRQ(ierr);
295 }
296 /* Create index set for upper bounded variables */
297 if (N_isu > 0) {
298 ierr = ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper);CHKERRQ(ierr);
299 } else {
300 ierr = PetscFree(isu);CHKERRQ(ierr);
301 }
302 /* Create index set for fixed variables */
303 if (N_isf > 0) {
304 ierr = ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed);CHKERRQ(ierr);
305 } else {
306 ierr = PetscFree(isf);CHKERRQ(ierr);
307 }
308 /* Create index set for all actively bounded variables */
309 if (N_isa > 0) {
310 ierr = ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active);CHKERRQ(ierr);
311 } else {
312 ierr = PetscFree(isa);CHKERRQ(ierr);
313 }
314 /* Create index set for all inactive variables */
315 if (N_isi > 0) {
316 ierr = ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive);CHKERRQ(ierr);
317 } else {
318 ierr = PetscFree(isi);CHKERRQ(ierr);
319 }
320
321 /* Clean up and exit */
322 PetscFunctionReturn(0);
323 }
324
325 /*@C
326 TaoBoundStep - Ensures the correct zero or adjusted step direction
327 values for active variables.
328
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
337
338 Output Parameters:
339 . S - step direction to be modified
340
341 Level: developer
342 @*/
TaoBoundStep(Vec X,Vec XL,Vec XU,IS active_lower,IS active_upper,IS active_fixed,PetscReal scale,Vec S)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;
346
347 Vec step_lower, step_upper, step_fixed;
348 Vec x_lower, x_upper;
349 Vec bound_lower, bound_upper;
350
351 PetscFunctionBegin;
352 /* Adjust step for variables at the estimated lower bound */
353 if (active_lower) {
354 ierr = VecGetSubVector(S, active_lower, &step_lower);CHKERRQ(ierr);
355 ierr = VecGetSubVector(X, active_lower, &x_lower);CHKERRQ(ierr);
356 ierr = VecGetSubVector(XL, active_lower, &bound_lower);CHKERRQ(ierr);
357 ierr = VecCopy(bound_lower, step_lower);CHKERRQ(ierr);
358 ierr = VecAXPY(step_lower, -1.0, x_lower);CHKERRQ(ierr);
359 ierr = VecScale(step_lower, scale);CHKERRQ(ierr);
360 ierr = VecRestoreSubVector(S, active_lower, &step_lower);CHKERRQ(ierr);
361 ierr = VecRestoreSubVector(X, active_lower, &x_lower);CHKERRQ(ierr);
362 ierr = VecRestoreSubVector(XL, active_lower, &bound_lower);CHKERRQ(ierr);
363 }
364
365 /* Adjust step for the variables at the estimated upper bound */
366 if (active_upper) {
367 ierr = VecGetSubVector(S, active_upper, &step_upper);CHKERRQ(ierr);
368 ierr = VecGetSubVector(X, active_upper, &x_upper);CHKERRQ(ierr);
369 ierr = VecGetSubVector(XU, active_upper, &bound_upper);CHKERRQ(ierr);
370 ierr = VecCopy(bound_upper, step_upper);CHKERRQ(ierr);
371 ierr = VecAXPY(step_upper, -1.0, x_upper);CHKERRQ(ierr);
372 ierr = VecScale(step_upper, scale);CHKERRQ(ierr);
373 ierr = VecRestoreSubVector(S, active_upper, &step_upper);CHKERRQ(ierr);
374 ierr = VecRestoreSubVector(X, active_upper, &x_upper);CHKERRQ(ierr);
375 ierr = VecRestoreSubVector(XU, active_upper, &bound_upper);CHKERRQ(ierr);
376 }
377
378 /* Zero out step for fixed variables */
379 if (active_fixed) {
380 ierr = VecGetSubVector(S, active_fixed, &step_fixed);CHKERRQ(ierr);
381 ierr = VecSet(step_fixed, 0.0);CHKERRQ(ierr);
382 ierr = VecRestoreSubVector(S, active_fixed, &step_fixed);CHKERRQ(ierr);
383 }
384 PetscFunctionReturn(0);
385 }
386
387 /*@C
388 TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance.
389
390 Collective on Vec
391
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
397
398 Output Parameters:
399 + nDiff - total number of vector entries that have been bounded
400 - Xout - modified solution vector satisfying bounds to bound_tol
401
402 Level: developer
403
404 .seealso: TAOBNCG, TAOBNTL, TAOBNTR
405 @*/
TaoBoundSolution(Vec X,Vec XL,Vec XU,PetscReal bound_tol,PetscInt * nDiff,Vec Xout)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;
412
413 PetscFunctionBegin;
414 PetscValidHeaderSpecific(X,VEC_CLASSID,1);
415 PetscValidHeaderSpecific(XL,VEC_CLASSID,2);
416 PetscValidHeaderSpecific(XU,VEC_CLASSID,3);
417 PetscValidHeaderSpecific(Xout,VEC_CLASSID,4);
418
419 PetscValidType(X,1);
420 PetscValidType(XL,2);
421 PetscValidType(XU,3);
422 PetscValidType(Xout,4);
423 PetscCheckSameType(X,1,XL,2);
424 PetscCheckSameType(X,1,XU,3);
425 PetscCheckSameType(X,1,Xout,4);
426 PetscCheckSameComm(X,1,XL,2);
427 PetscCheckSameComm(X,1,XU,3);
428 PetscCheckSameComm(X,1,Xout,4);
429 VecCheckSameSize(X,1,XL,2);
430 VecCheckSameSize(X,1,XU,3);
431 VecCheckSameSize(X,1,Xout,4);
432
433 ierr = VecGetOwnershipRange(X,&low,&high);CHKERRQ(ierr);
434 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
435 if (n>0){
436 ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr);
437 ierr = VecGetArrayRead(XL, &xl);CHKERRQ(ierr);
438 ierr = VecGetArrayRead(XU, &xu);CHKERRQ(ierr);
439 ierr = VecGetArray(Xout, &xout);CHKERRQ(ierr);
440
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 }
448
449 ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr);
450 ierr = VecRestoreArrayRead(XL, &xl);CHKERRQ(ierr);
451 ierr = VecRestoreArrayRead(XU, &xu);CHKERRQ(ierr);
452 ierr = VecRestoreArray(Xout, &xout);CHKERRQ(ierr);
453 }
454 ierr = MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X));CHKERRQ(ierr);
455 PetscFunctionReturn(0);
456 }
457