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, NULL, 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, PetscCtx 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, PetscCtx ctx)
    226: {
    227:   PetscFunctionBegin;
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, PetscCtx 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, PetscCtx 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, PetscCtx ctx)
    251: {
    252:   PetscFunctionBegin;
    253:   PetscFunctionReturn(PETSC_SUCCESS);
    254: }

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

    262: /*TEST

    264:    build:
    265:       requires: !complex

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

    271: TEST*/