Actual source code: maros.c

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

  3: /* ----------------------------------------------------------------------
  4: TODO Explain maros example
  5: ---------------------------------------------------------------------- */

  7: #include <petsctao.h>

  9: static char help[] = "";

 11: /*
 12:    User-defined application context - contains data needed by the
 13:    application-provided call-back routines, FormFunction(),
 14:    FormGradient(), and FormHessian().
 15: */

 17: /*
 18:    x,d in R^n
 19:    f in R
 20:    bin in R^mi
 21:    beq in R^me
 22:    Aeq in R^(me x n)
 23:    Ain in R^(mi x n)
 24:    H in R^(n x n)
 25:    min f=(1/2)*x'*H*x + d'*x
 26:    s.t.  Aeq*x == beq
 27:          Ain*x >= bin
 28: */
 29: typedef struct {
 30:   char     name[32];
 31:   PetscInt n;  /* Length x */
 32:   PetscInt me; /* number of equality constraints */
 33:   PetscInt mi; /* number of inequality constraints */
 34:   PetscInt m;  /* me+mi */
 35:   Mat      Aeq, Ain, H;
 36:   Vec      beq, bin, d;
 37: } AppCtx;

 39: /* -------- User-defined Routines --------- */

 41: PetscErrorCode InitializeProblem(AppCtx *);
 42: PetscErrorCode DestroyProblem(AppCtx *);
 43: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 44: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
 45: PetscErrorCode FormInequalityConstraints(Tao, Vec, Vec, void *);
 46: PetscErrorCode FormEqualityConstraints(Tao, Vec, Vec, void *);
 47: PetscErrorCode FormInequalityJacobian(Tao, Vec, Mat, Mat, void *);
 48: PetscErrorCode FormEqualityJacobian(Tao, Vec, Mat, Mat, void *);

 50: PetscErrorCode main(int argc, char **argv)
 51: {
 52:   PetscMPIInt        size;
 53:   Vec                x; /* solution */
 54:   KSP                ksp;
 55:   PC                 pc;
 56:   Vec                ceq, cin;
 57:   PetscBool          flg; /* A return value when checking for use options */
 58:   Tao                tao; /* Tao solver context */
 59:   TaoConvergedReason reason;
 60:   AppCtx             user; /* application context */

 62:   /* Initialize TAO,PETSc */
 64:   PetscInitialize(&argc, &argv, (char *)0, help);
 65:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 66:   /* Specify default parameters for the problem, check for command-line overrides */
 67:   PetscStrncpy(user.name, "HS21", sizeof(user.name));
 68:   PetscOptionsGetString(NULL, NULL, "-cutername", user.name, sizeof(user.name), &flg);

 70:   PetscPrintf(PETSC_COMM_WORLD, "\n---- MAROS Problem %s -----\n", user.name);
 71:   InitializeProblem(&user);
 72:   VecDuplicate(user.d, &x);
 73:   VecDuplicate(user.beq, &ceq);
 74:   VecDuplicate(user.bin, &cin);
 75:   VecSet(x, 1.0);

 77:   TaoCreate(PETSC_COMM_WORLD, &tao);
 78:   TaoSetType(tao, TAOIPM);
 79:   TaoSetSolution(tao, x);
 80:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
 81:   TaoSetEqualityConstraintsRoutine(tao, ceq, FormEqualityConstraints, (void *)&user);
 82:   TaoSetInequalityConstraintsRoutine(tao, cin, FormInequalityConstraints, (void *)&user);
 83:   TaoSetInequalityBounds(tao, user.bin, NULL);
 84:   TaoSetJacobianEqualityRoutine(tao, user.Aeq, user.Aeq, FormEqualityJacobian, (void *)&user);
 85:   TaoSetJacobianInequalityRoutine(tao, user.Ain, user.Ain, FormInequalityJacobian, (void *)&user);
 86:   TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);
 87:   TaoGetKSP(tao, &ksp);
 88:   KSPGetPC(ksp, &pc);
 89:   PCSetType(pc, PCLU);
 90:   /*
 91:       This algorithm produces matrices with zeros along the diagonal therefore we need to use
 92:     SuperLU which does partial pivoting
 93:   */
 94:   PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU);
 95:   KSPSetType(ksp, KSPPREONLY);
 96:   TaoSetTolerances(tao, 0, 0, 0);

 98:   TaoSetFromOptions(tao);
 99:   TaoSolve(tao);
100:   TaoGetConvergedReason(tao, &reason);
101:   if (reason < 0) {
102:     PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n", TaoConvergedReasons[reason]);
103:   } else {
104:     PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n", TaoConvergedReasons[reason]);
105:   }

107:   DestroyProblem(&user);
108:   VecDestroy(&x);
109:   VecDestroy(&ceq);
110:   VecDestroy(&cin);
111:   TaoDestroy(&tao);

113:   PetscFinalize();
114:   return 0;
115: }

