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");

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;

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));
256:   PetscCall(VecGetArray(Da, &da));
257:   PetscCall(VecGetArray(Db, &db));
258:   PetscCall(VecGetArray(T1, &t1));

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));
281:   PetscCall(MatMult(jac, T1, 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));
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));
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:     }

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