Actual source code: rosenbrock2.c

  1: /* Program usage: mpiexec -n 1 rosenbrock2 [-help] [all TAO options] */

  3: /*  Include "petsctao.h" so we can use TAO solvers.  */
  4: #include <petsctao.h>

  6: static char help[] = "This example demonstrates use of the TAO package to \n\
  7: solve an unconstrained minimization problem on a single processor.  We \n\
  8: minimize the extended Rosenbrock function: \n\
  9:    sum_{i=0}^{n/2-1} (alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2) \n\
 10: or the chained Rosenbrock function:\n\
 11:    sum_{i=0}^{n-1} alpha*(x_{i+1} - x_i^2)^2 + (1 - x_i)^2\n";

 13: /*
 14:    User-defined application context - contains data needed by the
 15:    application-provided call-back routines that evaluate the function,
 16:    gradient, and hessian.
 17: */
 18: typedef struct {
 19:   PetscInt  n;     /* dimension */
 20:   PetscReal alpha; /* condition parameter */
 21:   PetscBool chained;
 22: } AppCtx;

 24: /* -------------- User-defined routines ---------- */
 25: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 26: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);

 28: int main(int argc, char **argv)
 29: {
 30:   PetscReal          zero = 0.0;
 31:   Vec                x; /* solution vector */
 32:   Mat                H;
 33:   Tao                tao; /* Tao solver context */
 34:   PetscBool          flg, test_lmvm = PETSC_FALSE;
 35:   PetscMPIInt        size; /* number of processes running */
 36:   AppCtx             user; /* user-defined application context */
 37:   TaoConvergedReason reason;
 38:   PetscInt           its, recycled_its = 0, oneshot_its = 0;

 40:   /* Initialize TAO and PETSc */
 42:   PetscInitialize(&argc, &argv, (char *)0, help);
 43:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 46:   /* Initialize problem parameters */
 47:   user.n       = 2;
 48:   user.alpha   = 99.0;
 49:   user.chained = PETSC_FALSE;
 50:   /* Check for command line arguments to override defaults */
 51:   PetscOptionsGetInt(NULL, NULL, "-n", &user.n, &flg);
 52:   PetscOptionsGetReal(NULL, NULL, "-alpha", &user.alpha, &flg);
 53:   PetscOptionsGetBool(NULL, NULL, "-chained", &user.chained, &flg);
 54:   PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg);

 56:   /* Allocate vectors for the solution and gradient */
 57:   VecCreateSeq(PETSC_COMM_SELF, user.n, &x);
 58:   MatCreateSeqBAIJ(PETSC_COMM_SELF, 2, user.n, user.n, 1, NULL, &H);

 60:   /* The TAO code begins here */

 62:   /* Create TAO solver with desired solution method */
 63:   TaoCreate(PETSC_COMM_SELF, &tao);
 64:   TaoSetType(tao, TAOLMVM);

 66:   /* Set solution vec and an initial guess */
 67:   VecSet(x, zero);
 68:   TaoSetSolution(tao, x);

 70:   /* Set routines for function, gradient, hessian evaluation */
 71:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, &user);
 72:   TaoSetHessian(tao, H, H, FormHessian, &user);

 74:   /* Check for TAO command line options */
 75:   TaoSetFromOptions(tao);

 77:   /* Solve the problem */
 78:   TaoSetTolerances(tao, 1.e-5, 0.0, 0.0);
 79:   TaoSetMaximumIterations(tao, 5);
 80:   TaoLMVMRecycle(tao, PETSC_TRUE);
 81:   reason = TAO_CONTINUE_ITERATING;
 82:   while (reason != TAO_CONVERGED_GATOL) {
 83:     TaoSolve(tao);
 84:     TaoGetConvergedReason(tao, &reason);
 85:     TaoGetIterationNumber(tao, &its);
 86:     recycled_its += its;
 87:     PetscPrintf(PETSC_COMM_SELF, "-----------------------\n");
 88:   }

 90:   /* Disable recycling and solve again! */
 91:   TaoSetMaximumIterations(tao, 100);
 92:   TaoLMVMRecycle(tao, PETSC_FALSE);
 93:   VecSet(x, zero);
 94:   TaoSolve(tao);
 95:   TaoGetConvergedReason(tao, &reason);
 97:   TaoGetIterationNumber(tao, &oneshot_its);
 98:   PetscPrintf(PETSC_COMM_SELF, "-----------------------\n");
 99:   PetscPrintf(PETSC_COMM_SELF, "recycled its: %" PetscInt_FMT " | oneshot its: %" PetscInt_FMT "\n", recycled_its, oneshot_its);

102:   TaoDestroy(&tao);
103:   VecDestroy(&x);
104:   MatDestroy(&H);

106:   PetscFinalize();
107:   return 0;
108: }

