xref: /petsc/src/tao/util/tao_util.c (revision 7f296bb328fcd4c99f2da7bfe8ba7ed8a4ebceee)
1 #include <petsc/private/petscimpl.h>
2 #include <petsctao.h> /*I "petsctao.h" I*/
3 #include <petscsys.h>
4 
5 static inline PetscReal Fischer(PetscReal a, PetscReal b)
6 {
7   /* Method suggested by Bob Vanderbei */
8   if (a + b <= 0) return PetscSqrtReal(a * a + b * b) - (a + b);
9   return -2.0 * a * b / (PetscSqrtReal(a * a + b * b) + (a + b));
10 }
11 
12 /*@
13   VecFischer - Evaluates the Fischer-Burmeister function for complementarity
14   problems.
15 
16   Logically Collective
17 
18   Input Parameters:
19 + X - current point
20 . F - function evaluated at x
21 . L - lower bounds
22 - U - upper bounds
23 
24   Output Parameter:
25 . FB - The Fischer-Burmeister function vector
26 
27   Level: developer
28 
29   Notes:
30   The Fischer-Burmeister function is defined as
31 $        phi(a,b) := sqrt(a*a + b*b) - a - b
32   and is used reformulate a complementarity problem as a semismooth
33   system of equations.
34 
35   The result of this function is done by cases\:
36 +  l[i] == -infinity, u[i] == infinity  -- fb[i] = -f[i]
37 .  l[i] == -infinity, u[i] finite       -- fb[i] = phi(u[i]-x[i], -f[i])
38 .  l[i] finite,       u[i] == infinity  -- fb[i] = phi(x[i]-l[i],  f[i])
39 .  l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u]))
40 -  otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
41 
42 .seealso: `Vec`, `VecSFischer()`, `MatDFischer()`, `MatDSFischer()`
43 @*/
44 PetscErrorCode VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FB)
45 {
46   const PetscScalar *x, *f, *l, *u;
47   PetscScalar       *fb;
48   PetscReal          xval, fval, lval, uval;
49   PetscInt           low[5], high[5], n, i;
50 
51   PetscFunctionBegin;
52   PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
53   PetscValidHeaderSpecific(F, VEC_CLASSID, 2);
54   if (L) PetscValidHeaderSpecific(L, VEC_CLASSID, 3);
55   if (U) PetscValidHeaderSpecific(U, VEC_CLASSID, 4);
56   PetscValidHeaderSpecific(FB, VEC_CLASSID, 5);
57 
58   if (!L && !U) {
59     PetscCall(VecAXPBY(FB, -1.0, 0.0, F));
60     PetscFunctionReturn(PETSC_SUCCESS);
61   }
62 
63   PetscCall(VecGetOwnershipRange(X, low, high));
64   PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
65   PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
66   PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
67   PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));
68 
69   for (i = 1; i < 4; ++i) PetscCheck(low[0] == low[i] && high[0] == high[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Vectors must be identically loaded over processors");
70 
71   PetscCall(VecGetArrayRead(X, &x));
72   PetscCall(VecGetArrayRead(F, &f));
73   PetscCall(VecGetArrayRead(L, &l));
74   PetscCall(VecGetArrayRead(U, &u));
75   PetscCall(VecGetArray(FB, &fb));
76 
77   PetscCall(VecGetLocalSize(X, &n));
78 
79   for (i = 0; i < n; ++i) {
80     xval = PetscRealPart(x[i]);
81     fval = PetscRealPart(f[i]);
82     lval = PetscRealPart(l[i]);
83     uval = PetscRealPart(u[i]);
84 
85     if (lval <= -PETSC_INFINITY && uval >= PETSC_INFINITY) {
86       fb[i] = -fval;
87     } else if (lval <= -PETSC_INFINITY) {
88       fb[i] = -Fischer(uval - xval, -fval);
89     } else if (uval >= PETSC_INFINITY) {
90       fb[i] = Fischer(xval - lval, fval);
91     } else if (lval == uval) {
92       fb[i] = lval - xval;
93     } else {
94       fval  = Fischer(uval - xval, -fval);
95       fb[i] = Fischer(xval - lval, fval);
96     }
97   }
98 
99   PetscCall(VecRestoreArrayRead(X, &x));
100   PetscCall(VecRestoreArrayRead(F, &f));
101   PetscCall(VecRestoreArrayRead(L, &l));
102   PetscCall(VecRestoreArrayRead(U, &u));
103   PetscCall(VecRestoreArray(FB, &fb));
104   PetscFunctionReturn(PETSC_SUCCESS);
105 }
106 
107 static inline PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c)
108 {
109   /* Method suggested by Bob Vanderbei */
110   if (a + b <= 0) return PetscSqrtReal(a * a + b * b + 2.0 * c * c) - (a + b);
111   return 2.0 * (c * c - a * b) / (PetscSqrtReal(a * a + b * b + 2.0 * c * c) + (a + b));
112 }
113 
114 /*@
115   VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for
116   complementarity problems.
117 
118   Logically Collective
119 
120   Input Parameters:
121 + X  - current point
122 . F  - function evaluated at x
123 . L  - lower bounds
124 . U  - upper bounds
125 - mu - smoothing parameter
126 
127   Output Parameter:
128 . FB - The Smoothed Fischer-Burmeister function vector
129 
130   Notes:
131   The Smoothed Fischer-Burmeister function is defined as
132 $        phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b
133   and is used reformulate a complementarity problem as a semismooth
134   system of equations.
135 
136   The result of this function is done by cases\:
137 +  l[i] == -infinity, u[i] == infinity  -- fb[i] = -f[i] - 2*mu*x[i]
138 .  l[i] == -infinity, u[i] finite       -- fb[i] = phi(u[i]-x[i], -f[i], mu)
139 .  l[i] finite,       u[i] == infinity  -- fb[i] = phi(x[i]-l[i],  f[i], mu)
140 .  l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
141 -  otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
142 
143   Level: developer
144 
145 .seealso: `Vec`, `VecFischer()`, `MatDFischer()`, `MatDSFischer()`
146 @*/
147 PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB)
148 {
149   const PetscScalar *x, *f, *l, *u;
150   PetscScalar       *fb;
151   PetscReal          xval, fval, lval, uval;
152   PetscInt           low[5], high[5], n, i;
153 
154   PetscFunctionBegin;
155   PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
156   PetscValidHeaderSpecific(F, VEC_CLASSID, 2);
157   PetscValidHeaderSpecific(L, VEC_CLASSID, 3);
158   PetscValidHeaderSpecific(U, VEC_CLASSID, 4);
159   PetscValidHeaderSpecific(FB, VEC_CLASSID, 6);
160 
161   PetscCall(VecGetOwnershipRange(X, low, high));
162   PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
163   PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
164   PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
165   PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));
166 
167   for (i = 1; i < 4; ++i) PetscCheck(low[0] == low[i] && high[0] == high[i], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Vectors must be identically loaded over processors");
168 
169   PetscCall(VecGetArrayRead(X, &x));
170   PetscCall(VecGetArrayRead(F, &f));
171   PetscCall(VecGetArrayRead(L, &l));
172   PetscCall(VecGetArrayRead(U, &u));
173   PetscCall(VecGetArray(FB, &fb));
174 
175   PetscCall(VecGetLocalSize(X, &n));
176 
177   for (i = 0; i < n; ++i) {
178     xval = PetscRealPart(*x++);
179     fval = PetscRealPart(*f++);
180     lval = PetscRealPart(*l++);
181     uval = PetscRealPart(*u++);
182 
183     if ((lval <= -PETSC_INFINITY) && (uval >= PETSC_INFINITY)) {
184       (*fb++) = -fval - mu * xval;
185     } else if (lval <= -PETSC_INFINITY) {
186       (*fb++) = -SFischer(uval - xval, -fval, mu);
187     } else if (uval >= PETSC_INFINITY) {
188       (*fb++) = SFischer(xval - lval, fval, mu);
189     } else if (lval == uval) {
190       (*fb++) = lval - xval;
191     } else {
192       fval    = SFischer(uval - xval, -fval, mu);
193       (*fb++) = SFischer(xval - lval, fval, mu);
194     }
195   }
196   x -= n;
197   f -= n;
198   l -= n;
199   u -= n;
200   fb -= n;
201 
202   PetscCall(VecRestoreArrayRead(X, &x));
203   PetscCall(VecRestoreArrayRead(F, &f));
204   PetscCall(VecRestoreArrayRead(L, &l));
205   PetscCall(VecRestoreArrayRead(U, &u));
206   PetscCall(VecRestoreArray(FB, &fb));
207   PetscFunctionReturn(PETSC_SUCCESS);
208 }
209 
210 static inline PetscReal fischnorm(PetscReal a, PetscReal b)
211 {
212   return PetscSqrtReal(a * a + b * b);
213 }
214 
215 static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
216 {
217   return PetscSqrtReal(a * a + b * b + 2.0 * c * c);
218 }
219 
220 /*@
221   MatDFischer - Calculates an element of the B-subdifferential of the
222   Fischer-Burmeister function for complementarity problems.
223 
224   Collective
225 
226   Input Parameters:
227 + jac - the jacobian of `f` at `X`
228 . X   - current point
229 . Con - constraints function evaluated at `X`
230 . XL  - lower bounds
231 . XU  - upper bounds
232 . T1  - work vector
233 - T2  - work vector
234 
235   Output Parameters:
236 + Da - diagonal perturbation component of the result
237 - Db - row scaling component of the result
238 
239   Level: developer
240 
241 .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDSFischer()`
242 @*/
243 PetscErrorCode MatDFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db)
244 {
245   PetscInt           i, nn;
246   const PetscScalar *x, *f, *l, *u, *t2;
247   PetscScalar       *da, *db, *t1;
248   PetscReal          ai, bi, ci, di, ei;
249 
250   PetscFunctionBegin;
251   PetscCall(VecGetLocalSize(X, &nn));
252   PetscCall(VecGetArrayRead(X, &x));
253   PetscCall(VecGetArrayRead(Con, &f));
254   PetscCall(VecGetArrayRead(XL, &l));
255   PetscCall(VecGetArrayRead(XU, &u));
256   PetscCall(VecGetArray(Da, &da));
257   PetscCall(VecGetArray(Db, &db));
258   PetscCall(VecGetArray(T1, &t1));
259   PetscCall(VecGetArrayRead(T2, &t2));
260 
261   for (i = 0; i < nn; i++) {
262     da[i] = 0.0;
263     db[i] = 0.0;
264     t1[i] = 0.0;
265 
266     if (PetscAbsScalar(f[i]) <= PETSC_MACHINE_EPSILON) {
267       if (PetscRealPart(l[i]) > PETSC_NINFINITY && PetscAbsScalar(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) {
268         t1[i] = 1.0;
269         da[i] = 1.0;
270       }
271 
272       if (PetscRealPart(u[i]) < PETSC_INFINITY && PetscAbsScalar(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) {
273         t1[i] = 1.0;
274         db[i] = 1.0;
275       }
276     }
277   }
278 
279   PetscCall(VecRestoreArray(T1, &t1));
280   PetscCall(VecRestoreArrayRead(T2, &t2));
281   PetscCall(MatMult(jac, T1, T2));
282   PetscCall(VecGetArrayRead(T2, &t2));
283 
284   for (i = 0; i < nn; i++) {
285     if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
286       da[i] = 0.0;
287       db[i] = -1.0;
288     } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
289       if (PetscRealPart(db[i]) >= 1) {
290         ai = fischnorm(1.0, PetscRealPart(t2[i]));
291 
292         da[i] = -1.0 / ai - 1.0;
293         db[i] = -t2[i] / ai - 1.0;
294       } else {
295         bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
296         ai = fischnorm(bi, PetscRealPart(f[i]));
297         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
298 
299         da[i] = bi / ai - 1.0;
300         db[i] = -f[i] / ai - 1.0;
301       }
302     } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
303       if (PetscRealPart(da[i]) >= 1) {
304         ai = fischnorm(1.0, PetscRealPart(t2[i]));
305 
306         da[i] = 1.0 / ai - 1.0;
307         db[i] = t2[i] / ai - 1.0;
308       } else {
309         bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
310         ai = fischnorm(bi, PetscRealPart(f[i]));
311         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
312 
313         da[i] = bi / ai - 1.0;
314         db[i] = f[i] / ai - 1.0;
315       }
316     } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
317       da[i] = -1.0;
318       db[i] = 0.0;
319     } else {
320       if (PetscRealPart(db[i]) >= 1) {
321         ai = fischnorm(1.0, PetscRealPart(t2[i]));
322 
323         ci = 1.0 / ai + 1.0;
324         di = PetscRealPart(t2[i]) / ai + 1.0;
325       } else {
326         bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
327         ai = fischnorm(bi, PetscRealPart(f[i]));
328         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
329 
330         ci = bi / ai + 1.0;
331         di = PetscRealPart(f[i]) / ai + 1.0;
332       }
333 
334       if (PetscRealPart(da[i]) >= 1) {
335         bi = ci + di * PetscRealPart(t2[i]);
336         ai = fischnorm(1.0, bi);
337 
338         bi = bi / ai - 1.0;
339         ai = 1.0 / ai - 1.0;
340       } else {
341         ei = Fischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]));
342         ai = fischnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei);
343         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
344 
345         bi = ei / ai - 1.0;
346         ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
347       }
348 
349       da[i] = ai + bi * ci;
350       db[i] = bi * di;
351     }
352   }
353 
354   PetscCall(VecRestoreArray(Da, &da));
355   PetscCall(VecRestoreArray(Db, &db));
356   PetscCall(VecRestoreArrayRead(X, &x));
357   PetscCall(VecRestoreArrayRead(Con, &f));
358   PetscCall(VecRestoreArrayRead(XL, &l));
359   PetscCall(VecRestoreArrayRead(XU, &u));
360   PetscCall(VecRestoreArrayRead(T2, &t2));
361   PetscFunctionReturn(PETSC_SUCCESS);
362 }
363 
364 /*@
365   MatDSFischer - Calculates an element of the B-subdifferential of the
366   smoothed Fischer-Burmeister function for complementarity problems.
367 
368   Collective
369 
370   Input Parameters:
371 + jac - the jacobian of f at X
372 . X   - current point
373 . Con - constraint function evaluated at X
374 . XL  - lower bounds
375 . XU  - upper bounds
376 . mu  - smoothing parameter
377 . T1  - work vector
378 - T2  - work vector
379 
380   Output Parameters:
381 + Da - diagonal perturbation component of the result
382 . Db - row scaling component of the result
383 - Dm - derivative with respect to scaling parameter
384 
385   Level: developer
386 
387 .seealso: `Mat`, `VecFischer()`, `VecDFischer()`, `MatDFischer()`
388 @*/
389 PetscErrorCode MatDSFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, PetscReal mu, Vec T1, Vec T2, Vec Da, Vec Db, Vec Dm)
390 {
391   PetscInt           i, nn;
392   const PetscScalar *x, *f, *l, *u;
393   PetscScalar       *da, *db, *dm;
394   PetscReal          ai, bi, ci, di, ei, fi;
395 
396   PetscFunctionBegin;
397   if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) {
398     PetscCall(VecZeroEntries(Dm));
399     PetscCall(MatDFischer(jac, X, Con, XL, XU, T1, T2, Da, Db));
400   } else {
401     PetscCall(VecGetLocalSize(X, &nn));
402     PetscCall(VecGetArrayRead(X, &x));
403     PetscCall(VecGetArrayRead(Con, &f));
404     PetscCall(VecGetArrayRead(XL, &l));
405     PetscCall(VecGetArrayRead(XU, &u));
406     PetscCall(VecGetArray(Da, &da));
407     PetscCall(VecGetArray(Db, &db));
408     PetscCall(VecGetArray(Dm, &dm));
409 
410     for (i = 0; i < nn; ++i) {
411       if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
412         da[i] = -mu;
413         db[i] = -1.0;
414         dm[i] = -x[i];
415       } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
416         bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
417         ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
418         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
419 
420         da[i] = bi / ai - 1.0;
421         db[i] = -PetscRealPart(f[i]) / ai - 1.0;
422         dm[i] = 2.0 * mu / ai;
423       } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
424         bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
425         ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
426         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
427 
428         da[i] = bi / ai - 1.0;
429         db[i] = PetscRealPart(f[i]) / ai - 1.0;
430         dm[i] = 2.0 * mu / ai;
431       } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
432         da[i] = -1.0;
433         db[i] = 0.0;
434         dm[i] = 0.0;
435       } else {
436         bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
437         ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
438         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
439 
440         ci = bi / ai + 1.0;
441         di = PetscRealPart(f[i]) / ai + 1.0;
442         fi = 2.0 * mu / ai;
443 
444         ei = SFischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]), mu);
445         ai = fischsnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei, mu);
446         ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
447 
448         bi = ei / ai - 1.0;
449         ei = 2.0 * mu / ei;
450         ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
451 
452         da[i] = ai + bi * ci;
453         db[i] = bi * di;
454         dm[i] = ei + bi * fi;
455       }
456     }
457 
458     PetscCall(VecRestoreArrayRead(X, &x));
459     PetscCall(VecRestoreArrayRead(Con, &f));
460     PetscCall(VecRestoreArrayRead(XL, &l));
461     PetscCall(VecRestoreArrayRead(XU, &u));
462     PetscCall(VecRestoreArray(Da, &da));
463     PetscCall(VecRestoreArray(Db, &db));
464     PetscCall(VecRestoreArray(Dm, &dm));
465   }
466   PetscFunctionReturn(PETSC_SUCCESS);
467 }
468 
469 static inline PetscReal ST_InternalPN(PetscScalar in, PetscReal lb, PetscReal ub)
470 {
471   return PetscMax(0, PetscRealPart(in) - ub) - PetscMax(0, -PetscRealPart(in) - PetscAbsReal(lb));
472 }
473 
474 static inline PetscReal ST_InternalNN(PetscScalar in, PetscReal lb, PetscReal ub)
475 {
476   return PetscMax(0, PetscRealPart(in) + PetscAbsReal(ub)) - PetscMax(0, -PetscRealPart(in) - PetscAbsReal(lb));
477 }
478 
479 static inline PetscReal ST_InternalPP(PetscScalar in, PetscReal lb, PetscReal ub)
480 {
481   return PetscMax(0, PetscRealPart(in) - ub) + PetscMin(0, PetscRealPart(in) - lb);
482 }
483 
484 /*@
485   TaoSoftThreshold - Calculates soft thresholding routine with input vector
486   and given lower and upper bound and returns it to output vector.
487 
488   Input Parameters:
489 + in - input vector to be thresholded
490 . lb - lower bound
491 - ub - upper bound
492 
493   Output Parameter:
494 . out - Soft thresholded output vector
495 
496   Notes:
497   Soft thresholding is defined as
498   \[ S(input,lb,ub) =
499   \begin{cases}
500   input - ub  \text{input > ub} \\
501   0           \text{lb =< input <= ub} \\
502   input + lb  \text{input < lb} \\
503   \]
504 
505   Level: developer
506 
507 .seealso: `Tao`, `Vec`
508 @*/
509 PetscErrorCode TaoSoftThreshold(Vec in, PetscReal lb, PetscReal ub, Vec out)
510 {
511   PetscInt     i, nlocal, mlocal;
512   PetscScalar *inarray, *outarray;
513 
514   PetscFunctionBegin;
515   PetscCall(VecGetArrayPair(in, out, &inarray, &outarray));
516   PetscCall(VecGetLocalSize(in, &nlocal));
517   PetscCall(VecGetLocalSize(out, &mlocal));
518 
519   PetscCheck(nlocal == mlocal, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Input and output vectors need to be of same size.");
520   if (lb == ub) {
521     PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray));
522     PetscFunctionReturn(PETSC_SUCCESS);
523   }
524   PetscCheck(lb < ub, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Lower bound needs to be lower than upper bound.");
525 
526   if (ub >= 0 && lb < 0) {
527     for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalPN(inarray[i], lb, ub);
528   } else if (ub < 0 && lb < 0) {
529     for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalNN(inarray[i], lb, ub);
530   } else {
531     for (i = 0; i < nlocal; i++) outarray[i] = ST_InternalPP(inarray[i], lb, ub);
532   }
533 
534   PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray));
535   PetscFunctionReturn(PETSC_SUCCESS);
536 }
537