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: }