110: /* -------------------------------------------------------------------- */
111: /*
112:     FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X).

114:     Input Parameters:
115: .   tao  - the Tao context
116: .   X    - input vector
117: .   ptr  - optional user-defined context, as set by TaoSetFunctionGradient()

119:     Output Parameters:
120: .   G - vector containing the newly evaluated gradient
121: .   f - function value

123:     Note:
124:     Some optimization methods ask for the function and the gradient evaluation
125:     at the same time.  Evaluating both at once may be more efficient that
126:     evaluating each separately.
127: */
128: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
129: {
130:   AppCtx            *user = (AppCtx *)ptr;
131:   PetscInt           i, nn = user->n / 2;
132:   PetscReal          ff = 0, t1, t2, alpha = user->alpha;
133:   PetscScalar       *g;
134:   const PetscScalar *x;

137:   /* Get pointers to vector data */
138:   VecGetArrayRead(X, &x);
139:   VecGetArray(G, &g);

141:   /* Compute G(X) */
142:   if (user->chained) {
143:     g[0] = 0;
144:     for (i = 0; i < user->n - 1; i++) {
145:       t1 = x[i + 1] - x[i] * x[i];
146:       ff += PetscSqr(1 - x[i]) + alpha * t1 * t1;
147:       g[i] += -2 * (1 - x[i]) + 2 * alpha * t1 * (-2 * x[i]);
148:       g[i + 1] = 2 * alpha * t1;
149:     }
150:   } else {
151:     for (i = 0; i < nn; i++) {
152:       t1 = x[2 * i + 1] - x[2 * i] * x[2 * i];
153:       t2 = 1 - x[2 * i];
154:       ff += alpha * t1 * t1 + t2 * t2;
155:       g[2 * i]     = -4 * alpha * t1 * x[2 * i] - 2.0 * t2;
156:       g[2 * i + 1] = 2 * alpha * t1;
157:     }
158:   }

160:   /* Restore vectors */
161:   VecRestoreArrayRead(X, &x);
162:   VecRestoreArray(G, &g);
163:   *f = ff;

165:   PetscLogFlops(15.0 * nn);
166:   return 0;
167: }

169: /* ------------------------------------------------------------------- */
170: /*
171:    FormHessian - Evaluates Hessian matrix.

173:    Input Parameters:
174: .  tao   - the Tao context
175: .  x     - input vector
176: .  ptr   - optional user-defined context, as set by TaoSetHessian()

178:    Output Parameters:
179: .  H     - Hessian matrix

181:    Note:  Providing the Hessian may not be necessary.  Only some solvers
182:    require this matrix.
183: */
184: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
185: {
186:   AppCtx            *user = (AppCtx *)ptr;
187:   PetscInt           i, ind[2];
188:   PetscReal          alpha = user->alpha;
189:   PetscReal          v[2][2];
190:   const PetscScalar *x;
191:   PetscBool          assembled;

194:   /* Zero existing matrix entries */
195:   MatAssembled(H, &assembled);
196:   if (assembled) MatZeroEntries(H);

198:   /* Get a pointer to vector data */
199:   VecGetArrayRead(X, &x);

201:   /* Compute H(X) entries */
202:   if (user->chained) {
203:     MatZeroEntries(H);
204:     for (i = 0; i < user->n - 1; i++) {
205:       PetscScalar t1 = x[i + 1] - x[i] * x[i];
206:       v[0][0]        = 2 + 2 * alpha * (t1 * (-2) - 2 * x[i]);
207:       v[0][1]        = 2 * alpha * (-2 * x[i]);
208:       v[1][0]        = 2 * alpha * (-2 * x[i]);
209:       v[1][1]        = 2 * alpha * t1;
210:       ind[0]         = i;
211:       ind[1]         = i + 1;
212:       MatSetValues(H, 2, ind, 2, ind, v[0], ADD_VALUES);
213:     }
214:   } else {
215:     for (i = 0; i < user->n / 2; i++) {
216:       v[1][1] = 2 * alpha;
217:       v[0][0] = -4 * alpha * (x[2 * i + 1] - 3 * x[2 * i] * x[2 * i]) + 2;
218:       v[1][0] = v[0][1] = -4.0 * alpha * x[2 * i];
219:       ind[0]            = 2 * i;
220:       ind[1]            = 2 * i + 1;
221:       MatSetValues(H, 2, ind, 2, ind, v[0], INSERT_VALUES);
222:     }
223:   }
224:   VecRestoreArrayRead(X, &x);

226:   /* Assemble matrix */
227:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
228:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
229:   PetscLogFlops(9.0 * user->n / 2.0);
230:   return 0;
231: }

233: /*TEST

235:    build:
236:       requires: !complex

238:    test:
239:       args: -tao_type lmvm -tao_monitor
240:       requires: !single

242: TEST*/