Actual source code: tao_util.c
1: #include <petsc/private/petscimpl.h>
2: #include <petsctao.h>
4: static inline PetscReal Fischer(PetscReal a, PetscReal b)
5: {
6: /* Method suggested by Bob Vanderbei */
7: if (a + b <= 0) return PetscSqrtReal(a * a + b * b) - (a + b);
8: return -2.0 * a * b / (PetscSqrtReal(a * a + b * b) + (a + b));
9: }
11: /*@
12: VecFischer - Evaluates the Fischer-Burmeister function for complementarity
13: problems.
15: Logically Collective
17: Input Parameters:
18: + X - current point
19: . F - function evaluated at `X`
20: . L - lower bounds
21: - U - upper bounds
23: Output Parameter:
24: . FB - The Fischer-Burmeister function vector
26: Level: developer
28: Notes:
29: The Fischer-Burmeister function is defined as
31: $$
32: \phi(a,b) = \sqrt{a^2 + b^2} - a - b
33: $$
35: and is used reformulate a complementarity problem as a semismooth
36: system of equations.
38: The result of this function is done by cases\:
39: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i]
40: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i])
41: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i])
42: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u]))
43: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
45: .seealso: `Vec`, `VecSFischer()`, `MatDFischer()`, `MatDSFischer()`
46: @*/
47: PetscErrorCode VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FB)
48: {
49: const PetscScalar *x, *f, *l, *u;
50: PetscScalar *fb;
51: PetscReal xval, fval, lval, uval;
52: PetscInt low[5], high[5], n, i;
54: PetscFunctionBegin;
61: if (!L && !U) {
62: PetscCall(VecAXPBY(FB, -1.0, 0.0, F));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: PetscCall(VecGetOwnershipRange(X, low, high));
67: PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
68: PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
69: PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
70: PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));
72: 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: PetscCall(VecGetArrayRead(X, &x));
75: PetscCall(VecGetArrayRead(F, &f));
76: PetscCall(VecGetArrayRead(L, &l));
77: PetscCall(VecGetArrayRead(U, &u));
78: PetscCall(VecGetArray(FB, &fb));
80: PetscCall(VecGetLocalSize(X, &n));
82: for (i = 0; i < n; ++i) {
83: xval = PetscRealPart(x[i]);
84: fval = PetscRealPart(f[i]);
85: lval = PetscRealPart(l[i]);
86: uval = PetscRealPart(u[i]);
88: if (lval <= -PETSC_INFINITY && uval >= PETSC_INFINITY) {
89: fb[i] = -fval;
90: } else if (lval <= -PETSC_INFINITY) {
91: fb[i] = -Fischer(uval - xval, -fval);
92: } else if (uval >= PETSC_INFINITY) {
93: fb[i] = Fischer(xval - lval, fval);
94: } else if (lval == uval) {
95: fb[i] = lval - xval;
96: } else {
97: fval = Fischer(uval - xval, -fval);
98: fb[i] = Fischer(xval - lval, fval);
99: }
100: }
102: PetscCall(VecRestoreArrayRead(X, &x));
103: PetscCall(VecRestoreArrayRead(F, &f));
104: PetscCall(VecRestoreArrayRead(L, &l));
105: PetscCall(VecRestoreArrayRead(U, &u));
106: PetscCall(VecRestoreArray(FB, &fb));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static inline PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c)
111: {
112: /* Method suggested by Bob Vanderbei */
113: if (a + b <= 0) return PetscSqrtReal(a * a + b * b + 2.0 * c * c) - (a + b);
114: return 2.0 * (c * c - a * b) / (PetscSqrtReal(a * a + b * b + 2.0 * c * c) + (a + b));
115: }
117: /*@
118: VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for
119: complementarity problems.
121: Logically Collective
123: Input Parameters:
124: + X - current point
125: . F - function evaluated at `X`
126: . L - lower bounds
127: . U - upper bounds
128: - mu - smoothing parameter
130: Output Parameter:
131: . FB - The Smoothed Fischer-Burmeister function vector
133: Notes:
134: The Smoothed Fischer-Burmeister function is defined as $\phi(a,b) = \sqrt{a^2 + b^2 + 2 \mu^2} - a - b$
135: and is used reformulate a complementarity problem as a semismooth system of equations.
137: The result of this function is done by cases\:
138: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i]
139: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu)
140: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu)
141: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
142: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
144: Level: developer
146: .seealso: `Vec`, `VecFischer()`, `MatDFischer()`, `MatDSFischer()`
147: @*/
148: PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB)
149: {
150: const PetscScalar *x, *f, *l, *u;
151: PetscScalar *fb;
152: PetscReal xval, fval, lval, uval;
153: PetscInt low[5], high[5], n, i;
155: PetscFunctionBegin;
162: PetscCall(VecGetOwnershipRange(X, low, high));
163: PetscCall(VecGetOwnershipRange(F, low + 1, high + 1));
164: PetscCall(VecGetOwnershipRange(L, low + 2, high + 2));
165: PetscCall(VecGetOwnershipRange(U, low + 3, high + 3));
166: PetscCall(VecGetOwnershipRange(FB, low + 4, high + 4));
168: 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");
170: PetscCall(VecGetArrayRead(X, &x));
171: PetscCall(VecGetArrayRead(F, &f));
172: PetscCall(VecGetArrayRead(L, &l));
173: PetscCall(VecGetArrayRead(U, &u));
174: PetscCall(VecGetArray(FB, &fb));
176: PetscCall(VecGetLocalSize(X, &n));
178: for (i = 0; i < n; ++i) {
179: xval = PetscRealPart(*x++);
180: fval = PetscRealPart(*f++);
181: lval = PetscRealPart(*l++);
182: uval = PetscRealPart(*u++);
184: if (lval <= -PETSC_INFINITY && uval >= PETSC_INFINITY) {
185: (*fb++) = -fval - mu * xval;
186: } else if (lval <= -PETSC_INFINITY) {
187: (*fb++) = -SFischer(uval - xval, -fval, mu);
188: } else if (uval >= PETSC_INFINITY) {
189: (*fb++) = SFischer(xval - lval, fval, mu);
190: } else if (lval == uval) {
191: (*fb++) = lval - xval;
192: } else {
193: fval = SFischer(uval - xval, -fval, mu);
194: (*fb++) = SFischer(xval - lval, fval, mu);
195: }
196: }
197: x -= n;
198: f -= n;
199: l -= n;
200: u -= n;
201: fb -= n;
203: PetscCall(VecRestoreArrayRead(X, &x));
204: PetscCall(VecRestoreArrayRead(F, &f));
205: PetscCall(VecRestoreArrayRead(L, &l));
206: PetscCall(VecRestoreArrayRead(U, &u));
207: PetscCall(VecRestoreArray(FB, &fb));
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: static inline PetscReal fischnorm(PetscReal a, PetscReal b)
212: {
213: return PetscSqrtReal(a * a + b * b);
214: }
216: static inline PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
217: {
218: return PetscSqrtReal(a * a + b * b + 2.0 * c * c);
219: }
221: /*@
222: MatDFischer - Calculates an element of the B-subdifferential of the
223: Fischer-Burmeister function for complementarity problems.
225: Collective
227: Input Parameters:
228: + jac - the jacobian of `f` at `X`
229: . X - current point
230: . Con - constraints function evaluated at `X`
231: . XL - lower bounds
232: . XU - upper bounds
233: . T1 - work vector
234: - T2 - work vector
236: Output Parameters:
237: + Da - diagonal perturbation component of the result
238: - Db - row scaling component of the result
240: Level: developer
242: .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDSFischer()`
243: @*/
244: PetscErrorCode MatDFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db)
245: {
246: PetscInt nn;
247: const PetscScalar *x, *f, *l, *u, *t2;
248: PetscScalar *da, *db, *t1;
249: PetscReal ai, bi, ci, di, ei;
251: PetscFunctionBegin;
252: PetscCall(VecGetLocalSize(X, &nn));
253: PetscCall(VecGetArrayRead(X, &x));
254: PetscCall(VecGetArrayRead(Con, &f));
255: PetscCall(VecGetArrayRead(XL, &l));
256: PetscCall(VecGetArrayRead(XU, &u));
257: PetscCall(VecGetArray(Da, &da));
258: PetscCall(VecGetArray(Db, &db));
259: PetscCall(VecGetArray(T1, &t1));
260: PetscCall(VecGetArrayRead(T2, &t2));
262: for (PetscInt i = 0; i < nn; i++) {
263: da[i] = 0.0;
264: db[i] = 0.0;
265: t1[i] = 0.0;
267: if (PetscAbsScalar(f[i]) <= PETSC_MACHINE_EPSILON) {
268: if (PetscRealPart(l[i]) > PETSC_NINFINITY && PetscAbsScalar(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) {
269: t1[i] = 1.0;
270: da[i] = 1.0;
271: }
273: if (PetscRealPart(u[i]) < PETSC_INFINITY && PetscAbsScalar(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) {
274: t1[i] = 1.0;
275: db[i] = 1.0;
276: }
277: }
278: }
280: PetscCall(VecRestoreArray(T1, &t1));
281: PetscCall(VecRestoreArrayRead(T2, &t2));
282: PetscCall(MatMult(jac, T1, T2));
283: PetscCall(VecGetArrayRead(T2, &t2));
285: for (PetscInt i = 0; i < nn; i++) {
286: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
287: da[i] = 0.0;
288: db[i] = -1.0;
289: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
290: if (PetscRealPart(db[i]) >= 1) {
291: ai = fischnorm(1.0, PetscRealPart(t2[i]));
293: da[i] = -1.0 / ai - 1.0;
294: db[i] = -t2[i] / ai - 1.0;
295: } else {
296: bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
297: ai = fischnorm(bi, PetscRealPart(f[i]));
298: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
300: da[i] = bi / ai - 1.0;
301: db[i] = -f[i] / ai - 1.0;
302: }
303: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
304: if (PetscRealPart(da[i]) >= 1) {
305: ai = fischnorm(1.0, PetscRealPart(t2[i]));
307: da[i] = 1.0 / ai - 1.0;
308: db[i] = t2[i] / ai - 1.0;
309: } else {
310: bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
311: ai = fischnorm(bi, PetscRealPart(f[i]));
312: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
314: da[i] = bi / ai - 1.0;
315: db[i] = f[i] / ai - 1.0;
316: }
317: } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
318: da[i] = -1.0;
319: db[i] = 0.0;
320: } else {
321: if (PetscRealPart(db[i]) >= 1) {
322: ai = fischnorm(1.0, PetscRealPart(t2[i]));
324: ci = 1.0 / ai + 1.0;
325: di = PetscRealPart(t2[i]) / ai + 1.0;
326: } else {
327: bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
328: ai = fischnorm(bi, PetscRealPart(f[i]));
329: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
331: ci = bi / ai + 1.0;
332: di = PetscRealPart(f[i]) / ai + 1.0;
333: }
335: if (PetscRealPart(da[i]) >= 1) {
336: bi = ci + di * PetscRealPart(t2[i]);
337: ai = fischnorm(1.0, bi);
339: bi = bi / ai - 1.0;
340: ai = 1.0 / ai - 1.0;
341: } else {
342: ei = Fischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]));
343: ai = fischnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei);
344: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
346: bi = ei / ai - 1.0;
347: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
348: }
350: da[i] = ai + bi * ci;
351: db[i] = bi * di;
352: }
353: }
355: PetscCall(VecRestoreArray(Da, &da));
356: PetscCall(VecRestoreArray(Db, &db));
357: PetscCall(VecRestoreArrayRead(X, &x));
358: PetscCall(VecRestoreArrayRead(Con, &f));
359: PetscCall(VecRestoreArrayRead(XL, &l));
360: PetscCall(VecRestoreArrayRead(XU, &u));
361: PetscCall(VecRestoreArrayRead(T2, &t2));
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: MatDSFischer - Calculates an element of the B-subdifferential of the
367: smoothed Fischer-Burmeister function for complementarity problems.
369: Collective
371: Input Parameters:
372: + jac - the Jacobian of the function at `X`
373: . X - current point
374: . Con - constraint function evaluated at `X`
375: . XL - lower bounds
376: . XU - upper bounds
377: . mu - smoothing parameter
378: . T1 - work vector
379: - T2 - work vector
381: Output Parameters:
382: + Da - diagonal perturbation component of the result
383: . Db - row scaling component of the result
384: - Dm - derivative with respect to scaling parameter
386: Level: developer
388: .seealso: `Mat`, `VecFischer()`, `VecSFischer()`, `MatDFischer()`
389: @*/
390: PetscErrorCode MatDSFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, PetscReal mu, Vec T1, Vec T2, Vec Da, Vec Db, Vec Dm)
391: {
392: PetscInt nn;
393: const PetscScalar *x, *f, *l, *u;
394: PetscScalar *da, *db, *dm;
395: PetscReal ai, bi, ci, di, ei, fi;
397: PetscFunctionBegin;
398: if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) {
399: PetscCall(VecZeroEntries(Dm));
400: PetscCall(MatDFischer(jac, X, Con, XL, XU, T1, T2, Da, Db));
401: } else {
402: PetscCall(VecGetLocalSize(X, &nn));
403: PetscCall(VecGetArrayRead(X, &x));
404: PetscCall(VecGetArrayRead(Con, &f));
405: PetscCall(VecGetArrayRead(XL, &l));
406: PetscCall(VecGetArrayRead(XU, &u));
407: PetscCall(VecGetArray(Da, &da));
408: PetscCall(VecGetArray(Db, &db));
409: PetscCall(VecGetArray(Dm, &dm));
411: for (PetscInt i = 0; i < nn; ++i) {
412: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
413: da[i] = -mu;
414: db[i] = -1.0;
415: dm[i] = -x[i];
416: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
417: bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
418: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
419: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
421: da[i] = bi / ai - 1.0;
422: db[i] = -PetscRealPart(f[i]) / ai - 1.0;
423: dm[i] = 2.0 * mu / ai;
424: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
425: bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
426: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
427: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
429: da[i] = bi / ai - 1.0;
430: db[i] = PetscRealPart(f[i]) / ai - 1.0;
431: dm[i] = 2.0 * mu / ai;
432: } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
433: da[i] = -1.0;
434: db[i] = 0.0;
435: dm[i] = 0.0;
436: } else {
437: bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
438: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
439: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
441: ci = bi / ai + 1.0;
442: di = PetscRealPart(f[i]) / ai + 1.0;
443: fi = 2.0 * mu / ai;
445: ei = SFischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]), mu);
446: ai = fischsnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei, mu);
447: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
449: bi = ei / ai - 1.0;
450: ei = 2.0 * mu / ei;
451: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
453: da[i] = ai + bi * ci;
454: db[i] = bi * di;
455: dm[i] = ei + bi * fi;
456: }
457: }
459: PetscCall(VecRestoreArrayRead(X, &x));
460: PetscCall(VecRestoreArrayRead(Con, &f));
461: PetscCall(VecRestoreArrayRead(XL, &l));
462: PetscCall(VecRestoreArrayRead(XU, &u));
463: PetscCall(VecRestoreArray(Da, &da));
464: PetscCall(VecRestoreArray(Db, &db));
465: PetscCall(VecRestoreArray(Dm, &dm));
466: }
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }