1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 #include <petscdm.h>
3
4 /*@C
5 SNESVISetComputeVariableBounds - Sets a function that is called to compute the variable bounds
6
7 Input parameter:
8 + snes - the SNES context
9 - compute - computes the bounds
10
11 Level: advanced
12
13 .seealso: SNESVISetVariableBounds()
14
15 @*/
SNESVISetComputeVariableBounds(SNES snes,PetscErrorCode (* compute)(SNES,Vec,Vec))16 PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES,Vec,Vec))
17 {
18 PetscErrorCode ierr,(*f)(SNES,PetscErrorCode (*)(SNES,Vec,Vec));
19
20 PetscFunctionBegin;
21 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
22 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",&f);CHKERRQ(ierr);
23 if (!f) {
24 ierr = SNESVISetComputeVariableBounds_VI(snes,compute);CHKERRQ(ierr);
25 } else {
26 ierr = PetscUseMethod(snes,"SNESVISetComputeVariableBounds_C",(SNES,PetscErrorCode (*)(SNES,Vec,Vec)),(snes,compute));CHKERRQ(ierr);
27 }
28 PetscFunctionReturn(0);
29 }
30
SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)31 PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes,SNESVIComputeVariableBoundsFunction compute)
32 {
33 PetscFunctionBegin;
34 snes->ops->computevariablebounds = compute;
35 PetscFunctionReturn(0);
36 }
37
38 /* --------------------------------------------------------------------------------------------------------*/
39
SNESVIMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void * dummy)40 PetscErrorCode SNESVIMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
41 {
42 PetscErrorCode ierr;
43 Vec X, F, Finactive;
44 IS isactive;
45 PetscViewer viewer = (PetscViewer) dummy;
46
47 PetscFunctionBegin;
48 ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr);
49 ierr = SNESGetSolution(snes,&X);CHKERRQ(ierr);
50 ierr = SNESVIGetActiveSetIS(snes,X,F,&isactive);CHKERRQ(ierr);
51 ierr = VecDuplicate(F,&Finactive);CHKERRQ(ierr);
52 ierr = VecCopy(F,Finactive);CHKERRQ(ierr);
53 ierr = VecISSet(Finactive,isactive,0.0);CHKERRQ(ierr);
54 ierr = ISDestroy(&isactive);CHKERRQ(ierr);
55 ierr = VecView(Finactive,viewer);CHKERRQ(ierr);
56 ierr = VecDestroy(&Finactive);CHKERRQ(ierr);
57 PetscFunctionReturn(0);
58 }
59
SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void * dummy)60 PetscErrorCode SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
61 {
62 PetscErrorCode ierr;
63 PetscViewer viewer = (PetscViewer) dummy;
64 const PetscScalar *x,*xl,*xu,*f;
65 PetscInt i,n,act[2] = {0,0},fact[2],N;
66 /* Number of components that actually hit the bounds (c.f. active variables) */
67 PetscInt act_bound[2] = {0,0},fact_bound[2];
68 PetscReal rnorm,fnorm,zerotolerance = snes->vizerotolerance;
69 double tmp;
70
71 PetscFunctionBegin;
72 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
73 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr);
74 ierr = VecGetSize(snes->vec_sol,&N);CHKERRQ(ierr);
75 ierr = VecGetArrayRead(snes->xl,&xl);CHKERRQ(ierr);
76 ierr = VecGetArrayRead(snes->xu,&xu);CHKERRQ(ierr);
77 ierr = VecGetArrayRead(snes->vec_sol,&x);CHKERRQ(ierr);
78 ierr = VecGetArrayRead(snes->vec_func,&f);CHKERRQ(ierr);
79
80 rnorm = 0.0;
81 for (i=0; i<n; i++) {
82 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]);
83 else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
84 else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
85 else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_PLIB,"Can never get here");
86 }
87
88 for (i=0; i<n; i++) {
89 if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
90 else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
91 }
92 ierr = VecRestoreArrayRead(snes->vec_func,&f);CHKERRQ(ierr);
93 ierr = VecRestoreArrayRead(snes->xl,&xl);CHKERRQ(ierr);
94 ierr = VecRestoreArrayRead(snes->xu,&xu);CHKERRQ(ierr);
95 ierr = VecRestoreArrayRead(snes->vec_sol,&x);CHKERRQ(ierr);
96 ierr = MPIU_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
97 ierr = MPIU_Allreduce(act,fact,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
98 ierr = MPIU_Allreduce(act_bound,fact_bound,2,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
99 fnorm = PetscSqrtReal(fnorm);
100
101 ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
102 if (snes->ntruebounds) tmp = ((double)(fact[0]+fact[1]))/((double)snes->ntruebounds);
103 else tmp = 0.0;
104 ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES VI Function norm %g Active lower constraints %D/%D upper constraints %D/%D Percent of total %g Percent of bounded %g\n",its,(double)fnorm,fact[0],fact_bound[0],fact[1],fact_bound[1],((double)(fact[0]+fact[1]))/((double)N),tmp);CHKERRQ(ierr);
105
106 ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
107 PetscFunctionReturn(0);
108 }
109
110 /*
111 Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
112 || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that
113 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
114 for this trick. One assumes that the probability that W is in the null space of J is very, very small.
115 */
SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool * ismin)116 PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin)
117 {
118 PetscReal a1;
119 PetscErrorCode ierr;
120 PetscBool hastranspose;
121
122 PetscFunctionBegin;
123 *ismin = PETSC_FALSE;
124 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr);
125 if (hastranspose) {
126 /* Compute || J^T F|| */
127 ierr = MatMultTranspose(A,F,W);CHKERRQ(ierr);
128 ierr = VecNorm(W,NORM_2,&a1);CHKERRQ(ierr);
129 ierr = PetscInfo1(snes,"|| J^T F|| %g near zero implies found a local minimum\n",(double)(a1/fnorm));CHKERRQ(ierr);
130 if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE;
131 } else {
132 Vec work;
133 PetscScalar result;
134 PetscReal wnorm;
135
136 ierr = VecSetRandom(W,NULL);CHKERRQ(ierr);
137 ierr = VecNorm(W,NORM_2,&wnorm);CHKERRQ(ierr);
138 ierr = VecDuplicate(W,&work);CHKERRQ(ierr);
139 ierr = MatMult(A,W,work);CHKERRQ(ierr);
140 ierr = VecDot(F,work,&result);CHKERRQ(ierr);
141 ierr = VecDestroy(&work);CHKERRQ(ierr);
142 a1 = PetscAbsScalar(result)/(fnorm*wnorm);
143 ierr = PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n",(double)a1);CHKERRQ(ierr);
144 if (a1 < 1.e-4) *ismin = PETSC_TRUE;
145 }
146 PetscFunctionReturn(0);
147 }
148
149 /*
150 Checks if J^T(F - J*X) = 0
151 */
SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)152 PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2)
153 {
154 PetscReal a1,a2;
155 PetscErrorCode ierr;
156 PetscBool hastranspose;
157
158 PetscFunctionBegin;
159 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr);
160 if (hastranspose) {
161 ierr = MatMult(A,X,W1);CHKERRQ(ierr);
162 ierr = VecAXPY(W1,-1.0,F);CHKERRQ(ierr);
163
164 /* Compute || J^T W|| */
165 ierr = MatMultTranspose(A,W1,W2);CHKERRQ(ierr);
166 ierr = VecNorm(W1,NORM_2,&a1);CHKERRQ(ierr);
167 ierr = VecNorm(W2,NORM_2,&a2);CHKERRQ(ierr);
168 if (a1 != 0.0) {
169 ierr = PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %g near zero implies inconsistent rhs\n",(double)(a2/a1));CHKERRQ(ierr);
170 }
171 }
172 PetscFunctionReturn(0);
173 }
174
175 /*
176 SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
177
178 Notes:
179 The convergence criterion currently implemented is
180 merit < abstol
181 merit < rtol*merit_initial
182 */
SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason * reason,void * dummy)183 PetscErrorCode SNESConvergedDefault_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
184 {
185 PetscErrorCode ierr;
186
187 PetscFunctionBegin;
188 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
189 PetscValidPointer(reason,6);
190
191 *reason = SNES_CONVERGED_ITERATING;
192
193 if (!it) {
194 /* set parameter for default relative tolerance convergence test */
195 snes->ttol = fnorm*snes->rtol;
196 }
197 if (fnorm != fnorm) {
198 ierr = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr);
199 *reason = SNES_DIVERGED_FNORM_NAN;
200 } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
201 ierr = PetscInfo2(snes,"Converged due to function norm %g < %g\n",(double)fnorm,(double)snes->abstol);CHKERRQ(ierr);
202 *reason = SNES_CONVERGED_FNORM_ABS;
203 } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
204 ierr = PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);CHKERRQ(ierr);
205 *reason = SNES_DIVERGED_FUNCTION_COUNT;
206 }
207
208 if (it && !*reason) {
209 if (fnorm < snes->ttol) {
210 ierr = PetscInfo2(snes,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)snes->ttol);CHKERRQ(ierr);
211 *reason = SNES_CONVERGED_FNORM_RELATIVE;
212 }
213 }
214 PetscFunctionReturn(0);
215 }
216
217
218 /* -------------------------------------------------------------------------- */
219 /*
220 SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
221
222 Input Parameters:
223 . SNES - nonlinear solver context
224
225 Output Parameters:
226 . X - Bound projected X
227
228 */
229
SNESVIProjectOntoBounds(SNES snes,Vec X)230 PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X)
231 {
232 PetscErrorCode ierr;
233 const PetscScalar *xl,*xu;
234 PetscScalar *x;
235 PetscInt i,n;
236
237 PetscFunctionBegin;
238 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
239 ierr = VecGetArray(X,&x);CHKERRQ(ierr);
240 ierr = VecGetArrayRead(snes->xl,&xl);CHKERRQ(ierr);
241 ierr = VecGetArrayRead(snes->xu,&xu);CHKERRQ(ierr);
242
243 for (i = 0; i<n; i++) {
244 if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
245 else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
246 }
247 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
248 ierr = VecRestoreArrayRead(snes->xl,&xl);CHKERRQ(ierr);
249 ierr = VecRestoreArrayRead(snes->xu,&xu);CHKERRQ(ierr);
250 PetscFunctionReturn(0);
251 }
252
253
254 /*
255 SNESVIGetActiveSetIndices - Gets the global indices for the active set variables
256
257 Input parameter:
258 . snes - the SNES context
259 . X - the snes solution vector
260 . F - the nonlinear function vector
261
262 Output parameter:
263 . ISact - active set index set
264 */
SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS * ISact)265 PetscErrorCode SNESVIGetActiveSetIS(SNES snes,Vec X,Vec F,IS *ISact)
266 {
267 PetscErrorCode ierr;
268 Vec Xl=snes->xl,Xu=snes->xu;
269 const PetscScalar *x,*f,*xl,*xu;
270 PetscInt *idx_act,i,nlocal,nloc_isact=0,ilow,ihigh,i1=0;
271 PetscReal zerotolerance = snes->vizerotolerance;
272
273 PetscFunctionBegin;
274 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr);
275 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr);
276 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
277 ierr = VecGetArrayRead(Xl,&xl);CHKERRQ(ierr);
278 ierr = VecGetArrayRead(Xu,&xu);CHKERRQ(ierr);
279 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr);
280 /* Compute active set size */
281 for (i=0; i < nlocal;i++) {
282 if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) nloc_isact++;
283 }
284
285 ierr = PetscMalloc1(nloc_isact,&idx_act);CHKERRQ(ierr);
286
287 /* Set active set indices */
288 for (i=0; i < nlocal; i++) {
289 if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) idx_act[i1++] = ilow+i;
290 }
291
292 /* Create active set IS */
293 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),nloc_isact,idx_act,PETSC_OWN_POINTER,ISact);CHKERRQ(ierr);
294
295 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
296 ierr = VecRestoreArrayRead(Xl,&xl);CHKERRQ(ierr);
297 ierr = VecRestoreArrayRead(Xu,&xu);CHKERRQ(ierr);
298 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr);
299 PetscFunctionReturn(0);
300 }
301
SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS * ISact,IS * ISinact)302 PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
303 {
304 PetscErrorCode ierr;
305 PetscInt rstart,rend;
306
307 PetscFunctionBegin;
308 ierr = SNESVIGetActiveSetIS(snes,X,F,ISact);CHKERRQ(ierr);
309 ierr = VecGetOwnershipRange(X,&rstart,&rend);CHKERRQ(ierr);
310 ierr = ISComplement(*ISact,rstart,rend,ISinact);CHKERRQ(ierr);
311 PetscFunctionReturn(0);
312 }
313
SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X,PetscReal * fnorm)314 PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X, PetscReal *fnorm)
315 {
316 PetscErrorCode ierr;
317 const PetscScalar *x,*xl,*xu,*f;
318 PetscInt i,n;
319 PetscReal rnorm,zerotolerance = snes->vizerotolerance;
320
321 PetscFunctionBegin;
322 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
323 ierr = VecGetArrayRead(snes->xl,&xl);CHKERRQ(ierr);
324 ierr = VecGetArrayRead(snes->xu,&xu);CHKERRQ(ierr);
325 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
326 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr);
327 rnorm = 0.0;
328 for (i=0; i<n; i++) {
329 if (((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) rnorm += PetscRealPart(PetscConj(f[i])*f[i]);
330 }
331 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr);
332 ierr = VecRestoreArrayRead(snes->xl,&xl);CHKERRQ(ierr);
333 ierr = VecRestoreArrayRead(snes->xu,&xu);CHKERRQ(ierr);
334 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
335 ierr = MPIU_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
336 *fnorm = PetscSqrtReal(*fnorm);
337 PetscFunctionReturn(0);
338 }
339
SNESVIDMComputeVariableBounds(SNES snes,Vec xl,Vec xu)340 PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes,Vec xl, Vec xu)
341 {
342 PetscErrorCode ierr;
343
344 PetscFunctionBegin;
345 ierr = DMComputeVariableBounds(snes->dm, xl, xu);CHKERRQ(ierr);
346 PetscFunctionReturn(0);
347 }
348
349
350 /* -------------------------------------------------------------------------- */
351 /*
352 SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
353 of the SNESVI nonlinear solver.
354
355 Input Parameter:
356 . snes - the SNES context
357
358 Application Interface Routine: SNESSetUp()
359
360 Notes:
361 For basic use of the SNES solvers, the user need not explicitly call
362 SNESSetUp(), since these actions will automatically occur during
363 the call to SNESSolve().
364 */
SNESSetUp_VI(SNES snes)365 PetscErrorCode SNESSetUp_VI(SNES snes)
366 {
367 PetscErrorCode ierr;
368 PetscInt i_start[3],i_end[3];
369
370 PetscFunctionBegin;
371 ierr = SNESSetWorkVecs(snes,1);CHKERRQ(ierr);
372 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
373
374 if (!snes->ops->computevariablebounds && snes->dm) {
375 PetscBool flag;
376 ierr = DMHasVariableBounds(snes->dm, &flag);CHKERRQ(ierr);
377 if (flag) {
378 snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
379 }
380 }
381 if (!snes->usersetbounds) {
382 if (snes->ops->computevariablebounds) {
383 if (!snes->xl) {ierr = VecDuplicate(snes->vec_sol,&snes->xl);CHKERRQ(ierr);}
384 if (!snes->xu) {ierr = VecDuplicate(snes->vec_sol,&snes->xu);CHKERRQ(ierr);}
385 ierr = (*snes->ops->computevariablebounds)(snes,snes->xl,snes->xu);CHKERRQ(ierr);
386 } else if (!snes->xl && !snes->xu) {
387 /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
388 ierr = VecDuplicate(snes->vec_sol, &snes->xl);CHKERRQ(ierr);
389 ierr = VecSet(snes->xl,PETSC_NINFINITY);CHKERRQ(ierr);
390 ierr = VecDuplicate(snes->vec_sol, &snes->xu);CHKERRQ(ierr);
391 ierr = VecSet(snes->xu,PETSC_INFINITY);CHKERRQ(ierr);
392 } else {
393 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
394 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr);
395 ierr = VecGetOwnershipRange(snes->xl,i_start+1,i_end+1);CHKERRQ(ierr);
396 ierr = VecGetOwnershipRange(snes->xu,i_start+2,i_end+2);CHKERRQ(ierr);
397 if ((i_start[0] != i_start[1]) || (i_start[0] != i_start[2]) || (i_end[0] != i_end[1]) || (i_end[0] != i_end[2]))
398 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Distribution of lower bound, upper bound and the solution vector should be identical across all the processors.");
399 }
400 }
401 PetscFunctionReturn(0);
402 }
403 /* -------------------------------------------------------------------------- */
SNESReset_VI(SNES snes)404 PetscErrorCode SNESReset_VI(SNES snes)
405 {
406 PetscErrorCode ierr;
407
408 PetscFunctionBegin;
409 ierr = VecDestroy(&snes->xl);CHKERRQ(ierr);
410 ierr = VecDestroy(&snes->xu);CHKERRQ(ierr);
411 snes->usersetbounds = PETSC_FALSE;
412 PetscFunctionReturn(0);
413 }
414
415 /*
416 SNESDestroy_VI - Destroys the private SNES_VI context that was created
417 with SNESCreate_VI().
418
419 Input Parameter:
420 . snes - the SNES context
421
422 Application Interface Routine: SNESDestroy()
423 */
SNESDestroy_VI(SNES snes)424 PetscErrorCode SNESDestroy_VI(SNES snes)
425 {
426 PetscErrorCode ierr;
427
428 PetscFunctionBegin;
429 ierr = PetscFree(snes->data);CHKERRQ(ierr);
430
431 /* clear composed functions */
432 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSet_C",NULL);CHKERRQ(ierr);
433 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESLineSearchSetDefaultMonitor_C",NULL);CHKERRQ(ierr);
434 PetscFunctionReturn(0);
435 }
436
437 /*@
438 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu.
439
440 Input Parameters:
441 + snes - the SNES context.
442 . xl - lower bound.
443 - xu - upper bound.
444
445 Notes:
446 If this routine is not called then the lower and upper bounds are set to
447 PETSC_NINFINITY and PETSC_INFINITY respectively during SNESSetUp().
448
449 Level: advanced
450
451 @*/
SNESVISetVariableBounds(SNES snes,Vec xl,Vec xu)452 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
453 {
454 PetscErrorCode ierr,(*f)(SNES,Vec,Vec);
455
456 PetscFunctionBegin;
457 PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
458 PetscValidHeaderSpecific(xl,VEC_CLASSID,2);
459 PetscValidHeaderSpecific(xu,VEC_CLASSID,3);
460 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESVISetVariableBounds_C",&f);CHKERRQ(ierr);
461 if (!f) {
462 ierr = SNESVISetVariableBounds_VI(snes, xl, xu);CHKERRQ(ierr);
463 } else {
464 ierr = PetscUseMethod(snes,"SNESVISetVariableBounds_C",(SNES,Vec,Vec),(snes,xl,xu));CHKERRQ(ierr);
465 }
466 snes->usersetbounds = PETSC_TRUE;
467 PetscFunctionReturn(0);
468 }
469
SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)470 PetscErrorCode SNESVISetVariableBounds_VI(SNES snes,Vec xl,Vec xu)
471 {
472 PetscErrorCode ierr;
473 const PetscScalar *xxl,*xxu;
474 PetscInt i,n, cnt = 0;
475
476 PetscFunctionBegin;
477 ierr = SNESGetFunction(snes,&snes->vec_func,NULL,NULL);CHKERRQ(ierr);
478 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
479 {
480 PetscInt xlN,xuN,N;
481 ierr = VecGetSize(xl,&xlN);CHKERRQ(ierr);
482 ierr = VecGetSize(xu,&xuN);CHKERRQ(ierr);
483 ierr = VecGetSize(snes->vec_func,&N);CHKERRQ(ierr);
484 if (xlN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xlN,N);
485 if (xuN != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xuN,N);
486 }
487 ierr = PetscObjectReference((PetscObject)xl);CHKERRQ(ierr);
488 ierr = PetscObjectReference((PetscObject)xu);CHKERRQ(ierr);
489 ierr = VecDestroy(&snes->xl);CHKERRQ(ierr);
490 ierr = VecDestroy(&snes->xu);CHKERRQ(ierr);
491 snes->xl = xl;
492 snes->xu = xu;
493 ierr = VecGetLocalSize(xl,&n);CHKERRQ(ierr);
494 ierr = VecGetArrayRead(xl,&xxl);CHKERRQ(ierr);
495 ierr = VecGetArrayRead(xu,&xxu);CHKERRQ(ierr);
496 for (i=0; i<n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
497
498 ierr = MPIU_Allreduce(&cnt,&snes->ntruebounds,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
499 ierr = VecRestoreArrayRead(xl,&xxl);CHKERRQ(ierr);
500 ierr = VecRestoreArrayRead(xu,&xxu);CHKERRQ(ierr);
501 PetscFunctionReturn(0);
502 }
503
SNESSetFromOptions_VI(PetscOptionItems * PetscOptionsObject,SNES snes)504 PetscErrorCode SNESSetFromOptions_VI(PetscOptionItems *PetscOptionsObject,SNES snes)
505 {
506 PetscErrorCode ierr;
507 PetscBool flg = PETSC_FALSE;
508
509 PetscFunctionBegin;
510 ierr = PetscOptionsHead(PetscOptionsObject,"SNES VI options");CHKERRQ(ierr);
511 ierr = PetscOptionsReal("-snes_vi_zero_tolerance","Tolerance for considering x[] value to be on a bound","None",snes->vizerotolerance,&snes->vizerotolerance,NULL);CHKERRQ(ierr);
512 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","SNESMonitorResidual",flg,&flg,NULL);CHKERRQ(ierr);
513 if (flg) {
514 ierr = SNESMonitorSet(snes,SNESMonitorVI,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)),NULL);CHKERRQ(ierr);
515 }
516 flg = PETSC_FALSE;
517 ierr = PetscOptionsBool("-snes_vi_monitor_residual","Monitor residual all non-active variables; using zero for active constraints","SNESMonitorVIResidual",flg,&flg,NULL);CHKERRQ(ierr);
518 if (flg) {
519 ierr = SNESMonitorSet(snes,SNESVIMonitorResidual,PETSC_VIEWER_DRAW_(PetscObjectComm((PetscObject)snes)),NULL);CHKERRQ(ierr);
520 }
521 ierr = PetscOptionsTail();CHKERRQ(ierr);
522 PetscFunctionReturn(0);
523 }
524