Actual source code: tao_util.c
1: #include <petsc/private/petscimpl.h>
2: #include <petsctao.h>
3: #include <petscsys.h>
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: }
12: /*@
13: VecFischer - Evaluates the Fischer-Burmeister function for complementarity
14: problems.
16: Logically Collective
18: Input Parameters:
19: + X - current point
20: . F - function evaluated at x
21: . L - lower bounds
22: - U - upper bounds
24: Output Parameter:
25: . FB - The Fischer-Burmeister function vector
27: Level: developer
29: Notes:
30: The Fischer-Burmeister function is defined as
32: $$
33: \phi(a,b) := \sqrt{a^2 + b^2} - a - b
34: $$
36: and is used reformulate a complementarity problem as a semismooth
37: system of equations.
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]
46: .seealso: `Vec`, `VecSFischer()`, `MatDFischer()`, `MatDSFischer()`
47: @*/
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;
55: PetscFunctionBegin;
62: if (!L && !U) {
63: PetscCall(VecAXPBY(FB, -1.0, 0.0, F));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
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));
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");
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));
81: PetscCall(VecGetLocalSize(X, &n));
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]);
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: }
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: }
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: }
118: /*@
119: VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for
120: complementarity problems.
122: Logically Collective
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
131: Output Parameter:
132: . FB - The Smoothed Fischer-Burmeister function vector
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.
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]
147: Level: developer
149: .seealso: `Vec`, `VecFischer()`, `MatDFischer()`, `MatDSFischer()`
150: @*/
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;
158: PetscFunctionBegin;
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));
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");
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));
179: PetscCall(VecGetLocalSize(X, &n));
181: for (i = 0; i < n; ++i) {
182: xval = PetscRealPart(*x++);
183: fval = PetscRealPart(*f++);
184: lval = PetscRealPart(*l++);
185: uval = PetscRealPart(*u++);
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;
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: }
214: static inline PetscReal fischnorm(PetscReal a, PetscReal b)
215: {
216: return PetscSqrtReal(a * a + b * b);
217: }
219: static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
220: {
221: return PetscSqrtReal(a * a + b * b + 2.0 * c * c);
222: }
224: /*@
225: MatDFischer - Calculates an element of the B-subdifferential of the
226: Fischer-Burmeister function for complementarity problems.
228: Collective
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
239: Output Parameters:
240: + Da - diagonal perturbation component of the result
241: - Db - row scaling component of the result
243: Level: developer
245: .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDSFischer()`
246: @*/
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;
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));
265: for (i = 0; i < nn; i++) {
266: da[i] = 0.0;
267: db[i] = 0.0;
268: t1[i] = 0.0;
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: }
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: }
283: PetscCall(VecRestoreArray(T1, &t1));
284: PetscCall(VecRestoreArrayRead(T2, &t2));
285: PetscCall(MatMult(jac, T1, T2));
286: PetscCall(VecGetArrayRead(T2, &t2));
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]));
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);
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]));
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);
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]));
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);
334: ci = bi / ai + 1.0;
335: di = PetscRealPart(f[i]) / ai + 1.0;
336: }
338: if (PetscRealPart(da[i]) >= 1) {
339: bi = ci + di * PetscRealPart(t2[i]);
340: ai = fischnorm(1.0, bi);
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);
349: bi = ei / ai - 1.0;
350: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
351: }
353: da[i] = ai + bi * ci;
354: db[i] = bi * di;
355: }
356: }
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: }
368: /*@
369: MatDSFischer - Calculates an element of the B-subdifferential of the
370: smoothed Fischer-Burmeister function for complementarity problems.
372: Collective
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
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
389: Level: developer
391: .seealso: `Mat`, `VecFischer()`, `VecDFischer()`, `MatDFischer()`
392: @*/
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;
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));
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);
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);
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);
444: ci = bi / ai + 1.0;
445: di = PetscRealPart(f[i]) / ai + 1.0;
446: fi = 2.0 * mu / ai;
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);
452: bi = ei / ai - 1.0;
453: ei = 2.0 * mu / ei;
454: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
456: da[i] = ai + bi * ci;
457: db[i] = bi * di;
458: dm[i] = ei + bi * fi;
459: }
460: }
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: }
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: }
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: }
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: }
488: /*@
489: TaoSoftThreshold - Calculates soft thresholding routine with input vector
490: and given lower and upper bound and returns it to output vector.
492: Input Parameters:
493: + in - input vector to be thresholded
494: . lb - lower bound
495: - ub - upper bound
497: Output Parameter:
498: . out - Soft thresholded output vector
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: \]
509: Level: developer
511: .seealso: `Tao`, `Vec`
512: @*/
513: PetscErrorCode TaoSoftThreshold(Vec in, PetscReal lb, PetscReal ub, Vec out)
514: {
515: PetscInt i, nlocal, mlocal;
516: PetscScalar *inarray, *outarray;
518: PetscFunctionBegin;
519: PetscCall(VecGetArrayPair(in, out, &inarray, &outarray));
520: PetscCall(VecGetLocalSize(in, &nlocal));
521: PetscCall(VecGetLocalSize(out, &mlocal));
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.");
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: }
538: PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }