Actual source code: tao_util.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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: }