Actual source code: lcl.c

      1: #include <../src/tao/pde_constrained/impls/lcl/lcl.h>
      2: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch, Vec, PetscReal *, Vec, void *);
      3: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch, Vec, PetscReal *, Vec, void *);
      4: static PetscErrorCode LCLScatter(TAO_LCL *, Vec, Vec, Vec);
      5: static PetscErrorCode LCLGather(TAO_LCL *, Vec, Vec, Vec);

      7: static PetscErrorCode TaoDestroy_LCL(Tao tao)
      8: {
      9:   TAO_LCL *lclP = (TAO_LCL *)tao->data;

     11:   PetscFunctionBegin;
     12:   if (tao->setupcalled) {
     13:     PetscCall(MatDestroy(&lclP->R));
     14:     PetscCall(VecDestroy(&lclP->lambda));
     15:     PetscCall(VecDestroy(&lclP->lambda0));
     16:     PetscCall(VecDestroy(&lclP->WL));
     17:     PetscCall(VecDestroy(&lclP->W));
     18:     PetscCall(VecDestroy(&lclP->X0));
     19:     PetscCall(VecDestroy(&lclP->G0));
     20:     PetscCall(VecDestroy(&lclP->GL));
     21:     PetscCall(VecDestroy(&lclP->GAugL));
     22:     PetscCall(VecDestroy(&lclP->dbar));
     23:     PetscCall(VecDestroy(&lclP->U));
     24:     PetscCall(VecDestroy(&lclP->U0));
     25:     PetscCall(VecDestroy(&lclP->V));
     26:     PetscCall(VecDestroy(&lclP->V0));
     27:     PetscCall(VecDestroy(&lclP->V1));
     28:     PetscCall(VecDestroy(&lclP->GU));
     29:     PetscCall(VecDestroy(&lclP->GV));
     30:     PetscCall(VecDestroy(&lclP->GU0));
     31:     PetscCall(VecDestroy(&lclP->GV0));
     32:     PetscCall(VecDestroy(&lclP->GL_U));
     33:     PetscCall(VecDestroy(&lclP->GL_V));
     34:     PetscCall(VecDestroy(&lclP->GAugL_U));
     35:     PetscCall(VecDestroy(&lclP->GAugL_V));
     36:     PetscCall(VecDestroy(&lclP->GL_U0));
     37:     PetscCall(VecDestroy(&lclP->GL_V0));
     38:     PetscCall(VecDestroy(&lclP->GAugL_U0));
     39:     PetscCall(VecDestroy(&lclP->GAugL_V0));
     40:     PetscCall(VecDestroy(&lclP->DU));
     41:     PetscCall(VecDestroy(&lclP->DV));
     42:     PetscCall(VecDestroy(&lclP->WU));
     43:     PetscCall(VecDestroy(&lclP->WV));
     44:     PetscCall(VecDestroy(&lclP->g1));
     45:     PetscCall(VecDestroy(&lclP->g2));
     46:     PetscCall(VecDestroy(&lclP->con1));

     48:     PetscCall(VecDestroy(&lclP->r));
     49:     PetscCall(VecDestroy(&lclP->s));

     51:     PetscCall(ISDestroy(&tao->state_is));
     52:     PetscCall(ISDestroy(&tao->design_is));

     54:     PetscCall(VecScatterDestroy(&lclP->state_scatter));
     55:     PetscCall(VecScatterDestroy(&lclP->design_scatter));
     56:   }
     57:   PetscCall(MatDestroy(&lclP->R));
     58:   PetscCall(KSPDestroy(&tao->ksp));
     59:   PetscCall(PetscFree(tao->data));
     60:   PetscFunctionReturn(PETSC_SUCCESS);
     61: }

     63: static PetscErrorCode TaoSetFromOptions_LCL(Tao tao, PetscOptionItems PetscOptionsObject)
     64: {
     65:   TAO_LCL *lclP = (TAO_LCL *)tao->data;

     67:   PetscFunctionBegin;
     68:   PetscOptionsHeadBegin(PetscOptionsObject, "Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");
     69:   PetscCall(PetscOptionsReal("-tao_lcl_eps1", "epsilon 1 tolerance", "", lclP->eps1, &lclP->eps1, NULL));
     70:   PetscCall(PetscOptionsReal("-tao_lcl_eps2", "epsilon 2 tolerance", "", lclP->eps2, &lclP->eps2, NULL));
     71:   PetscCall(PetscOptionsReal("-tao_lcl_rho0", "init value for rho", "", lclP->rho0, &lclP->rho0, NULL));
     72:   PetscCall(PetscOptionsReal("-tao_lcl_rhomax", "max value for rho", "", lclP->rhomax, &lclP->rhomax, NULL));
     73:   lclP->phase2_niter = 1;
     74:   PetscCall(PetscOptionsInt("-tao_lcl_phase2_niter", "Number of phase 2 iterations in LCL algorithm", "", lclP->phase2_niter, &lclP->phase2_niter, NULL));
     75:   lclP->verbose = PETSC_FALSE;
     76:   PetscCall(PetscOptionsBool("-tao_lcl_verbose", "Print verbose output", "", lclP->verbose, &lclP->verbose, NULL));
     77:   lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
     78:   PetscCall(PetscOptionsReal("-tao_lcl_tola", "Tolerance for first forward solve", "", lclP->tau[0], &lclP->tau[0], NULL));
     79:   PetscCall(PetscOptionsReal("-tao_lcl_tolb", "Tolerance for first adjoint solve", "", lclP->tau[1], &lclP->tau[1], NULL));
     80:   PetscCall(PetscOptionsReal("-tao_lcl_tolc", "Tolerance for second forward solve", "", lclP->tau[2], &lclP->tau[2], NULL));
     81:   PetscCall(PetscOptionsReal("-tao_lcl_told", "Tolerance for second adjoint solve", "", lclP->tau[3], &lclP->tau[3], NULL));
     82:   PetscOptionsHeadEnd();
     83:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
     84:   PetscCall(MatSetFromOptions(lclP->R));
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: static PetscErrorCode TaoView_LCL(Tao tao, PetscViewer viewer)
     89: {
     90:   return PETSC_SUCCESS;
     91: }

     93: static PetscErrorCode TaoSetup_LCL(Tao tao)
     94: {
     95:   TAO_LCL *lclP = (TAO_LCL *)tao->data;
     96:   PetscInt lo, hi, nlocalstate, nlocaldesign;
     97:   IS       is_state, is_design;

     99:   PetscFunctionBegin;
    100:   PetscCheck(tao->state_is, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "LCL Solver requires an initial state index set -- use TaoSetStateIS()");
    101:   PetscCall(VecDuplicate(tao->solution, &tao->gradient));
    102:   PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
    103:   PetscCall(VecDuplicate(tao->solution, &lclP->W));
    104:   PetscCall(VecDuplicate(tao->solution, &lclP->X0));
    105:   PetscCall(VecDuplicate(tao->solution, &lclP->G0));
    106:   PetscCall(VecDuplicate(tao->solution, &lclP->GL));
    107:   PetscCall(VecDuplicate(tao->solution, &lclP->GAugL));

    109:   PetscCall(VecDuplicate(tao->constraints, &lclP->lambda));
    110:   PetscCall(VecDuplicate(tao->constraints, &lclP->WL));
    111:   PetscCall(VecDuplicate(tao->constraints, &lclP->lambda0));
    112:   PetscCall(VecDuplicate(tao->constraints, &lclP->con1));

    114:   PetscCall(VecSet(lclP->lambda, 0.0));

    116:   PetscCall(VecGetSize(tao->solution, &lclP->n));
    117:   PetscCall(VecGetSize(tao->constraints, &lclP->m));

    119:   PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->U));
    120:   PetscCall(VecCreate(((PetscObject)tao)->comm, &lclP->V));
    121:   PetscCall(ISGetLocalSize(tao->state_is, &nlocalstate));
    122:   PetscCall(ISGetLocalSize(tao->design_is, &nlocaldesign));
    123:   PetscCall(VecSetSizes(lclP->U, nlocalstate, lclP->m));
    124:   PetscCall(VecSetSizes(lclP->V, nlocaldesign, lclP->n - lclP->m));
    125:   PetscCall(VecSetType(lclP->U, ((PetscObject)tao->solution)->type_name));
    126:   PetscCall(VecSetType(lclP->V, ((PetscObject)tao->solution)->type_name));
    127:   PetscCall(VecSetFromOptions(lclP->U));
    128:   PetscCall(VecSetFromOptions(lclP->V));
    129:   PetscCall(VecDuplicate(lclP->U, &lclP->DU));
    130:   PetscCall(VecDuplicate(lclP->U, &lclP->U0));
    131:   PetscCall(VecDuplicate(lclP->U, &lclP->GU));
    132:   PetscCall(VecDuplicate(lclP->U, &lclP->GU0));
    133:   PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U));
    134:   PetscCall(VecDuplicate(lclP->U, &lclP->GL_U));
    135:   PetscCall(VecDuplicate(lclP->U, &lclP->GAugL_U0));
    136:   PetscCall(VecDuplicate(lclP->U, &lclP->GL_U0));
    137:   PetscCall(VecDuplicate(lclP->U, &lclP->WU));
    138:   PetscCall(VecDuplicate(lclP->U, &lclP->r));
    139:   PetscCall(VecDuplicate(lclP->V, &lclP->V0));
    140:   PetscCall(VecDuplicate(lclP->V, &lclP->V1));
    141:   PetscCall(VecDuplicate(lclP->V, &lclP->DV));
    142:   PetscCall(VecDuplicate(lclP->V, &lclP->s));
    143:   PetscCall(VecDuplicate(lclP->V, &lclP->GV));
    144:   PetscCall(VecDuplicate(lclP->V, &lclP->GV0));
    145:   PetscCall(VecDuplicate(lclP->V, &lclP->dbar));
    146:   PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V));
    147:   PetscCall(VecDuplicate(lclP->V, &lclP->GL_V));
    148:   PetscCall(VecDuplicate(lclP->V, &lclP->GAugL_V0));
    149:   PetscCall(VecDuplicate(lclP->V, &lclP->GL_V0));
    150:   PetscCall(VecDuplicate(lclP->V, &lclP->WV));
    151:   PetscCall(VecDuplicate(lclP->V, &lclP->g1));
    152:   PetscCall(VecDuplicate(lclP->V, &lclP->g2));

    154:   /* create scatters for state, design subvecs */
    155:   PetscCall(VecGetOwnershipRange(lclP->U, &lo, &hi));
    156:   PetscCall(ISCreateStride(((PetscObject)lclP->U)->comm, hi - lo, lo, 1, &is_state));
    157:   PetscCall(VecGetOwnershipRange(lclP->V, &lo, &hi));
    158:   if (0) {
    159:     PetscInt sizeU, sizeV;
    160:     PetscCall(VecGetSize(lclP->U, &sizeU));
    161:     PetscCall(VecGetSize(lclP->V, &sizeV));
    162:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "size(U)=%" PetscInt_FMT ", size(V)=%" PetscInt_FMT "\n", sizeU, sizeV));
    163:   }
    164:   PetscCall(ISCreateStride(((PetscObject)lclP->V)->comm, hi - lo, lo, 1, &is_design));
    165:   PetscCall(VecScatterCreate(tao->solution, tao->state_is, lclP->U, is_state, &lclP->state_scatter));
    166:   PetscCall(VecScatterCreate(tao->solution, tao->design_is, lclP->V, is_design, &lclP->design_scatter));
    167:   PetscCall(ISDestroy(&is_state));
    168:   PetscCall(ISDestroy(&is_design));
    169:   PetscFunctionReturn(PETSC_SUCCESS);
    170: }

    172: static PetscErrorCode TaoSolve_LCL(Tao tao)
    173: {
    174:   TAO_LCL                     *lclP = (TAO_LCL *)tao->data;
    175:   PetscInt                     phase2_iter, nlocal, its;
    176:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
    177:   PetscReal                    step      = 1.0, f, descent, aldescent;
    178:   PetscReal                    cnorm, mnorm;
    179:   PetscReal                    adec, r2, rGL_U, rWU;
    180:   PetscBool                    set, pset, flag, pflag, symmetric;

    182:   PetscFunctionBegin;
    183:   lclP->rho = lclP->rho0;
    184:   PetscCall(VecGetLocalSize(lclP->U, &nlocal));
    185:   PetscCall(VecGetLocalSize(lclP->V, &nlocal));
    186:   PetscCall(MatSetSizes(lclP->R, nlocal, nlocal, lclP->n - lclP->m, lclP->n - lclP->m));
    187:   PetscCall(MatLMVMAllocate(lclP->R, lclP->V, lclP->V));
    188:   lclP->recompute_jacobian_flag = PETSC_TRUE;

    190:   /* Scatter to U,V */
    191:   PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));

    193:   /* Evaluate Function, Gradient, Constraints, and Jacobian */
    194:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
    195:   PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
    196:   PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
    197:   PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));

    199:   /* Scatter gradient to GU,GV */
    200:   PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));

    202:   /* Evaluate Lagrangian function and gradient */
    203:   /* p0 */
    204:   PetscCall(VecSet(lclP->lambda, 0.0)); /*  Initial guess in CG */
    205:   PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
    206:   if (tao->jacobian_state_pre) {
    207:     PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
    208:   } else {
    209:     pset = pflag = PETSC_TRUE;
    210:   }
    211:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
    212:   else symmetric = PETSC_FALSE;

    214:   lclP->solve_type = LCL_ADJOINT2;
    215:   if (tao->jacobian_state_inv) {
    216:     if (symmetric) {
    217:       PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
    218:     } else {
    219:       PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
    220:     }
    221:   } else {
    222:     PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
    223:     if (symmetric) {
    224:       PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
    225:     } else {
    226:       PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
    227:     }
    228:     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    229:     tao->ksp_its += its;
    230:     tao->ksp_tot_its += its;
    231:   }
    232:   PetscCall(VecCopy(lclP->lambda, lclP->lambda0));
    233:   PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));

    235:   PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
    236:   PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));

    238:   /* Evaluate constraint norm */
    239:   PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
    240:   PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));

    242:   /* Monitor convergence */
    243:   tao->reason = TAO_CONTINUE_ITERATING;
    244:   PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
    245:   PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
    246:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

    248:   while (tao->reason == TAO_CONTINUE_ITERATING) {
    249:     /* Call general purpose update function */
    250:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
    251:     tao->ksp_its = 0;
    252:     /* Compute a descent direction for the linearly constrained subproblem
    253:        minimize f(u+du, v+dv)
    254:        s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */

    256:     /* Store the points around the linearization */
    257:     PetscCall(VecCopy(lclP->U, lclP->U0));
    258:     PetscCall(VecCopy(lclP->V, lclP->V0));
    259:     PetscCall(VecCopy(lclP->GU, lclP->GU0));
    260:     PetscCall(VecCopy(lclP->GV, lclP->GV0));
    261:     PetscCall(VecCopy(lclP->GAugL_U, lclP->GAugL_U0));
    262:     PetscCall(VecCopy(lclP->GAugL_V, lclP->GAugL_V0));
    263:     PetscCall(VecCopy(lclP->GL_U, lclP->GL_U0));
    264:     PetscCall(VecCopy(lclP->GL_V, lclP->GL_V0));

    266:     lclP->aug0 = lclP->aug;
    267:     lclP->lgn0 = lclP->lgn;

    269:     /* Given the design variables, we need to project the current iterate
    270:        onto the linearized constraint.  We choose to fix the design variables
    271:        and solve the linear system for the state variables.  The resulting
    272:        point is the Newton direction */

    274:     /* Solve r = A\con */
    275:     lclP->solve_type = LCL_FORWARD1;
    276:     PetscCall(VecSet(lclP->r, 0.0)); /*  Initial guess in CG */

    278:     if (tao->jacobian_state_inv) {
    279:       PetscCall(MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r));
    280:     } else {
    281:       PetscCall(KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre));
    282:       PetscCall(KSPSolve(tao->ksp, tao->constraints, lclP->r));
    283:       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    284:       tao->ksp_its += its;
    285:       tao->ksp_tot_its += tao->ksp_its;
    286:     }

    288:     /* Set design step direction dv to zero */
    289:     PetscCall(VecSet(lclP->s, 0.0));

    291:     /*
    292:        Check sufficient descent for constraint merit function .5*||con||^2
    293:        con' Ak r >= eps1 ||r||^(2+eps2)
    294:     */

    296:     /* Compute WU= Ak' * con */
    297:     if (symmetric) {
    298:       PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->WU));
    299:     } else {
    300:       PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->WU));
    301:     }
    302:     /* Compute r * Ak' * con */
    303:     PetscCall(VecDot(lclP->r, lclP->WU, &rWU));

    305:     /* compute ||r||^(2+eps2) */
    306:     PetscCall(VecNorm(lclP->r, NORM_2, &r2));
    307:     r2   = PetscPowScalar(r2, 2.0 + lclP->eps2);
    308:     adec = lclP->eps1 * r2;

    310:     if (rWU < adec) {
    311:       PetscCall(PetscInfo(tao, "Newton direction not descent for constraint, feasibility phase required\n"));
    312:       if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Newton direction not descent for constraint: %g -- using steepest descent\n", (double)descent));

    314:       PetscCall(PetscInfo(tao, "Using steepest descent direction instead.\n"));
    315:       PetscCall(VecSet(lclP->r, 0.0));
    316:       PetscCall(VecAXPY(lclP->r, -1.0, lclP->WU));
    317:       PetscCall(VecDot(lclP->r, lclP->r, &rWU));
    318:       PetscCall(VecNorm(lclP->r, NORM_2, &r2));
    319:       r2 = PetscPowScalar(r2, 2.0 + lclP->eps2);
    320:       PetscCall(VecDot(lclP->r, lclP->GAugL_U, &descent));
    321:       adec = lclP->eps1 * r2;
    322:     }

    324:     /*
    325:        Check descent for aug. lagrangian
    326:        r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2)
    327:           GL_U = GUk - Ak'*yk
    328:           WU   = Ak'*con
    329:                                          adec=eps1||r||^(2+eps2)

    331:        ==>
    332:        Check r'GL_U - rho*r'WU <= adec
    333:     */

    335:     PetscCall(VecDot(lclP->r, lclP->GL_U, &rGL_U));
    336:     aldescent = rGL_U - lclP->rho * rWU;
    337:     if (aldescent > -adec) {
    338:       if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Newton direction not descent for augmented Lagrangian: %g", (double)aldescent));
    339:       PetscCall(PetscInfo(tao, "Newton direction not descent for augmented Lagrangian: %g\n", (double)aldescent));
    340:       lclP->rho = (rGL_U - adec) / rWU;
    341:       PetscCheck(lclP->rho <= lclP->rhomax, PETSC_COMM_WORLD, PETSC_ERR_SUP, "rho=%g > rhomax, case not implemented.  Increase rhomax (-tao_lcl_rhomax)", (double)lclP->rho);
    342:       if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Increasing penalty parameter to %g\n", (double)lclP->rho));
    343:       PetscCall(PetscInfo(tao, "  Increasing penalty parameter to %g\n", (double)lclP->rho));
    344:     }

    346:     PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));
    347:     PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
    348:     PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));

    350:     /* We now minimize the augmented Lagrangian along the Newton direction */
    351:     PetscCall(VecScale(lclP->r, -1.0));
    352:     PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
    353:     PetscCall(VecScale(lclP->r, -1.0));
    354:     PetscCall(LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL));
    355:     PetscCall(LCLGather(lclP, lclP->U0, lclP->V0, lclP->X0));

    357:     lclP->recompute_jacobian_flag = PETSC_TRUE;

    359:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
    360:     PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
    361:     PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
    362:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
    363:     if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steplength = %g\n", (double)step));

    365:     PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
    366:     PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
    367:     PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));

    369:     PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));

    371:     /* Check convergence */
    372:     PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));
    373:     PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));
    374:     PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
    375:     PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
    376:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    377:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

    379:     /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */
    380:     for (phase2_iter = 0; phase2_iter < lclP->phase2_niter; phase2_iter++) {
    381:       /* We now minimize the objective function starting from the fraction of
    382:          the Newton point accepted by applying one step of a reduced-space
    383:          method.  The optimization problem is:

    385:          minimize f(u+du, v+dv)
    386:          s. t.    A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0)

    388:          In particular, we have that
    389:          du = -inv(A)*(Bdv + alpha g) */

    391:       PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
    392:       PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));

    394:       /* Store V and constraints */
    395:       PetscCall(VecCopy(lclP->V, lclP->V1));
    396:       PetscCall(VecCopy(tao->constraints, lclP->con1));

    398:       /* Compute multipliers */
    399:       /* p1 */
    400:       PetscCall(VecSet(lclP->lambda, 0.0)); /*  Initial guess in CG */
    401:       lclP->solve_type = LCL_ADJOINT1;
    402:       PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
    403:       if (tao->jacobian_state_pre) {
    404:         PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
    405:       } else {
    406:         pset = pflag = PETSC_TRUE;
    407:       }
    408:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
    409:       else symmetric = PETSC_FALSE;

    411:       if (tao->jacobian_state_inv) {
    412:         if (symmetric) {
    413:           PetscCall(MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
    414:         } else {
    415:           PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lambda));
    416:         }
    417:       } else {
    418:         if (symmetric) {
    419:           PetscCall(KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lambda));
    420:         } else {
    421:           PetscCall(KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lambda));
    422:         }
    423:         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    424:         tao->ksp_its += its;
    425:         tao->ksp_tot_its += its;
    426:       }
    427:       PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g1));
    428:       PetscCall(VecAXPY(lclP->g1, -1.0, lclP->GAugL_V));

    430:       /* Compute the limited-memory quasi-newton direction */
    431:       if (tao->niter > 0) {
    432:         PetscCall(MatSolve(lclP->R, lclP->g1, lclP->s));
    433:         PetscCall(VecDot(lclP->s, lclP->g1, &descent));
    434:         if (descent <= 0) {
    435:           if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space direction not descent: %g\n", (double)descent));
    436:           PetscCall(VecCopy(lclP->g1, lclP->s));
    437:         }
    438:       } else {
    439:         PetscCall(VecCopy(lclP->g1, lclP->s));
    440:       }
    441:       PetscCall(VecScale(lclP->g1, -1.0));

    443:       /* Recover the full space direction */
    444:       PetscCall(MatMult(tao->jacobian_design, lclP->s, lclP->WU));
    445:       /* PetscCall(VecSet(lclP->r,0.0)); */ /*  Initial guess in CG */
    446:       lclP->solve_type = LCL_FORWARD2;
    447:       if (tao->jacobian_state_inv) {
    448:         PetscCall(MatMult(tao->jacobian_state_inv, lclP->WU, lclP->r));
    449:       } else {
    450:         PetscCall(KSPSolve(tao->ksp, lclP->WU, lclP->r));
    451:         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    452:         tao->ksp_its += its;
    453:         tao->ksp_tot_its += its;
    454:       }

    456:       /* We now minimize the augmented Lagrangian along the direction -r,s */
    457:       PetscCall(VecScale(lclP->r, -1.0));
    458:       PetscCall(LCLGather(lclP, lclP->r, lclP->s, tao->stepdirection));
    459:       PetscCall(VecScale(lclP->r, -1.0));
    460:       lclP->recompute_jacobian_flag = PETSC_TRUE;

    462:       PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
    463:       PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHMT));
    464:       PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
    465:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason));
    466:       if (lclP->verbose) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reduced-space steplength =  %g\n", (double)step));

    468:       PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
    469:       PetscCall(LCLScatter(lclP, lclP->GL, lclP->GL_U, lclP->GL_V));
    470:       PetscCall(LCLScatter(lclP, lclP->GAugL, lclP->GAugL_U, lclP->GAugL_V));
    471:       PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
    472:       PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));

    474:       /* Compute the reduced gradient at the new point */

    476:       PetscCall(TaoComputeJacobianState(tao, lclP->X0, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
    477:       PetscCall(TaoComputeJacobianDesign(tao, lclP->X0, tao->jacobian_design));

    479:       /* p2 */
    480:       /* Compute multipliers, using lambda-rho*con as an initial guess in PCG */
    481:       if (phase2_iter == 0) {
    482:         PetscCall(VecWAXPY(lclP->lambda, -lclP->rho, lclP->con1, lclP->lambda0));
    483:       } else {
    484:         PetscCall(VecAXPY(lclP->lambda, -lclP->rho, tao->constraints));
    485:       }

    487:       PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
    488:       if (tao->jacobian_state_pre) {
    489:         PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
    490:       } else {
    491:         pset = pflag = PETSC_TRUE;
    492:       }
    493:       if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
    494:       else symmetric = PETSC_FALSE;

    496:       lclP->solve_type = LCL_ADJOINT2;
    497:       if (tao->jacobian_state_inv) {
    498:         if (symmetric) {
    499:           PetscCall(MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
    500:         } else {
    501:           PetscCall(MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lambda));
    502:         }
    503:       } else {
    504:         if (symmetric) {
    505:           PetscCall(KSPSolve(tao->ksp, lclP->GU, lclP->lambda));
    506:         } else {
    507:           PetscCall(KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lambda));
    508:         }
    509:         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    510:         tao->ksp_its += its;
    511:         tao->ksp_tot_its += its;
    512:       }

    514:       PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda, lclP->g2));
    515:       PetscCall(VecAXPY(lclP->g2, -1.0, lclP->GV));

    517:       PetscCall(VecScale(lclP->g2, -1.0));

    519:       /* Update the quasi-newton approximation */
    520:       PetscCall(MatLMVMUpdate(lclP->R, lclP->V, lclP->g2));
    521:       /* Use "-tao_ls_type gpcg -tao_ls_ftol 0 -tao_lmm_broyden_phi 0.0 -tao_lmm_scale_type scalar" to obtain agreement with MATLAB code */
    522:     }

    524:     PetscCall(VecCopy(lclP->lambda, lclP->lambda0));

    526:     /* Evaluate Function, Gradient, Constraints, and Jacobian */
    527:     PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
    528:     PetscCall(LCLScatter(lclP, tao->solution, lclP->U, lclP->V));
    529:     PetscCall(LCLScatter(lclP, tao->gradient, lclP->GU, lclP->GV));

    531:     PetscCall(TaoComputeJacobianState(tao, tao->solution, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
    532:     PetscCall(TaoComputeJacobianDesign(tao, tao->solution, tao->jacobian_design));
    533:     PetscCall(TaoComputeConstraints(tao, tao->solution, tao->constraints));

    535:     PetscCall(LCLComputeAugmentedLagrangianAndGradient(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao));

    537:     PetscCall(VecNorm(lclP->GAugL, NORM_2, &mnorm));

    539:     /* Evaluate constraint norm */
    540:     PetscCall(VecNorm(tao->constraints, NORM_2, &cnorm));

    542:     /* Monitor convergence */
    543:     tao->niter++;
    544:     PetscCall(TaoLogConvergenceHistory(tao, f, mnorm, cnorm, tao->ksp_its));
    545:     PetscCall(TaoMonitor(tao, tao->niter, f, mnorm, cnorm, step));
    546:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    547:   }
    548:   PetscFunctionReturn(PETSC_SUCCESS);
    549: }

    551: /*MC
    552:  TAOLCL - linearly constrained lagrangian method for pde-constrained optimization

    554: + -tao_lcl_eps1 - epsilon 1 tolerance
    555: . -tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,NULL);
    556: . -tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,NULL);
    557: . -tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,NULL);
    558: . -tao_lcl_phase2_niter - Number of phase 2 iterations in LCL algorithm
    559: . -tao_lcl_verbose - Print verbose output if True
    560: . -tao_lcl_tola - Tolerance for first forward solve
    561: . -tao_lcl_tolb - Tolerance for first adjoint solve
    562: . -tao_lcl_tolc - Tolerance for second forward solve
    563: - -tao_lcl_told - Tolerance for second adjoint solve

    565:   Level: beginner
    566: M*/
    567: PETSC_EXTERN PetscErrorCode TaoCreate_LCL(Tao tao)
    568: {
    569:   TAO_LCL    *lclP;
    570:   const char *morethuente_type = TAOLINESEARCHMT;

    572:   PetscFunctionBegin;
    573:   tao->ops->setup          = TaoSetup_LCL;
    574:   tao->ops->solve          = TaoSolve_LCL;
    575:   tao->ops->view           = TaoView_LCL;
    576:   tao->ops->setfromoptions = TaoSetFromOptions_LCL;
    577:   tao->ops->destroy        = TaoDestroy_LCL;
    578:   PetscCall(PetscNew(&lclP));
    579:   tao->data = (void *)lclP;

    581:   /* Override default settings (unless already changed) */
    582:   PetscCall(TaoParametersInitialize(tao));
    583:   PetscObjectParameterSetDefault(tao, max_it, 200);
    584:   PetscObjectParameterSetDefault(tao, catol, 1.0e-4);
    585:   PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
    586:   PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);
    587:   PetscObjectParameterSetDefault(tao, gttol, 1.0e-4);

    589:   lclP->rho0       = 1.0e-4;
    590:   lclP->rhomax     = 1e5;
    591:   lclP->eps1       = 1.0e-8;
    592:   lclP->eps2       = 0.0;
    593:   lclP->solve_type = 2;
    594:   lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4;
    595:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    596:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    597:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
    598:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));

    600:   PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, LCLComputeAugmentedLagrangianAndGradient, tao));
    601:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
    602:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
    603:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
    604:   PetscCall(KSPSetFromOptions(tao->ksp));

    606:   PetscCall(MatCreate(((PetscObject)tao)->comm, &lclP->R));
    607:   PetscCall(MatSetType(lclP->R, MATLMVMBFGS));
    608:   PetscFunctionReturn(PETSC_SUCCESS);
    609: }

    611: static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
    612: {
    613:   Tao       tao  = (Tao)ptr;
    614:   TAO_LCL  *lclP = (TAO_LCL *)tao->data;
    615:   PetscBool set, pset, flag, pflag, symmetric;
    616:   PetscReal cdotl;

    618:   PetscFunctionBegin;
    619:   PetscCall(TaoComputeObjectiveAndGradient(tao, X, f, G));
    620:   PetscCall(LCLScatter(lclP, G, lclP->GU, lclP->GV));
    621:   if (lclP->recompute_jacobian_flag) {
    622:     PetscCall(TaoComputeJacobianState(tao, X, tao->jacobian_state, tao->jacobian_state_pre, tao->jacobian_state_inv));
    623:     PetscCall(TaoComputeJacobianDesign(tao, X, tao->jacobian_design));
    624:   }
    625:   PetscCall(TaoComputeConstraints(tao, X, tao->constraints));
    626:   PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
    627:   if (tao->jacobian_state_pre) {
    628:     PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
    629:   } else {
    630:     pset = pflag = PETSC_TRUE;
    631:   }
    632:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
    633:   else symmetric = PETSC_FALSE;

    635:   PetscCall(VecDot(lclP->lambda0, tao->constraints, &cdotl));
    636:   lclP->lgn = *f - cdotl;

    638:   /* Gradient of Lagrangian GL = G - J' * lambda */
    639:   /*      WU = A' * WL
    640:           WV = B' * WL */
    641:   if (symmetric) {
    642:     PetscCall(MatMult(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
    643:   } else {
    644:     PetscCall(MatMultTranspose(tao->jacobian_state, lclP->lambda0, lclP->GL_U));
    645:   }
    646:   PetscCall(MatMultTranspose(tao->jacobian_design, lclP->lambda0, lclP->GL_V));
    647:   PetscCall(VecScale(lclP->GL_U, -1.0));
    648:   PetscCall(VecScale(lclP->GL_V, -1.0));
    649:   PetscCall(VecAXPY(lclP->GL_U, 1.0, lclP->GU));
    650:   PetscCall(VecAXPY(lclP->GL_V, 1.0, lclP->GV));
    651:   PetscCall(LCLGather(lclP, lclP->GL_U, lclP->GL_V, lclP->GL));

    653:   f[0] = lclP->lgn;
    654:   PetscCall(VecCopy(lclP->GL, G));
    655:   PetscFunctionReturn(PETSC_SUCCESS);
    656: }

    658: static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr)
    659: {
    660:   Tao       tao  = (Tao)ptr;
    661:   TAO_LCL  *lclP = (TAO_LCL *)tao->data;
    662:   PetscReal con2;
    663:   PetscBool flag, pflag, set, pset, symmetric;

    665:   PetscFunctionBegin;
    666:   PetscCall(LCLComputeLagrangianAndGradient(tao->linesearch, X, f, G, tao));
    667:   PetscCall(LCLScatter(lclP, G, lclP->GL_U, lclP->GL_V));
    668:   PetscCall(VecDot(tao->constraints, tao->constraints, &con2));
    669:   lclP->aug = lclP->lgn + 0.5 * lclP->rho * con2;

    671:   /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */
    672:   /*      WU = A' * c
    673:           WV = B' * c */
    674:   PetscCall(MatIsSymmetricKnown(tao->jacobian_state, &set, &flag));
    675:   if (tao->jacobian_state_pre) {
    676:     PetscCall(MatIsSymmetricKnown(tao->jacobian_state_pre, &pset, &pflag));
    677:   } else {
    678:     pset = pflag = PETSC_TRUE;
    679:   }
    680:   if (set && pset && flag && pflag) symmetric = PETSC_TRUE;
    681:   else symmetric = PETSC_FALSE;

    683:   if (symmetric) {
    684:     PetscCall(MatMult(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
    685:   } else {
    686:     PetscCall(MatMultTranspose(tao->jacobian_state, tao->constraints, lclP->GAugL_U));
    687:   }

    689:   PetscCall(MatMultTranspose(tao->jacobian_design, tao->constraints, lclP->GAugL_V));
    690:   PetscCall(VecAYPX(lclP->GAugL_U, lclP->rho, lclP->GL_U));
    691:   PetscCall(VecAYPX(lclP->GAugL_V, lclP->rho, lclP->GL_V));
    692:   PetscCall(LCLGather(lclP, lclP->GAugL_U, lclP->GAugL_V, lclP->GAugL));

    694:   f[0] = lclP->aug;
    695:   PetscCall(VecCopy(lclP->GAugL, G));
    696:   PetscFunctionReturn(PETSC_SUCCESS);
    697: }

    699: static PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x)
    700: {
    701:   PetscFunctionBegin;
    702:   PetscCall(VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
    703:   PetscCall(VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE));
    704:   PetscCall(VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
    705:   PetscCall(VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE));
    706:   PetscFunctionReturn(PETSC_SUCCESS);
    707: }
    708: static PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v)
    709: {
    710:   PetscFunctionBegin;
    711:   PetscCall(VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
    712:   PetscCall(VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD));
    713:   PetscCall(VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
    714:   PetscCall(VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD));
    715:   PetscFunctionReturn(PETSC_SUCCESS);
    716: }