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