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
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.
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]
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;
51: PetscFunctionBegin;
58: if (!L && !U) {
59: PetscCall(VecAXPBY(FB, -1.0, 0.0, F));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
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));
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");
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));
77: PetscCall(VecGetLocalSize(X, &n));
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]);
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: }
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: }
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: }
114: /*@
115: VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for
116: complementarity problems.
118: Logically Collective
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
127: Output Parameter:
128: . FB - The Smoothed Fischer-Burmeister function vector
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.
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]
143: Level: developer
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;
154: PetscFunctionBegin;
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));
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");
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));
175: PetscCall(VecGetLocalSize(X, &n));
177: for (i = 0; i < n; ++i) {
178: xval = PetscRealPart(*x++);
179: fval = PetscRealPart(*f++);
180: lval = PetscRealPart(*l++);
181: uval = PetscRealPart(*u++);
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;
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: }
210: static inline PetscReal fischnorm(PetscReal a, PetscReal b)
211: {
212: return PetscSqrtReal(a * a + b * b);
213: }
215: static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
216: {
217: return PetscSqrtReal(a * a + b * b + 2.0 * c * c);
218: }
220: /*@
221: MatDFischer - Calculates an element of the B-subdifferential of the
222: Fischer-Burmeister function for complementarity problems.
224: Collective
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
235: Output Parameters:
236: + Da - diagonal perturbation component of the result
237: - Db - row scaling component of the result
239: Level: developer
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;
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));
261: for (i = 0; i < nn; i++) {
262: da[i] = 0.0;
263: db[i] = 0.0;
264: t1[i] = 0.0;
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: }
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: }
279: PetscCall(VecRestoreArray(T1, &t1));
280: PetscCall(VecRestoreArrayRead(T2, &t2));
281: PetscCall(MatMult(jac, T1, T2));
282: PetscCall(VecGetArrayRead(T2, &t2));
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]));
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);
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]));
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);
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]));
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);
330: ci = bi / ai + 1.0;
331: di = PetscRealPart(f[i]) / ai + 1.0;
332: }
334: if (PetscRealPart(da[i]) >= 1) {
335: bi = ci + di * PetscRealPart(t2[i]);
336: ai = fischnorm(1.0, bi);
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);
345: bi = ei / ai - 1.0;
346: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
347: }
349: da[i] = ai + bi * ci;
350: db[i] = bi * di;
351: }
352: }
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: }
364: /*@
365: MatDSFischer - Calculates an element of the B-subdifferential of the
366: smoothed Fischer-Burmeister function for complementarity problems.
368: Collective
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
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
385: Level: developer
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;
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));
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);
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);
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);
440: ci = bi / ai + 1.0;
441: di = PetscRealPart(f[i]) / ai + 1.0;
442: fi = 2.0 * mu / ai;
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);
448: bi = ei / ai - 1.0;
449: ei = 2.0 * mu / ei;
450: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
452: da[i] = ai + bi * ci;
453: db[i] = bi * di;
454: dm[i] = ei + bi * fi;
455: }
456: }
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: }
469: static inline PetscReal ST_InternalPN(PetscScalar in, PetscReal lb, PetscReal ub)
470: {
471: return PetscMax(0, (PetscReal)PetscRealPart(in) - ub) - PetscMax(0, -(PetscReal)PetscRealPart(in) - PetscAbsReal(lb));
472: }
474: static inline PetscReal ST_InternalNN(PetscScalar in, PetscReal lb, PetscReal ub)
475: {
476: return PetscMax(0, (PetscReal)PetscRealPart(in) + PetscAbsReal(ub)) - PetscMax(0, -(PetscReal)PetscRealPart(in) - PetscAbsReal(lb));
477: }
479: static inline PetscReal ST_InternalPP(PetscScalar in, PetscReal lb, PetscReal ub)
480: {
481: return PetscMax(0, (PetscReal)PetscRealPart(in) - ub) + PetscMin(0, (PetscReal)PetscRealPart(in) - lb);
482: }
484: /*@
485: TaoSoftThreshold - Calculates soft thresholding routine with input vector
486: and given lower and upper bound and returns it to output vector.
488: Input Parameters:
489: + in - input vector to be thresholded
490: . lb - lower bound
491: - ub - upper bound
493: Output Parameter:
494: . out - Soft thresholded output vector
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: \]
505: Level: developer
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;
514: PetscFunctionBegin;
515: PetscCall(VecGetArrayPair(in, out, &inarray, &outarray));
516: PetscCall(VecGetLocalSize(in, &nlocal));
517: PetscCall(VecGetLocalSize(out, &mlocal));
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.");
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: }
534: PetscCall(VecRestoreArrayPair(in, out, &inarray, &outarray));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }