Newer
Older
/* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */
/* Include "petsctao.h" so we can use TAO solvers. */
#include <petsctao.h>
static char help[] = "This example demonstrates use of the TAO package to \n\
solve an unconstrained minimization problem. This example is based on a \n\
problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
boundary values along the edges of the domain, the objective is to find the\n\
surface with the minimal area that satisfies the boundary conditions.\n\
The command line options are:\n\
-mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
-my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
-start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n";
/*
User-defined application context - contains data needed by the
application-provided call-back routines, FormFunctionGradient()
and FormHessian().
*/
typedef struct {
PetscInt mx, my; /* discretization in x, y directions */
PetscReal *bottom, *top, *left, *right; /* boundary values */
Mat H;
} AppCtx;
/* -------- User-defined Routines --------- */
static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
int main(int argc, char **argv)
{
PetscInt N; /* Size of vector */
PetscMPIInt size; /* Number of processors */
Vec x; /* solution, gradient vectors */
KSP ksp; /* PETSc Krylov subspace method */
PetscBool flg; /* A return value when checking for user options */
Tao tao; /* Tao solver context */
AppCtx user; /* user-defined work context */
PetscFunctionBeginUser;
PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");
/* Specify default dimension of the problem */
user.mx = 4;
user.my = 4;
/* Check for any command line arguments that override defaults */
PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n"));
PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my));
/* Calculate any derived values from parameters */
N = user.mx * user.my;
/* Create TAO solver and set desired solution method */
PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
PetscCall(TaoSetType(tao, TAOLMVM));
/* Initialize minsurf application data structure for use in the function evaluations */
PetscCall(MSA_BoundaryConditions(&user)); /* Application specific routine */
/*
Create a vector to hold the variables. Compute an initial solution.
Set this vector, which will be used by TAO.
*/
PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
PetscCall(MSA_InitialPoint(&user, x)); /* Application specific routine */
PetscCall(TaoSetSolution(tao, x)); /* A TAO routine */
/* Provide TAO routines for function, gradient, and Hessian evaluation */
PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
/* Create a matrix data structure to store the Hessian. This structure will be used by TAO */
PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 7, NULL, &user.H));
PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
/* Check for any TAO command line options */
/* Limit the number of iterations in the KSP linear solver */
PetscCall(TaoGetKSP(tao, &ksp));
if (ksp) PetscCall(KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, user.mx * user.my));
PetscCall(TaoDestroy(&tao));
PetscCall(VecDestroy(&x));
PetscCall(MatDestroy(&user.H));
PetscCall(PetscFree(user.bottom));
PetscCall(PetscFree(user.top));
PetscCall(PetscFree(user.left));
PetscCall(PetscFree(user.right));
}
/* -------------------------------------------------------------------- */
/* FormFunctionGradient - Evaluates function and corresponding gradient.
. X - input vector
. userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
Output Parameters:
. fcn - the newly evaluated function
. G - vector containing the newly evaluated gradient
*/
PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
{
AppCtx *user = (AppCtx *)userCtx;
PetscInt i, j, row;
PetscInt mx = user->mx, my = user->my;
PetscReal rhx = mx + 1, rhy = my + 1;
PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy, ft = 0;
PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
PetscReal zero = 0.0;
Stefano Zampini
committed
PetscScalar *g;
const PetscScalar *x;
PetscFunctionBeginUser;
PetscCall(VecGetArrayRead(X, &x));
PetscCall(VecGetArray(G, &g));
/* Compute function over the locally owned part of the mesh */
for (j = 0; j < my; j++) {
for (i = 0; i < mx; i++) {
row = (j)*mx + (i);
xc = x[row];
xlt = xrb = xl = xr = xb = xt = xc;
if (i == 0) { /* left side */
xl = user->left[j + 1];
xlt = user->left[j + 2];
xl = x[row - 1];
if (j == 0) { /* bottom side */
xb = user->bottom[i + 1];
xrb = user->bottom[i + 2];
xb = x[row - mx];
if (i + 1 == mx) { /* right side */
xr = user->right[j + 1];
xr = x[row + 1];
if (j + 1 == 0 + my) { /* top side */
xt = user->top[i + 1];
} else {
xt = x[row + mx];
if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];
d1 = (xc - xl);
d2 = (xc - xr);
d3 = (xc - xt);
d4 = (xc - xb);
d5 = (xr - xrb);
d6 = (xrb - xb);
d7 = (xlt - xl);
d8 = (xt - xlt);
df1dxc = d1 * hydhx;
df2dxc = (d1 * hydhx + d4 * hxdhy);
df3dxc = d3 * hxdhy;
df4dxc = (d2 * hydhx + d3 * hxdhy);
df5dxc = d2 * hydhx;
df6dxc = d4 * hxdhy;
d1 *= rhx;
d2 *= rhx;
d3 *= rhy;
d4 *= rhy;
d5 *= rhy;
d6 *= rhx;
d7 *= rhy;
d8 *= rhx;
f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
ft = ft + (f2 + f4);
df1dxc /= f1;
df2dxc /= f2;
df3dxc /= f3;
df4dxc /= f4;
df5dxc /= f5;
df6dxc /= f6;
g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
for (j = 0; j < my; j++) { /* left side */
d3 = (user->left[j + 1] - user->left[j + 2]) * rhy;
d2 = (user->left[j + 1] - x[j * mx]) * rhx;
ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
for (i = 0; i < mx; i++) { /* bottom */
d2 = (user->bottom[i + 1] - user->bottom[i + 2]) * rhx;
d3 = (user->bottom[i + 1] - x[i]) * rhy;
ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
for (j = 0; j < my; j++) { /* right side */
d1 = (x[(j + 1) * mx - 1] - user->right[j + 1]) * rhx;
d4 = (user->right[j] - user->right[j + 1]) * rhy;
ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
for (i = 0; i < mx; i++) { /* top side */
d1 = (x[(my - 1) * mx + i] - user->top[i + 1]) * rhy;
d4 = (user->top[i + 1] - user->top[i]) * rhx;
ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
d1 = (user->left[0] - user->left[1]) * rhy;
d2 = (user->bottom[0] - user->bottom[1]) * rhx;
ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
d1 = (user->right[my + 1] - user->right[my]) * rhy;
d2 = (user->top[mx + 1] - user->top[mx]) * rhx;
ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
(*fcn) = ft * area;
PetscCall(VecRestoreArrayRead(X, &x));
PetscCall(VecRestoreArray(G, &g));
PetscCall(PetscLogFlops(67.0 * mx * my));
PetscFunctionReturn(PETSC_SUCCESS);
}
/* ------------------------------------------------------------------- */
/*
FormHessian - Evaluates the Hessian matrix.
Input Parameters:
. x - input vector
. ptr - optional user-defined context, as set by TaoSetHessian()
Output Parameters:
. H - Hessian matrix
. Hpre - optionally different preconditioning matrix
. flg - flag indicating matrix structure
*/
PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
{
AppCtx *user = (AppCtx *)ptr;
Stefano Zampini
committed
PetscFunctionBeginUser;
PetscCall(QuadraticH(user, X, H));
PetscFunctionReturn(PETSC_SUCCESS);
}
/* ------------------------------------------------------------------- */
/*
QuadraticH - Evaluates the Hessian matrix.
Input Parameters:
. user - user-defined context, as set by TaoSetHessian()
. X - input vector
Output Parameter:
. H - Hessian matrix
*/
PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
{
PetscInt i, j, k, row;
PetscInt mx = user->mx, my = user->my;
PetscInt col[7];
PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
PetscReal rhx = mx + 1, rhy = my + 1;
PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
PetscReal hl, hr, ht, hb, hc, htl, hbr;
Stefano Zampini
committed
const PetscScalar *x;
PetscReal v[7];
PetscCall(VecGetArrayRead(X, &x));
/* Initialize matrix entries to zero */
PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
/* Compute Hessian over the locally owned part of the mesh */
for (i = 0; i < mx; i++) {
for (j = 0; j < my; j++) {
row = (j)*mx + (i);
xc = x[row];
xlt = xrb = xl = xr = xb = xt = xc;
if (i == 0) {
xl = user->left[j + 1];
xlt = user->left[j + 2];
xl = x[row - 1];
if (j == 0) {
xb = user->bottom[i + 1];
xrb = user->bottom[i + 2];
xb = x[row - mx];
if (i + 1 == mx) {
xr = user->right[j + 1];
xr = x[row + 1];
if (j + 1 == my) {
xt = user->top[i + 1];
} else {
xt = x[row + mx];
if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
d1 = (xc - xl) * rhx;
d2 = (xc - xr) * rhx;
d3 = (xc - xt) * rhy;
d4 = (xc - xb) * rhy;
d5 = (xrb - xr) * rhy;
d6 = (xrb - xb) * rhx;
d7 = (xlt - xl) * rhy;
d8 = (xlt - xt) * rhx;
f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2 * d2 * d3) / (f4 * f4 * f4);
hl *= 0.5;
hr *= 0.5;
ht *= 0.5;
hb *= 0.5;
hbr *= 0.5;
htl *= 0.5;
hc *= 0.5;
k = 0;
if (j > 0) {
v[k] = hb;
col[k] = row - mx;
k++;
if (j > 0 && i < mx - 1) {
v[k] = hbr;
col[k] = row - mx + 1;
k++;
if (i > 0) {
v[k] = hl;
col[k] = row - 1;
k++;
v[k] = hc;
col[k] = row;
k++;
if (i < mx - 1) {
v[k] = hr;
col[k] = row + 1;
k++;
if (i > 0 && j < my - 1) {
v[k] = htl;
col[k] = row + mx - 1;
k++;
if (j < my - 1) {
v[k] = ht;
col[k] = row + mx;
k++;
Set matrix values using local numbering, which was defined
earlier, in the main routine.
PetscCall(MatSetValues(Hessian, 1, &row, k, col, v, INSERT_VALUES));
PetscCall(VecRestoreArrayRead(X, &x));
PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
PetscCall(PetscLogFlops(199.0 * mx * my));
PetscFunctionReturn(PETSC_SUCCESS);
}
/* ------------------------------------------------------------------- */
MSA_BoundaryConditions - Calculates the boundary conditions for
the region.
Input Parameter:
. user - user-defined application context
Output Parameter:
. user - user-defined application context
*/
static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
{
PetscInt i, j, k, limit = 0;
PetscInt maxits = 5;
PetscInt mx = user->mx, my = user->my;
PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
PetscReal *boundary;
Stefano Zampini
committed
PetscFunctionBeginUser;
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
bsize = mx + 2;
lsize = my + 2;
rsize = my + 2;
tsize = mx + 2;
PetscCall(PetscMalloc1(bsize, &user->bottom));
PetscCall(PetscMalloc1(tsize, &user->top));
PetscCall(PetscMalloc1(lsize, &user->left));
PetscCall(PetscMalloc1(rsize, &user->right));
hx = (r - l) / (mx + 1);
hy = (t - b) / (my + 1);
for (j = 0; j < 4; j++) {
if (j == 0) {
yt = b;
xt = l;
limit = bsize;
boundary = user->bottom;
} else if (j == 1) {
yt = t;
xt = l;
limit = tsize;
boundary = user->top;
} else if (j == 2) {
yt = b;
xt = l;
limit = lsize;
boundary = user->left;
} else { /* if (j==3) */
yt = b;
xt = r;
limit = rsize;
boundary = user->right;
for (i = 0; i < limit; i++) {
u1 = xt;
u2 = -yt;
for (k = 0; k < maxits; k++) {
nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
njac11 = one + u2 * u2 - u1 * u1;
njac12 = two * u1 * u2;
njac21 = -two * u1 * u2;
njac22 = -one - u1 * u1 + u2 * u2;
det = njac11 * njac22 - njac21 * njac12;
u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
boundary[i] = u1 * u1 - u2 * u2;
if (j == 0 || j == 1) {
xt = xt + hx;
yt = yt + hy;
PetscFunctionReturn(PETSC_SUCCESS);
}
/* ------------------------------------------------------------------- */
/*
MSA_InitialPoint - Calculates the initial guess in one of three ways.
Input Parameters:
. user - user-defined application context
. X - vector for initial guess
Output Parameters:
. X - newly computed initial guess
*/
static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
{
PetscInt start = -1, i, j;
PetscReal zero = 0.0;
PetscBool flg;
Stefano Zampini
committed
PetscFunctionBeginUser;
PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));
if (flg && start == 0) { /* The zero vector is reasonable */
PetscCall(VecSet(X, zero));
} else { /* Take an average of the boundary conditions */
PetscInt row;
PetscInt mx = user->mx, my = user->my;
PetscCall(VecGetArray(X, &x));
/* Perform local computations */
for (j = 0; j < my; j++) {
for (i = 0; i < mx; i++) {
row = (j)*mx + (i);
x[row] = (((j + 1) * user->bottom[i + 1] + (my - j + 1) * user->top[i + 1]) / (my + 2) + ((i + 1) * user->left[j + 1] + (mx - i + 1) * user->right[j + 1]) / (mx + 2)) / 2.0;
PetscCall(VecRestoreArray(X, &x));
PetscFunctionReturn(PETSC_SUCCESS);
/*TEST
build:
requires: !complex
test:
args: -tao_monitor_short -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4
requires: !single
test:
suffix: 2
args: -tao_monitor_short -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
requires: !single
test:
suffix: 3
args: -tao_monitor_short -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
Barry Smith
committed
test:
suffix: 4
args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4
test:
suffix: 4_ew
args: -tao_monitor_short -tao_type bntr -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
requires: !single
test:
suffix: 5
args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4
test:
suffix: 5_ew
args: -tao_monitor_short -tao_type bntl -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
requires: !single
test:
suffix: 6
args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4
test:
suffix: 6_ew
args: -tao_monitor_short -tao_type bnls -tao_ksp_ew -tao_bnk_ksp_ew_version 3 -mx 10 -my 8 -tao_gatol 1.e-4
requires: !single
test:
suffix: 7
args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
test:
suffix: 8
args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
test:
suffix: 9
args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
args: -tao_monitor_short -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian
args: -tao_monitor_short -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
requires: !single
test:
suffix: 12
args: -tao_monitor_short -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian