Actual source code: ex1.c

      1: /* Program usage: mpiexec -n 2 ./ex1 [-help] [all TAO options] */

      3: /* ----------------------------------------------------------------------
      4: min f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
      5: s.t.  x0^2 + x1 - 2 = 0
      6:       0  <= x0^2 - x1 <= 1
      7:       -1 <= x0, x1 <= 2
      8: -->
      9:       g(x)  = 0
     10:       h(x) >= 0
     11:       -1 <= x0, x1 <= 2
     12: where
     13:       g(x) = x0^2 + x1 - 2
     14:       h(x) = [x0^2 - x1
     15:               1 -(x0^2 - x1)]
     16: ---------------------------------------------------------------------- */

     18: #include <petsctao.h>

     20: static char help[] = "Solves constrained optimization problem using pdipm.\n\
     21: Input parameters include:\n\
     22:   -tao_type pdipm    : sets Tao solver\n\
     23:   -no_eq             : removes the equaility constraints from the problem\n\
     24:   -init_view         : view initial object setup\n\
     25:   -snes_fd           : snes with finite difference Jacobian (needed for pdipm)\n\
     26:   -snes_compare_explicit : compare user Jacobian with finite difference Jacobian \n\
     27:   -tao_monitor_constraint_norm : convergence monitor with constraint norm \n\
     28:   -tao_monitor_solution : view exact solution at each iteration\n\
     29:   Note: external package MUMPS is required to run pdipm in parallel. This is designed for a maximum of 2 processors, the code will error if size > 2.\n";

     31: /*
     32:    User-defined application context - contains data needed by the application
     33: */
     34: typedef struct {
     35:   PetscInt   n;  /* Global length of x */
     36:   PetscInt   ne; /* Global number of equality constraints */
     37:   PetscInt   ni; /* Global number of inequality constraints */
     38:   PetscBool  noeqflag, initview;
     39:   Vec        x, xl, xu;
     40:   Vec        ce, ci, bl, bu, Xseq;
     41:   Mat        Ae, Ai, H;
     42:   VecScatter scat;
     43: } AppCtx;

     45: /* -------- User-defined Routines --------- */
     46: PetscErrorCode InitializeProblem(AppCtx *);
     47: PetscErrorCode DestroyProblem(AppCtx *);
     48: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
     49: PetscErrorCode FormPDIPMHessian(Tao, Vec, Mat, Mat, void *);
     50: PetscErrorCode FormInequalityConstraints(Tao, Vec, Vec, void *);
     51: PetscErrorCode FormEqualityConstraints(Tao, Vec, Vec, void *);
     52: PetscErrorCode FormInequalityJacobian(Tao, Vec, Mat, Mat, void *);
     53: PetscErrorCode FormEqualityJacobian(Tao, Vec, Mat, Mat, void *);

     55: int main(int argc, char **argv)
     56: {
     57:   Tao         tao;
     58:   KSP         ksp;
     59:   PC          pc;
     60:   AppCtx      user; /* application context */
     61:   Vec         x, G, CI, CE;
     62:   PetscMPIInt size;
     63:   TaoType     type;
     64:   PetscReal   f;
     65:   PetscBool   pdipm, mumps;

     67:   PetscFunctionBeginUser;
     68:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     69:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     70:   PetscCheck(size <= 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "More than 2 processors detected. Example written to use max of 2 processors.");

     72:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---- Constrained Problem -----\n"));
     73:   PetscCall(InitializeProblem(&user)); /* sets up problem, function below */

     75:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
     76:   PetscCall(TaoSetType(tao, TAOALMM));
     77:   PetscCall(TaoSetSolution(tao, user.x));
     78:   PetscCall(TaoSetVariableBounds(tao, user.xl, user.xu));
     79:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
     80:   PetscCall(TaoSetTolerances(tao, 1.e-4, 0.0, 0.0));
     81:   PetscCall(TaoSetConstraintTolerances(tao, 1.e-4, 0.0));
     82:   PetscCall(TaoSetMaximumFunctionEvaluations(tao, 1e6));
     83:   PetscCall(TaoSetFromOptions(tao));

     85:   if (!user.noeqflag) {
     86:     PetscCall(TaoSetEqualityConstraintsRoutine(tao, user.ce, FormEqualityConstraints, (void *)&user));
     87:     PetscCall(TaoSetJacobianEqualityRoutine(tao, user.Ae, user.Ae, FormEqualityJacobian, (void *)&user));
     88:   }
     89:   PetscCall(TaoSetInequalityConstraintsRoutine(tao, user.ci, FormInequalityConstraints, (void *)&user));
     90:   PetscCall(TaoSetJacobianInequalityRoutine(tao, user.Ai, user.Ai, FormInequalityJacobian, (void *)&user));

     92:   PetscCall(TaoGetType(tao, &type));
     93:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOPDIPM, &pdipm));
     94:   if (pdipm) {
     95:     /*
     96:       PDIPM produces an indefinite KKT matrix with some zeros along the diagonal
     97:       Inverting this indefinite matrix requires PETSc to be configured with MUMPS
     98:     */
     99:     PetscCall(PetscHasExternalPackage("mumps", &mumps));
    100:     PetscCheck(mumps, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAOPDIPM requires PETSc to be configured with MUMPS (--download-mumps)");
    101:     PetscCall(TaoGetKSP(tao, &ksp));
    102:     PetscCall(KSPGetPC(ksp, &pc));
    103:     PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS));
    104:     PetscCall(KSPSetType(ksp, KSPPREONLY));
    105:     PetscCall(PCSetType(pc, PCCHOLESKY));
    106:     PetscCall(KSPSetFromOptions(ksp));
    107:     PetscCall(TaoSetHessian(tao, user.H, user.H, FormPDIPMHessian, (void *)&user));
    108:   }

    110:   /* Print out an initial view of the problem */
    111:   if (user.initview) {
    112:     PetscCall(TaoSetUp(tao));
    113:     PetscCall(VecDuplicate(user.x, &G));
    114:     PetscCall(FormFunctionGradient(tao, user.x, &f, G, (void *)&user));
    115:     PetscCall(FormPDIPMHessian(tao, user.x, user.H, user.H, (void *)&user));
    116:     PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_WORLD));
    117:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial point X:\n"));
    118:     PetscCall(VecView(user.x, PETSC_VIEWER_STDOUT_WORLD));
    119:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial objective f(x) = %g\n", (double)f));
    120:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial gradient and Hessian:\n"));
    121:     PetscCall(VecView(G, PETSC_VIEWER_STDOUT_WORLD));
    122:     PetscCall(MatView(user.H, PETSC_VIEWER_STDOUT_WORLD));
    123:     PetscCall(VecDestroy(&G));
    124:     PetscCall(FormInequalityJacobian(tao, user.x, user.Ai, user.Ai, (void *)&user));
    125:     PetscCall(MatCreateVecs(user.Ai, NULL, &CI));
    126:     PetscCall(FormInequalityConstraints(tao, user.x, CI, (void *)&user));
    127:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial inequality constraints and Jacobian:\n"));
    128:     PetscCall(VecView(CI, PETSC_VIEWER_STDOUT_WORLD));
    129:     PetscCall(MatView(user.Ai, PETSC_VIEWER_STDOUT_WORLD));
    130:     PetscCall(VecDestroy(&CI));
    131:     if (!user.noeqflag) {
    132:       PetscCall(FormEqualityJacobian(tao, user.x, user.Ae, user.Ae, (void *)&user));
    133:       PetscCall(MatCreateVecs(user.Ae, NULL, &CE));
    134:       PetscCall(FormEqualityConstraints(tao, user.x, CE, (void *)&user));
    135:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nInitial equality constraints and Jacobian:\n"));
    136:       PetscCall(VecView(CE, PETSC_VIEWER_STDOUT_WORLD));
    137:       PetscCall(MatView(user.Ae, PETSC_VIEWER_STDOUT_WORLD));
    138:       PetscCall(VecDestroy(&CE));
    139:     }
    140:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
    141:     PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_WORLD));
    142:   }

    144:   PetscCall(TaoSolve(tao));
    145:   PetscCall(TaoGetSolution(tao, &x));
    146:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Found solution:\n"));
    147:   PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));

    149:   /* Free objects */
    150:   PetscCall(DestroyProblem(&user));
    151:   PetscCall(TaoDestroy(&tao));
    152:   PetscCall(PetscFinalize());
    153:   return PETSC_SUCCESS;
    154: }

    156: PetscErrorCode InitializeProblem(AppCtx *user)
    157: {
    158:   PetscMPIInt size;
    159:   PetscMPIInt rank;
    160:   PetscInt    nloc, neloc, niloc;

    162:   PetscFunctionBegin;
    163:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    164:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    165:   user->noeqflag = PETSC_FALSE;
    166:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_eq", &user->noeqflag, NULL));
    167:   user->initview = PETSC_FALSE;
    168:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-init_view", &user->initview, NULL));

    170:   /* Tell user the correct solution, not an error checking */
    171:   if (!user->noeqflag) {
    172:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution: f(1, 1) = -2\n"));
    173:   } else {
    174:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected solution (-no_eq): f(1.73205, 2) = -7.3923\n"));
    175:   }

    177:   /* create vector x and set initial values */
    178:   user->n = 2; /* global length */
    179:   nloc    = (size == 1) ? user->n : 1;
    180:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->x));
    181:   PetscCall(VecSetSizes(user->x, nloc, user->n));
    182:   PetscCall(VecSetFromOptions(user->x));
    183:   PetscCall(VecSet(user->x, 0.0));

    185:   /* create and set lower and upper bound vectors */
    186:   PetscCall(VecDuplicate(user->x, &user->xl));
    187:   PetscCall(VecDuplicate(user->x, &user->xu));
    188:   PetscCall(VecSet(user->xl, -1.0));
    189:   PetscCall(VecSet(user->xu, 2.0));

    191:   /* create scater to zero */
    192:   PetscCall(VecScatterCreateToZero(user->x, &user->scat, &user->Xseq));

    194:   user->ne = 1;
    195:   user->ni = 2;
    196:   neloc    = (rank == 0) ? user->ne : 0;
    197:   niloc    = (size == 1) ? user->ni : 1;

    199:   if (!user->noeqflag) {
    200:     PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ce)); /* a 1x1 vec for equality constraints */
    201:     PetscCall(VecSetSizes(user->ce, neloc, user->ne));
    202:     PetscCall(VecSetFromOptions(user->ce));
    203:     PetscCall(VecSetUp(user->ce));
    204:   }

    206:   PetscCall(VecCreate(PETSC_COMM_WORLD, &user->ci)); /* a 2x1 vec for inequality constraints */
    207:   PetscCall(VecSetSizes(user->ci, niloc, user->ni));
    208:   PetscCall(VecSetFromOptions(user->ci));
    209:   PetscCall(VecSetUp(user->ci));

    211:   /* nexn & nixn matrices for equally and inequalty constraints */
    212:   if (!user->noeqflag) {
    213:     PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ae));
    214:     PetscCall(MatSetSizes(user->Ae, neloc, nloc, user->ne, user->n));
    215:     PetscCall(MatSetFromOptions(user->Ae));
    216:     PetscCall(MatSetUp(user->Ae));
    217:   }

    219:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->Ai));
    220:   PetscCall(MatSetSizes(user->Ai, niloc, nloc, user->ni, user->n));
    221:   PetscCall(MatSetFromOptions(user->Ai));
    222:   PetscCall(MatSetUp(user->Ai));

    224:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user->H));
    225:   PetscCall(MatSetSizes(user->H, nloc, nloc, user->n, user->n));
    226:   PetscCall(MatSetFromOptions(user->H));
    227:   PetscCall(MatSetUp(user->H));
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: PetscErrorCode DestroyProblem(AppCtx *user)
    232: {
    233:   PetscFunctionBegin;
    234:   if (!user->noeqflag) PetscCall(MatDestroy(&user->Ae));
    235:   PetscCall(MatDestroy(&user->Ai));
    236:   PetscCall(MatDestroy(&user->H));

    238:   PetscCall(VecDestroy(&user->x));
    239:   if (!user->noeqflag) PetscCall(VecDestroy(&user->ce));
    240:   PetscCall(VecDestroy(&user->ci));
    241:   PetscCall(VecDestroy(&user->xl));
    242:   PetscCall(VecDestroy(&user->xu));
    243:   PetscCall(VecDestroy(&user->Xseq));
    244:   PetscCall(VecScatterDestroy(&user->scat));
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }

    248: /* Evaluate
    249:    f(x) = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1)
    250:    G = grad f = [2*(x0 - 2) - 2, 2*(x1 - 2) - 2] = 2*X - 6
    251: */
    252: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ctx)
    253: {
    254:   const PetscScalar *x;
    255:   MPI_Comm           comm;
    256:   PetscMPIInt        rank;
    257:   PetscReal          fin;
    258:   AppCtx            *user = (AppCtx *)ctx;
    259:   Vec                Xseq = user->Xseq;
    260:   VecScatter         scat = user->scat;

    262:   PetscFunctionBegin;
    263:   /* f = (x0 - 2)^2 + (x1 - 2)^2 - 2*(x0 + x1) */
    264:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
    265:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    267:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
    268:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

    270:   fin = 0.0;
    271:   if (rank == 0) {
    272:     PetscCall(VecGetArrayRead(Xseq, &x));
    273:     fin = (x[0] - 2.0) * (x[0] - 2.0) + (x[1] - 2.0) * (x[1] - 2.0) - 2.0 * (x[0] + x[1]);
    274:     PetscCall(VecRestoreArrayRead(Xseq, &x));
    275:   }
    276:   PetscCallMPI(MPIU_Allreduce(&fin, f, 1, MPIU_REAL, MPIU_SUM, comm));

    278:   /* G = 2*X - 6 */
    279:   PetscCall(VecSet(G, -6.0));
    280:   PetscCall(VecAXPY(G, 2.0, X));
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }

    284: /* Evaluate PDIPM Hessian, see Eqn(22) in http://doi.org/10.1049/gtd2.12708
    285:    H = Wxx = fxx        + Jacobian(grad g^T*DE)   - Jacobian(grad h^T*DI)]
    286:            = fxx        + Jacobin([2*x0; 1]de[0]) - Jacobian([2*x0, -2*x0; -1, 1][di[0] di[1]]^T)
    287:            = [2 0; 0 2] + [2*de[0]  0;      0  0] - [2*di[0]-2*di[1], 0; 0, 0]
    288:            = [ 2*(1+de[0]-di[0]+di[1]), 0;
    289:                           0,            2]
    290: */
    291: PetscErrorCode FormPDIPMHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
    292: {
    293:   AppCtx            *user = (AppCtx *)ctx;
    294:   Vec                DE, DI;
    295:   const PetscScalar *de, *di;
    296:   PetscInt           zero = 0, one = 1;
    297:   PetscScalar        two = 2.0;
    298:   PetscScalar        val = 0.0;
    299:   Vec                Deseq, Diseq;
    300:   VecScatter         Descat, Discat;
    301:   PetscMPIInt        rank;
    302:   MPI_Comm           comm;

    304:   PetscFunctionBegin;
    305:   PetscCall(TaoGetDualVariables(tao, &DE, &DI));

    307:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
    308:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    310:   if (!user->noeqflag) {
    311:     PetscCall(VecScatterCreateToZero(DE, &Descat, &Deseq));
    312:     PetscCall(VecScatterBegin(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
    313:     PetscCall(VecScatterEnd(Descat, DE, Deseq, INSERT_VALUES, SCATTER_FORWARD));
    314:   }
    315:   PetscCall(VecScatterCreateToZero(DI, &Discat, &Diseq));
    316:   PetscCall(VecScatterBegin(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));
    317:   PetscCall(VecScatterEnd(Discat, DI, Diseq, INSERT_VALUES, SCATTER_FORWARD));

    319:   if (rank == 0) {
    320:     if (!user->noeqflag) { PetscCall(VecGetArrayRead(Deseq, &de)); /* places equality constraint dual into array */ }
    321:     PetscCall(VecGetArrayRead(Diseq, &di)); /* places inequality constraint dual into array */

    323:     if (!user->noeqflag) {
    324:       val = 2.0 * (1 + de[0] - di[0] + di[1]);
    325:       PetscCall(VecRestoreArrayRead(Deseq, &de));
    326:       PetscCall(VecRestoreArrayRead(Diseq, &di));
    327:     } else {
    328:       val = 2.0 * (1 - di[0] + di[1]);
    329:     }
    330:     PetscCall(VecRestoreArrayRead(Diseq, &di));
    331:     PetscCall(MatSetValues(H, 1, &zero, 1, &zero, &val, INSERT_VALUES));
    332:     PetscCall(MatSetValues(H, 1, &one, 1, &one, &two, INSERT_VALUES));
    333:   }
    334:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
    335:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
    336:   if (!user->noeqflag) {
    337:     PetscCall(VecScatterDestroy(&Descat));
    338:     PetscCall(VecDestroy(&Deseq));
    339:   }
    340:   PetscCall(VecScatterDestroy(&Discat));
    341:   PetscCall(VecDestroy(&Diseq));
    342:   PetscFunctionReturn(PETSC_SUCCESS);
    343: }

    345: /* Evaluate
    346:    h = [ x0^2 - x1;
    347:          1 -(x0^2 - x1)]
    348: */
    349: PetscErrorCode FormInequalityConstraints(Tao tao, Vec X, Vec CI, void *ctx)
    350: {
    351:   const PetscScalar *x;
    352:   PetscScalar        ci;
    353:   MPI_Comm           comm;
    354:   PetscMPIInt        rank;
    355:   AppCtx            *user = (AppCtx *)ctx;
    356:   Vec                Xseq = user->Xseq;
    357:   VecScatter         scat = user->scat;

    359:   PetscFunctionBegin;
    360:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
    361:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    363:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
    364:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

    366:   if (rank == 0) {
    367:     PetscCall(VecGetArrayRead(Xseq, &x));
    368:     ci = x[0] * x[0] - x[1];
    369:     PetscCall(VecSetValue(CI, 0, ci, INSERT_VALUES));
    370:     ci = -x[0] * x[0] + x[1] + 1.0;
    371:     PetscCall(VecSetValue(CI, 1, ci, INSERT_VALUES));
    372:     PetscCall(VecRestoreArrayRead(Xseq, &x));
    373:   }
    374:   PetscCall(VecAssemblyBegin(CI));
    375:   PetscCall(VecAssemblyEnd(CI));
    376:   PetscFunctionReturn(PETSC_SUCCESS);
    377: }

    379: /* Evaluate
    380:    g = [ x0^2 + x1 - 2]
    381: */
    382: PetscErrorCode FormEqualityConstraints(Tao tao, Vec X, Vec CE, void *ctx)
    383: {
    384:   const PetscScalar *x;
    385:   PetscScalar        ce;
    386:   MPI_Comm           comm;
    387:   PetscMPIInt        rank;
    388:   AppCtx            *user = (AppCtx *)ctx;
    389:   Vec                Xseq = user->Xseq;
    390:   VecScatter         scat = user->scat;

    392:   PetscFunctionBegin;
    393:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
    394:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    396:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
    397:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

    399:   if (rank == 0) {
    400:     PetscCall(VecGetArrayRead(Xseq, &x));
    401:     ce = x[0] * x[0] + x[1] - 2.0;
    402:     PetscCall(VecSetValue(CE, 0, ce, INSERT_VALUES));
    403:     PetscCall(VecRestoreArrayRead(Xseq, &x));
    404:   }
    405:   PetscCall(VecAssemblyBegin(CE));
    406:   PetscCall(VecAssemblyEnd(CE));
    407:   PetscFunctionReturn(PETSC_SUCCESS);
    408: }

    410: /*
    411:   grad h = [  2*x0, -1;
    412:              -2*x0,  1]
    413: */
    414: PetscErrorCode FormInequalityJacobian(Tao tao, Vec X, Mat JI, Mat JIpre, void *ctx)
    415: {
    416:   AppCtx            *user = (AppCtx *)ctx;
    417:   PetscInt           zero = 0, one = 1, cols[2];
    418:   PetscScalar        vals[2];
    419:   const PetscScalar *x;
    420:   Vec                Xseq = user->Xseq;
    421:   VecScatter         scat = user->scat;
    422:   MPI_Comm           comm;
    423:   PetscMPIInt        rank;

    425:   PetscFunctionBegin;
    426:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
    427:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    428:   PetscCall(VecScatterBegin(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));
    429:   PetscCall(VecScatterEnd(scat, X, Xseq, INSERT_VALUES, SCATTER_FORWARD));

    431:   PetscCall(VecGetArrayRead(Xseq, &x));
    432:   if (rank == 0) {
    433:     cols[0] = 0;
    434:     cols[1] = 1;
    435:     vals[0] = 2 * x[0];
    436:     vals[1] = -1.0;
    437:     PetscCall(MatSetValues(JI, 1, &zero, 2, cols, vals, INSERT_VALUES));
    438:     vals[0] = -2 * x[0];
    439:     vals[1] = 1.0;
    440:     PetscCall(MatSetValues(JI, 1, &one, 2, cols, vals, INSERT_VALUES));
    441:   }
    442:   PetscCall(VecRestoreArrayRead(Xseq, &x));
    443:   PetscCall(MatAssemblyBegin(JI, MAT_FINAL_ASSEMBLY));
    444:   PetscCall(MatAssemblyEnd(JI, MAT_FINAL_ASSEMBLY));
    445:   PetscFunctionReturn(PETSC_SUCCESS);
    446: }

    448: /*
    449:   grad g = [2*x0, 1.0]
    450: */
    451: PetscErrorCode FormEqualityJacobian(Tao tao, Vec X, Mat JE, Mat JEpre, void *ctx)
    452: {
    453:   PetscInt           zero = 0, cols[2];
    454:   PetscScalar        vals[2];
    455:   const PetscScalar *x;
    456:   PetscMPIInt        rank;
    457:   MPI_Comm           comm;

    459:   PetscFunctionBegin;
    460:   PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
    461:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    463:   if (rank == 0) {
    464:     PetscCall(VecGetArrayRead(X, &x));
    465:     cols[0] = 0;
    466:     cols[1] = 1;
    467:     vals[0] = 2 * x[0];
    468:     vals[1] = 1.0;
    469:     PetscCall(MatSetValues(JE, 1, &zero, 2, cols, vals, INSERT_VALUES));
    470:     PetscCall(VecRestoreArrayRead(X, &x));
    471:   }
    472:   PetscCall(MatAssemblyBegin(JE, MAT_FINAL_ASSEMBLY));
    473:   PetscCall(MatAssemblyEnd(JE, MAT_FINAL_ASSEMBLY));
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: /*TEST

    479:    build:
    480:       requires: !complex !defined(PETSC_USE_CXX)

    482:    test:
    483:       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm -tao_pdipm_kkt_shift_pd
    484:       requires: mumps !single
    485:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    487:    test:
    488:       suffix: pdipm_2
    489:       requires: mumps
    490:       nsize: 2
    491:       args: -tao_converged_reason -tao_gatol 1.e-6 -tao_type pdipm
    492:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    494:    test:
    495:       suffix: 2
    496:       args: -tao_converged_reason
    497:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    499:    test:
    500:       suffix: 3
    501:       args: -tao_converged_reason -no_eq
    502:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    504:    test:
    505:       suffix: 4
    506:       args: -tao_converged_reason -tao_almm_type classic
    507:       requires: !single
    508:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    510:    test:
    511:       suffix: 5
    512:       args: -tao_converged_reason -tao_almm_type classic -no_eq -tao_almm_subsolver_tao_max_it 100
    513:       requires: !single !defined(PETSCTEST_VALGRIND)
    514:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    516:    test:
    517:       suffix: 6
    518:       args: -tao_converged_reason -tao_almm_subsolver_tao_type bqnktr
    519:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    521:    test:
    522:       suffix: 7
    523:       args: -tao_converged_reason -tao_almm_subsolver_tao_type bncg
    524:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    526:    test:
    527:       suffix: 8
    528:       nsize: 2
    529:       args: -tao_converged_reason
    530:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    532:    test:
    533:       suffix: 9
    534:       nsize: 2
    535:       args: -tao_converged_reason -vec_type cuda -mat_type aijcusparse
    536:       requires: cuda
    537:       filter: sed  -e "s/CONVERGED_GATOL iterations *[0-9]\{1,\}/CONVERGED_GATOL/g"

    539: TEST*/