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 $\phi(a,b) = \sqrt{a^2 + b^2 + 2 \mu^2} - a - b$
136: and is used reformulate a complementarity problem as a semismooth system of equations.
138: The result of this function is done by cases\:
139: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i]
140: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu)
141: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu)
142: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
143: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
145: Level: developer
147: .seealso: `Vec`, `VecFischer()`, `MatDFischer()`, `MatDSFischer()`
148: @*/
149: PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB)
150: {
151: const PetscScalar *x, *f, *l, *u;
152: PetscScalar *fb;
153: PetscReal xval, fval, lval, uval;
154: PetscInt low[5], high[5], n, i;
156: PetscFunctionBegin;
163: PetscCall(VecGetOwnershipRange(X, low, high));
164: PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
165: PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
166: PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
167: PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));
169: 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");
171: PetscCall(VecGetArrayRead(X, &x));
172: PetscCall(VecGetArrayRead(F, &f));
173: PetscCall(VecGetArrayRead(L, &l));
174: PetscCall(VecGetArrayRead(U, &u));
175: PetscCall(VecGetArray(FB, &fb));
177: PetscCall(VecGetLocalSize(X, &n));
179: for (i = 0; i < n; ++i) {
180: xval = PetscRealPart(*x++);
181: fval = PetscRealPart(*f++);
182: lval = PetscRealPart(*l++);
183: uval = PetscRealPart(*u++);
185: if (lval <= -PETSC_INFINITY && uval >= PETSC_INFINITY) {
186: (*fb++) = -fval - mu * xval;
187: } else if (lval <= -PETSC_INFINITY) {
188: (*fb++) = -SFischer(uval - xval, -fval, mu);
189: } else if (uval >= PETSC_INFINITY) {
190: (*fb++) = SFischer(xval - lval, fval, mu);
191: } else if (lval == uval) {
192: (*fb++) = lval - xval;
193: } else {
194: fval = SFischer(uval - xval, -fval, mu);
195: (*fb++) = SFischer(xval - lval, fval, mu);
196: }
197: }
198: x -= n;
199: f -= n;
200: l -= n;
201: u -= n;
202: fb -= n;
204: PetscCall(VecRestoreArrayRead(X, &x));
205: PetscCall(VecRestoreArrayRead(F, &f));
206: PetscCall(VecRestoreArrayRead(L, &l));
207: PetscCall(VecRestoreArrayRead(U, &u));
208: PetscCall(VecRestoreArray(FB, &fb));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static inline PetscReal fischnorm(PetscReal a, PetscReal b)
213: {
214: return PetscSqrtReal(a * a + b * b);
215: }
217: static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
218: {
219: return PetscSqrtReal(a * a + b * b + 2.0 * c * c);
220: }
222: /*@
223: MatDFischer - Calculates an element of the B-subdifferential of the
224: Fischer-Burmeister function for complementarity problems.
226: Collective
228: Input Parameters:
229: + jac - the jacobian of `f` at `X`
230: . X - current point
231: . Con - constraints function evaluated at `X`
232: . XL - lower bounds
233: . XU - upper bounds
234: . T1 - work vector
235: - T2 - work vector
237: Output Parameters:
238: + Da - diagonal perturbation component of the result
239: - Db - row scaling component of the result
241: Level: developer
243: .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDSFischer()`
244: @*/
245: PetscErrorCode MatDFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db)
246: {
247: PetscInt nn;
248: const PetscScalar *x, *f, *l, *u, *t2;
249: PetscScalar *da, *db, *t1;
250: PetscReal ai, bi, ci, di, ei;
252: PetscFunctionBegin;
253: PetscCall(VecGetLocalSize(X, &nn));
254: PetscCall(VecGetArrayRead(X, &x));
255: PetscCall(VecGetArrayRead(Con, &f));
256: PetscCall(VecGetArrayRead(XL, &l));
257: PetscCall(VecGetArrayRead(XU, &u));
258: PetscCall(VecGetArray(Da, &da));
259: PetscCall(VecGetArray(Db, &db));
260: PetscCall(VecGetArray(T1, &t1));
261: PetscCall(VecGetArrayRead(T2, &t2));
263: for (PetscInt i = 0; i < nn; i++) {
264: da[i] = 0.0;
265: db[i] = 0.0;
266: t1[i] = 0.0;
268: if (PetscAbsScalar(f[i]) <= PETSC_MACHINE_EPSILON) {
269: if (PetscRealPart(l[i]) > PETSC_NINFINITY && PetscAbsScalar(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) {
270: t1[i] = 1.0;
271: da[i] = 1.0;
272: }
274: if (PetscRealPart(u[i]) < PETSC_INFINITY && PetscAbsScalar(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) {
275: t1[i] = 1.0;
276: db[i] = 1.0;
277: }
278: }
279: }
281: PetscCall(VecRestoreArray(T1, &t1));
282: PetscCall(VecRestoreArrayRead(T2, &t2));
283: PetscCall(MatMult(jac, T1, T2));
284: PetscCall(VecGetArrayRead(T2, &t2));
286: for (PetscInt i = 0; i < nn; i++) {
287: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
288: da[i] = 0.0;
289: db[i] = -1.0;
290: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
291: if (PetscRealPart(db[i]) >= 1) {
292: ai = fischnorm(1.0, PetscRealPart(t2[i]));
294: da[i] = -1.0 / ai - 1.0;
295: db[i] = -t2[i] / ai - 1.0;
296: } else {
297: bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
298: ai = fischnorm(bi, PetscRealPart(f[i]));
299: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
301: da[i] = bi / ai - 1.0;
302: db[i] = -f[i] / ai - 1.0;
303: }
304: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
305: if (PetscRealPart(da[i]) >= 1) {
306: ai = fischnorm(1.0, PetscRealPart(t2[i]));
308: da[i] = 1.0 / ai - 1.0;
309: db[i] = t2[i] / ai - 1.0;
310: } else {
311: bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
312: ai = fischnorm(bi, PetscRealPart(f[i]));
313: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
315: da[i] = bi / ai - 1.0;
316: db[i] = f[i] / ai - 1.0;
317: }
318: } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
319: da[i] = -1.0;
320: db[i] = 0.0;
321: } else {
322: if (PetscRealPart(db[i]) >= 1) {
323: ai = fischnorm(1.0, PetscRealPart(t2[i]));
325: ci = 1.0 / ai + 1.0;
326: di = PetscRealPart(t2[i]) / ai + 1.0;
327: } else {
328: bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
329: ai = fischnorm(bi, PetscRealPart(f[i]));
330: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
332: ci = bi / ai + 1.0;
333: di = PetscRealPart(f[i]) / ai + 1.0;
334: }
336: if (PetscRealPart(da[i]) >= 1) {
337: bi = ci + di * PetscRealPart(t2[i]);
338: ai = fischnorm(1.0, bi);
340: bi = bi / ai - 1.0;
341: ai = 1.0 / ai - 1.0;
342: } else {
343: ei = Fischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]));
344: ai = fischnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei);
345: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
347: bi = ei / ai - 1.0;
348: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
349: }
351: da[i] = ai + bi * ci;
352: db[i] = bi * di;
353: }
354: }
356: PetscCall(VecRestoreArray(Da, &da));
357: PetscCall(VecRestoreArray(Db, &db));
358: PetscCall(VecRestoreArrayRead(X, &x));
359: PetscCall(VecRestoreArrayRead(Con, &f));
360: PetscCall(VecRestoreArrayRead(XL, &l));
361: PetscCall(VecRestoreArrayRead(XU, &u));
362: PetscCall(VecRestoreArrayRead(T2, &t2));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*@
367: MatDSFischer - Calculates an element of the B-subdifferential of the
368: smoothed Fischer-Burmeister function for complementarity problems.
370: Collective
372: Input Parameters:
373: + jac - the Jacobian of the function at `X`
374: . X - current point
375: . Con - constraint function evaluated at `X`
376: . XL - lower bounds
377: . XU - upper bounds
378: . mu - smoothing parameter
379: . T1 - work vector
380: - T2 - work vector
382: Output Parameters:
383: + Da - diagonal perturbation component of the result
384: . Db - row scaling component of the result
385: - Dm - derivative with respect to scaling parameter
387: Level: developer
389: .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDFischer()`
390: @*/
391: PetscErrorCode MatDSFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, PetscReal mu, Vec T1, Vec T2, Vec Da, Vec Db, Vec Dm)
392: {
393: PetscInt nn;
394: const PetscScalar *x, *f, *l, *u;
395: PetscScalar *da, *db, *dm;
396: PetscReal ai, bi, ci, di, ei, fi;
398: PetscFunctionBegin;
399: if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) {
400: PetscCall(VecZeroEntries(Dm));
401: PetscCall(MatDFischer(jac, X, Con, XL, XU, T1, T2, Da, Db));
402: } else {
403: PetscCall(VecGetLocalSize(X, &nn));
404: PetscCall(VecGetArrayRead(X, &x));
405: PetscCall(VecGetArrayRead(Con, &f));
406: PetscCall(VecGetArrayRead(XL, &l));
407: PetscCall(VecGetArrayRead(XU, &u));
408: PetscCall(VecGetArray(Da, &da));
409: PetscCall(VecGetArray(Db, &db));
410: PetscCall(VecGetArray(Dm, &dm));
412: for (PetscInt i = 0; i < nn; ++i) {
413: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
414: da[i] = -mu;
415: db[i] = -1.0;
416: dm[i] = -x[i];
417: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
418: bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
419: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
420: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
422: da[i] = bi / ai - 1.0;
423: db[i] = -PetscRealPart(f[i]) / ai - 1.0;
424: dm[i] = 2.0 * mu / ai;
425: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
426: bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
427: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
428: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
430: da[i] = bi / ai - 1.0;
431: db[i] = PetscRealPart(f[i]) / ai - 1.0;
432: dm[i] = 2.0 * mu / ai;
433: } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
434: da[i] = -1.0;
435: db[i] = 0.0;
436: dm[i] = 0.0;
437: } else {
438: bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
439: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
440: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
442: ci = bi / ai + 1.0;
443: di = PetscRealPart(f[i]) / ai + 1.0;
444: fi = 2.0 * mu / ai;
446: ei = SFischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]), mu);
447: ai = fischsnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei, mu);
448: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
450: bi = ei / ai - 1.0;
451: ei = 2.0 * mu / ei;
452: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
454: da[i] = ai + bi * ci;
455: db[i] = bi * di;
456: dm[i] = ei + bi * fi;
457: }
458: }
460: PetscCall(VecRestoreArrayRead(X, &x));
461: PetscCall(VecRestoreArrayRead(Con, &f));
462: PetscCall(VecRestoreArrayRead(XL, &l));
463: PetscCall(VecRestoreArrayRead(XU, &u));
464: PetscCall(VecRestoreArray(Da, &da));
465: PetscCall(VecRestoreArray(Db, &db));
466: PetscCall(VecRestoreArray(Dm, &dm));
467: }
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }