Actual source code: tao_util.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/petscimpl.h>
2: #include <petsctao.h> /*I "petsctao.h" I*/
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: }
16: /*@
17: VecFischer - Evaluates the Fischer-Burmeister function for complementarity
18: problems.
20: Logically Collective on vectors
22: Input Parameters:
23: + X - current point
24: . F - function evaluated at x
25: . L - lower bounds
26: - U - upper bounds
28: Output Parameters:
29: . FB - The Fischer-Burmeister function vector
31: Notes:
32: The Fischer-Burmeister function is defined as
33: $ phi(a,b) := sqrt(a*a + b*b) - a - b
34: and is used reformulate a complementarity problem as a semismooth
35: system of equations.
37: The result of this function is done by cases:
38: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i]
39: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i])
40: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i])
41: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u]))
42: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
44: Level: developer
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: PetscErrorCode ierr;
53: PetscInt low[5], high[5], n, i;
62: VecGetOwnershipRange(X, low, high);
63: VecGetOwnershipRange(F, low + 1, high + 1);
64: VecGetOwnershipRange(L, low + 2, high + 2);
65: VecGetOwnershipRange(U, low + 3, high + 3);
66: VecGetOwnershipRange(FB, low + 4, high + 4);
68: for (i = 1; i < 4; ++i) {
69: if (low[0] != low[i] || high[0] != high[i]) SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
70: }
72: VecGetArrayRead(X, &x);
73: VecGetArrayRead(F, &f);
74: VecGetArrayRead(L, &l);
75: VecGetArrayRead(U, &u);
76: VecGetArray(FB, &fb);
78: VecGetLocalSize(X, &n);
80: for (i = 0; i < n; ++i) {
81: xval = PetscRealPart(x[i]); fval = PetscRealPart(f[i]);
82: lval = PetscRealPart(l[i]); uval = PetscRealPart(u[i]);
84: if ((lval <= -PETSC_INFINITY) && (uval >= PETSC_INFINITY)) {
85: fb[i] = -fval;
86: } else if (lval <= -PETSC_INFINITY) {
87: fb[i] = -Fischer(uval - xval, -fval);
88: } else if (uval >= PETSC_INFINITY) {
89: fb[i] = Fischer(xval - lval, fval);
90: } else if (lval == uval) {
91: fb[i] = lval - xval;
92: } else {
93: fval = Fischer(uval - xval, -fval);
94: fb[i] = Fischer(xval - lval, fval);
95: }
96: }
98: VecRestoreArrayRead(X, &x);
99: VecRestoreArrayRead(F, &f);
100: VecRestoreArrayRead(L, &l);
101: VecRestoreArrayRead(U, &u);
102: VecRestoreArray(FB, &fb);
103: return(0);
104: }
106: PETSC_STATIC_INLINE PetscReal SFischer(PetscReal a, PetscReal b, PetscReal c)
107: {
108: /* Method suggested by Bob Vanderbei */
109: if (a + b <= 0) {
110: return PetscSqrtReal(a*a + b*b + 2.0*c*c) - (a + b);
111: }
112: return 2.0*(c*c - a*b) / (PetscSqrtReal(a*a + b*b + 2.0*c*c) + (a + b));
113: }
117: /*@
118: VecSFischer - Evaluates the Smoothed Fischer-Burmeister function for
119: complementarity problems.
121: Logically Collective on vectors
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 Parameters:
131: . FB - The Smoothed Fischer-Burmeister function vector
133: Notes:
134: The Smoothed Fischer-Burmeister function is defined as
135: $ phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b
136: and is used reformulate a complementarity problem as a semismooth
137: system of equations.
139: The result of this function is done by cases:
140: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i]
141: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu)
142: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu)
143: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
144: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
146: Level: developer
148: .seealso VecFischer()
149: @*/
150: PetscErrorCode VecSFischer(Vec X, Vec F, Vec L, Vec U, PetscReal mu, Vec FB)
151: {
152: const PetscScalar *x, *f, *l, *u;
153: PetscScalar *fb;
154: PetscReal xval, fval, lval, uval;
155: PetscErrorCode ierr;
156: PetscInt low[5], high[5], n, i;
165: VecGetOwnershipRange(X, low, high);
166: VecGetOwnershipRange(F, low + 1, high + 1);
167: VecGetOwnershipRange(L, low + 2, high + 2);
168: VecGetOwnershipRange(U, low + 3, high + 3);
169: VecGetOwnershipRange(FB, low + 4, high + 4);
171: for (i = 1; i < 4; ++i) {
172: if (low[0] != low[i] || high[0] != high[i]) SETERRQ(PETSC_COMM_SELF,1,"Vectors must be identically loaded over processors");
173: }
175: VecGetArrayRead(X, &x);
176: VecGetArrayRead(F, &f);
177: VecGetArrayRead(L, &l);
178: VecGetArrayRead(U, &u);
179: VecGetArray(FB, &fb);
181: VecGetLocalSize(X, &n);
183: for (i = 0; i < n; ++i) {
184: xval = PetscRealPart(*x++); fval = PetscRealPart(*f++);
185: lval = PetscRealPart(*l++); 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; f -= n; l -=n; u -= n; fb -= n;
202: VecRestoreArrayRead(X, &x);
203: VecRestoreArrayRead(F, &f);
204: VecRestoreArrayRead(L, &l);
205: VecRestoreArrayRead(U, &u);
206: VecRestoreArray(FB, &fb);
207: return(0);
208: }
210: PETSC_STATIC_INLINE PetscReal fischnorm(PetscReal a, PetscReal b)
211: {
212: return PetscSqrtReal(a*a + b*b);
213: }
215: PETSC_STATIC_INLINE PetscReal fischsnorm(PetscReal a, PetscReal b, PetscReal c)
216: {
217: return PetscSqrtReal(a*a + b*b + 2.0*c*c);
218: }
222: /*@
223: MatDFischer - Calculates an element of the B-subdifferential of the
224: Fischer-Burmeister function for complementarity problems.
226: Collective on jac
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: VecFischer()
244: @*/
245: PetscErrorCode MatDFischer(Mat jac, Vec X, Vec Con, Vec XL, Vec XU, Vec T1, Vec T2, Vec Da, Vec Db)
246: {
247: PetscErrorCode ierr;
248: PetscInt i,nn;
249: const PetscScalar *x,*f,*l,*u,*t2;
250: PetscScalar *da,*db,*t1;
251: PetscReal ai,bi,ci,di,ei;
254: VecGetLocalSize(X,&nn);
255: VecGetArrayRead(X,&x);
256: VecGetArrayRead(Con,&f);
257: VecGetArrayRead(XL,&l);
258: VecGetArrayRead(XU,&u);
259: VecGetArray(Da,&da);
260: VecGetArray(Db,&db);
261: VecGetArray(T1,&t1);
262: VecGetArrayRead(T2,&t2);
264: for (i = 0; i < nn; i++) {
265: da[i] = 0.0;
266: db[i] = 0.0;
267: t1[i] = 0.0;
269: if (PetscAbsScalar(f[i]) <= PETSC_MACHINE_EPSILON) {
270: if (PetscRealPart(l[i]) > PETSC_NINFINITY && PetscAbsScalar(x[i] - l[i]) <= PETSC_MACHINE_EPSILON) {
271: t1[i] = 1.0;
272: da[i] = 1.0;
273: }
275: if (PetscRealPart(u[i]) < PETSC_INFINITY && PetscAbsScalar(u[i] - x[i]) <= PETSC_MACHINE_EPSILON) {
276: t1[i] = 1.0;
277: db[i] = 1.0;
278: }
279: }
280: }
282: VecRestoreArray(T1,&t1);
283: VecRestoreArrayRead(T2,&t2);
284: MatMult(jac,T1,T2);
285: VecGetArrayRead(T2,&t2);
287: for (i = 0; i < nn; i++) {
288: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
289: da[i] = 0.0;
290: db[i] = -1.0;
291: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
292: if (PetscRealPart(db[i]) >= 1) {
293: ai = fischnorm(1.0, PetscRealPart(t2[i]));
295: da[i] = -1.0 / ai - 1.0;
296: db[i] = -t2[i] / ai - 1.0;
297: } else {
298: bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
299: ai = fischnorm(bi, PetscRealPart(f[i]));
300: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
302: da[i] = bi / ai - 1.0;
303: db[i] = -f[i] / ai - 1.0;
304: }
305: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) {
306: if (PetscRealPart(da[i]) >= 1) {
307: ai = fischnorm(1.0, PetscRealPart(t2[i]));
309: da[i] = 1.0 / ai - 1.0;
310: db[i] = t2[i] / ai - 1.0;
311: } else {
312: bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
313: ai = fischnorm(bi, PetscRealPart(f[i]));
314: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
316: da[i] = bi / ai - 1.0;
317: db[i] = f[i] / ai - 1.0;
318: }
319: } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
320: da[i] = -1.0;
321: db[i] = 0.0;
322: } else {
323: if (PetscRealPart(db[i]) >= 1) {
324: ai = fischnorm(1.0, PetscRealPart(t2[i]));
326: ci = 1.0 / ai + 1.0;
327: di = PetscRealPart(t2[i]) / ai + 1.0;
328: } else {
329: bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
330: ai = fischnorm(bi, PetscRealPart(f[i]));
331: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
333: ci = bi / ai + 1.0;
334: di = PetscRealPart(f[i]) / ai + 1.0;
335: }
337: if (PetscRealPart(da[i]) >= 1) {
338: bi = ci + di*PetscRealPart(t2[i]);
339: ai = fischnorm(1.0, bi);
341: bi = bi / ai - 1.0;
342: ai = 1.0 / ai - 1.0;
343: } else {
344: ei = Fischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]));
345: ai = fischnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei);
346: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
348: bi = ei / ai - 1.0;
349: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
350: }
352: da[i] = ai + bi*ci;
353: db[i] = bi*di;
354: }
355: }
357: VecRestoreArray(Da,&da);
358: VecRestoreArray(Db,&db);
359: VecRestoreArrayRead(X,&x);
360: VecRestoreArrayRead(Con,&f);
361: VecRestoreArrayRead(XL,&l);
362: VecRestoreArrayRead(XU,&u);
363: VecRestoreArrayRead(T2,&t2);
364: return(0);
365: }
369: /*@
370: MatDSFischer - Calculates an element of the B-subdifferential of the
371: smoothed Fischer-Burmeister function for complementarity problems.
373: Collective on jac
375: Input Parameters:
376: + jac - the jacobian of f at X
377: . X - current point
378: . F - constraint function evaluated at X
379: . XL - lower bounds
380: . XU - upper bounds
381: . mu - smoothing parameter
382: . T1 - work vector
383: - T2 - work vector
385: Output Parameter:
386: + Da - diagonal perturbation component of the result
387: . Db - row scaling component of the result
388: - Dm - derivative with respect to scaling parameter
390: Level: developer
392: .seealso MatDFischer()
393: @*/
394: PetscErrorCode MatDSFischer(Mat jac, Vec X, Vec Con,Vec XL, Vec XU, PetscReal mu,Vec T1, Vec T2,Vec Da, Vec Db, Vec Dm)
395: {
396: PetscErrorCode ierr;
397: PetscInt i,nn;
398: const PetscScalar *x, *f, *l, *u;
399: PetscScalar *da, *db, *dm;
400: PetscReal ai, bi, ci, di, ei, fi;
403: if (PetscAbsReal(mu) <= PETSC_MACHINE_EPSILON) {
404: VecZeroEntries(Dm);
405: MatDFischer(jac, X, Con, XL, XU, T1, T2, Da, Db);
406: } else {
407: VecGetLocalSize(X,&nn);
408: VecGetArrayRead(X,&x);
409: VecGetArrayRead(Con,&f);
410: VecGetArrayRead(XL,&l);
411: VecGetArrayRead(XU,&u);
412: VecGetArray(Da,&da);
413: VecGetArray(Db,&db);
414: VecGetArray(Dm,&dm);
416: for (i = 0; i < nn; ++i) {
417: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) {
418: da[i] = -mu;
419: db[i] = -1.0;
420: dm[i] = -x[i];
421: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) {
422: bi = PetscRealPart(u[i]) - PetscRealPart(x[i]);
423: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
424: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
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(u[i]) >= PETSC_INFINITY) {
430: bi = PetscRealPart(x[i]) - PetscRealPart(l[i]);
431: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
432: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
434: da[i] = bi / ai - 1.0;
435: db[i] = PetscRealPart(f[i]) / ai - 1.0;
436: dm[i] = 2.0 * mu / ai;
437: } else if (PetscRealPart(l[i]) == PetscRealPart(u[i])) {
438: da[i] = -1.0;
439: db[i] = 0.0;
440: dm[i] = 0.0;
441: } else {
442: bi = PetscRealPart(x[i]) - PetscRealPart(u[i]);
443: ai = fischsnorm(bi, PetscRealPart(f[i]), mu);
444: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
446: ci = bi / ai + 1.0;
447: di = PetscRealPart(f[i]) / ai + 1.0;
448: fi = 2.0 * mu / ai;
450: ei = SFischer(PetscRealPart(u[i]) - PetscRealPart(x[i]), -PetscRealPart(f[i]), mu);
451: ai = fischsnorm(PetscRealPart(x[i]) - PetscRealPart(l[i]), ei, mu);
452: ai = PetscMax(PETSC_MACHINE_EPSILON, ai);
454: bi = ei / ai - 1.0;
455: ei = 2.0 * mu / ei;
456: ai = (PetscRealPart(x[i]) - PetscRealPart(l[i])) / ai - 1.0;
458: da[i] = ai + bi*ci;
459: db[i] = bi*di;
460: dm[i] = ei + bi*fi;
461: }
462: }
464: VecRestoreArrayRead(X,&x);
465: VecRestoreArrayRead(Con,&f);
466: VecRestoreArrayRead(XL,&l);
467: VecRestoreArrayRead(XU,&u);
468: VecRestoreArray(Da,&da);
469: VecRestoreArray(Db,&db);
470: VecRestoreArray(Dm,&dm);
471: }
472: return(0);
473: }