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: int 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 */
 63:   PetscFunctionBeginUser;
 64:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 65:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 66:   /* Specify default parameters for the problem, check for command-line overrides */
 67:   PetscCall(PetscStrncpy(user.name, "HS21", sizeof(user.name)));
 68:   PetscCall(PetscOptionsGetString(NULL, NULL, "-cutername", user.name, sizeof(user.name), &flg));

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

 77:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
 78:   PetscCall(TaoSetType(tao, TAOIPM));
 79:   PetscCall(TaoSetSolution(tao, x));
 80:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
 81:   PetscCall(TaoSetEqualityConstraintsRoutine(tao, ceq, FormEqualityConstraints, (void *)&user));
 82:   PetscCall(TaoSetInequalityConstraintsRoutine(tao, cin, FormInequalityConstraints, (void *)&user));
 83:   PetscCall(TaoSetInequalityBounds(tao, user.bin, NULL));
 84:   PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Aeq, user.Aeq, FormEqualityJacobian, (void *)&user));
 85:   PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ain, user.Ain, FormInequalityJacobian, (void *)&user));
 86:   PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
 87:   PetscCall(TaoGetKSP(tao, &ksp));
 88:   PetscCall(KSPGetPC(ksp, &pc));
 89:   PetscCall(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:   PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU));
 95:   PetscCall(KSPSetType(ksp, KSPPREONLY));
 96:   PetscCall(TaoSetTolerances(tao, 0, 0, 0));

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

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

113:   PetscCall(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:   PetscFunctionBegin;
127:   comm = PETSC_COMM_WORLD;
128:   PetscCall(PetscStrncpy(filebase, user->name, sizeof(filebase)));
129:   PetscCall(PetscStrlcat(filebase, "/", sizeof(filebase)));
130:   PetscCall(PetscStrncpy(filename, filebase, sizeof(filename)));
131:   PetscCall(PetscStrlcat(filename, "f", sizeof(filename)));
132:   PetscCall(PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader));

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

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

145:   PetscCall(MatCreate(comm, &user->H));
146:   PetscCall(MatSetSizes(user->H, PETSC_DECIDE, PETSC_DECIDE, nrows, nrows));
147:   PetscCall(MatLoad(user->H, loader));
148:   PetscCall(PetscViewerDestroy(&loader));
149:   PetscCall(MatGetSize(user->H, &nrows, &ncols));
150:   PetscCheck(nrows == user->n, comm, PETSC_ERR_ARG_SIZ, "H: nrows != n");
151:   PetscCheck(ncols == user->n, comm, PETSC_ERR_ARG_SIZ, "H: ncols != n");
152:   PetscCall(MatSetFromOptions(user->H));

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

165:   PetscCall(PetscStrncpy(filename, filebase, sizeof(filename)));
166:   PetscCall(PetscStrlcat(filename, "Beq", sizeof(filename)));
167:   PetscCall(PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &loader));
168:   PetscCall(VecCreate(comm, &user->beq));
169:   PetscCall(VecLoad(user->beq, loader));
170:   PetscCall(PetscViewerDestroy(&loader));
171:   PetscCall(VecGetSize(user->beq, &nrows));
172:   PetscCheck(nrows == user->me, comm, PETSC_ERR_ARG_SIZ, "Aeq nrows != Beq n");
173:   PetscCall(VecSetFromOptions(user->beq));

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

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

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

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

199: PetscErrorCode DestroyProblem(AppCtx *user)
200: {
201:   PetscFunctionBegin;
202:   PetscCall(MatDestroy(&user->H));
203:   PetscCall(MatDestroy(&user->Aeq));
204:   PetscCall(MatDestroy(&user->Ain));
205:   PetscCall(VecDestroy(&user->beq));
206:   PetscCall(VecDestroy(&user->bin));
207:   PetscCall(VecDestroy(&user->d));
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

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

216:   PetscFunctionBegin;
217:   PetscCall(MatMult(user->H, x, g));
218:   PetscCall(VecDot(x, g, &xtHx));
219:   PetscCall(VecDot(x, user->d, f));
220:   *f += 0.5 * xtHx;
221:   PetscCall(VecAXPY(g, 1.0, user->d));
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

225: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
226: {
227:   PetscFunctionBegin;
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
232: {
233:   AppCtx *user = (AppCtx *)ctx;

235:   PetscFunctionBegin;
236:   PetscCall(MatMult(user->Ain, x, ci));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce, void *ctx)
241: {
242:   AppCtx *user = (AppCtx *)ctx;

244:   PetscFunctionBegin;
245:   PetscCall(MatMult(user->Aeq, x, ce));
246:   PetscCall(VecAXPY(ce, -1.0, user->beq));
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

250: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx)
251: {
252:   PetscFunctionBegin;
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

256: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
257: {
258:   PetscFunctionBegin;
259:   PetscFunctionReturn(PETSC_SUCCESS);
260: }

262: /*TEST

264:    build:
265:       requires: !complex

267:    test:
268:       requires: superlu
269:       localrunfiles: HS21

271: TEST*/