117: PetscErrorCode InitializeProblem(AppCtx *user)
118: {
119:   PetscViewer loader;
120:   MPI_Comm    comm;
121:   PetscInt    nrows, ncols, i;
122:   PetscScalar one = 1.0;
123:   char        filebase[128];
124:   char        filename[128];

126:   comm = PETSC_COMM_WORLD;
127:   PetscStrncpy(filebase, user->name, sizeof(filebase));
128:   PetscStrlcat(filebase, "/", sizeof(filebase));
129:   PetscStrncpy(filename, filebase, sizeof(filename));
130:   PetscStrlcat(filename, "f", sizeof(filename));
131:   PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);

133:   VecCreate(comm, &user->d);
134:   VecLoad(user->d, loader);
135:   PetscViewerDestroy(&loader);
136:   VecGetSize(user->d, &nrows);
137:   VecSetFromOptions(user->d);
138:   user->n = nrows;

140:   PetscStrncpy(filename, filebase, sizeof(filename));
141:   PetscStrlcat(filename, "H", sizeof(filename));
142:   PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);

144:   MatCreate(comm, &user->H);
145:   MatSetSizes(user->H, PETSC_DECIDE, PETSC_DECIDE, nrows, nrows);
146:   MatLoad(user->H, loader);
147:   PetscViewerDestroy(&loader);
148:   MatGetSize(user->H, &nrows, &ncols);
151:   MatSetFromOptions(user->H);

153:   PetscStrncpy(filename, filebase, sizeof(filename));
154:   PetscStrlcat(filename, "Aeq", sizeof(filename));
155:   PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);
156:   MatCreate(comm, &user->Aeq);
157:   MatLoad(user->Aeq, loader);
158:   PetscViewerDestroy(&loader);
159:   MatGetSize(user->Aeq, &nrows, &ncols);
161:   MatSetFromOptions(user->Aeq);
162:   user->me = nrows;

164:   PetscStrncpy(filename, filebase, sizeof(filename));
165:   PetscStrlcat(filename, "Beq", sizeof(filename));
166:   PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader);
167:   VecCreate(comm, &user->beq);
168:   VecLoad(user->beq, loader);
169:   PetscViewerDestroy(&loader);
170:   VecGetSize(user->beq, &nrows);
172:   VecSetFromOptions(user->beq);

174:   user->mi = user->n;
175:   /* Ain = eye(n,n) */
176:   MatCreate(comm, &user->Ain);
177:   MatSetType(user->Ain, MATAIJ);
178:   MatSetSizes(user->Ain, PETSC_DECIDE, PETSC_DECIDE, user->mi, user->mi);

180:   MatMPIAIJSetPreallocation(user->Ain, 1, NULL, 0, NULL);
181:   MatSeqAIJSetPreallocation(user->Ain, 1, NULL);

183:   for (i = 0; i < user->mi; i++) MatSetValues(user->Ain, 1, &i, 1, &i, &one, INSERT_VALUES);
184:   MatAssemblyBegin(user->Ain, MAT_FINAL_ASSEMBLY);
185:   MatAssemblyEnd(user->Ain, MAT_FINAL_ASSEMBLY);
186:   MatSetFromOptions(user->Ain);

188:   /* bin = [0,0 ... 0]' */
189:   VecCreate(comm, &user->bin);
190:   VecSetType(user->bin, VECMPI);
191:   VecSetSizes(user->bin, PETSC_DECIDE, user->mi);
192:   VecSet(user->bin, 0.0);
193:   VecSetFromOptions(user->bin);
194:   user->m = user->me + user->mi;
195:   return 0;
196: }

198: PetscErrorCode DestroyProblem(AppCtx *user)
199: {
200:   MatDestroy(&user->H);
201:   MatDestroy(&user->Aeq);
202:   MatDestroy(&user->Ain);
203:   VecDestroy(&user->beq);
204:   VecDestroy(&user->bin);
205:   VecDestroy(&user->d);
206:   return 0;
207: }

209: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
210: {
211:   AppCtx     *user = (AppCtx *)ctx;
212:   PetscScalar xtHx;

214:   MatMult(user->H, x, g);
215:   VecDot(x, g, &xtHx);
216:   VecDot(x, user->d, f);
217:   *f += 0.5 * xtHx;
218:   VecAXPY(g, 1.0, user->d);
219:   return 0;
220: }

222: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
223: {
224:   return 0;
225: }

227: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
228: {
229:   AppCtx *user = (AppCtx *)ctx;

231:   MatMult(user->Ain, x, ci);
232:   return 0;
233: }

235: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce, void *ctx)
236: {
237:   AppCtx *user = (AppCtx *)ctx;

239:   MatMult(user->Aeq, x, ce);
240:   VecAXPY(ce, -1.0, user->beq);
241:   return 0;
242: }

244: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
245: {
246:   return 0;
247: }

249: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
250: {
251:   return 0;
252: }

254: /*TEST

256:    build:
257:       requires: !complex

259:    test:
260:       requires: superlu
261:       localrunfiles: HS21

263: TEST*/