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