1 #include "lcl.h" 2 #include "../src/tao/matrix/lmvmmat.h" 3 static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*); 4 static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch,Vec,PetscReal*,Vec,void*); 5 static PetscErrorCode LCLScatter(TAO_LCL*,Vec,Vec,Vec); 6 static PetscErrorCode LCLGather(TAO_LCL*,Vec,Vec,Vec); 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "TaoDestroy_LCL" 10 static PetscErrorCode TaoDestroy_LCL(TaoSolver tao) 11 { 12 TAO_LCL *lclP = (TAO_LCL*)tao->data; 13 PetscErrorCode ierr; 14 PetscFunctionBegin; 15 if (tao->setupcalled) { 16 ierr = MatDestroy(&lclP->R); CHKERRQ(ierr); 17 18 ierr = VecDestroy(&lclP->lamda); CHKERRQ(ierr); 19 ierr = VecDestroy(&lclP->lamda0); CHKERRQ(ierr); 20 ierr = VecDestroy(&lclP->WL); CHKERRQ(ierr); 21 ierr = VecDestroy(&lclP->W); CHKERRQ(ierr); 22 ierr = VecDestroy(&lclP->X0); CHKERRQ(ierr); 23 ierr = VecDestroy(&lclP->G0); CHKERRQ(ierr); 24 ierr = VecDestroy(&lclP->GL); CHKERRQ(ierr); 25 ierr = VecDestroy(&lclP->GAugL); CHKERRQ(ierr); 26 ierr = VecDestroy(&lclP->dbar); CHKERRQ(ierr); 27 ierr = VecDestroy(&lclP->U); CHKERRQ(ierr); 28 ierr = VecDestroy(&lclP->U0); CHKERRQ(ierr); 29 ierr = VecDestroy(&lclP->V); CHKERRQ(ierr); 30 ierr = VecDestroy(&lclP->V0); CHKERRQ(ierr); 31 ierr = VecDestroy(&lclP->V1); CHKERRQ(ierr); 32 ierr = VecDestroy(&lclP->GU); CHKERRQ(ierr); 33 ierr = VecDestroy(&lclP->GV); CHKERRQ(ierr); 34 ierr = VecDestroy(&lclP->GU0); CHKERRQ(ierr); 35 ierr = VecDestroy(&lclP->GV0); CHKERRQ(ierr); 36 ierr = VecDestroy(&lclP->GL_U); CHKERRQ(ierr); 37 ierr = VecDestroy(&lclP->GL_V); CHKERRQ(ierr); 38 ierr = VecDestroy(&lclP->GAugL_U); CHKERRQ(ierr); 39 ierr = VecDestroy(&lclP->GAugL_V); CHKERRQ(ierr); 40 ierr = VecDestroy(&lclP->GL_U0); CHKERRQ(ierr); 41 ierr = VecDestroy(&lclP->GL_V0); CHKERRQ(ierr); 42 ierr = VecDestroy(&lclP->GAugL_U0); CHKERRQ(ierr); 43 ierr = VecDestroy(&lclP->GAugL_V0); CHKERRQ(ierr); 44 ierr = VecDestroy(&lclP->DU); CHKERRQ(ierr); 45 ierr = VecDestroy(&lclP->DV); CHKERRQ(ierr); 46 ierr = VecDestroy(&lclP->WU); CHKERRQ(ierr); 47 ierr = VecDestroy(&lclP->WV); CHKERRQ(ierr); 48 ierr = VecDestroy(&lclP->g1); CHKERRQ(ierr); 49 ierr = VecDestroy(&lclP->g2); CHKERRQ(ierr); 50 ierr = VecDestroy(&lclP->con1); CHKERRQ(ierr); 51 52 ierr = VecDestroy(&lclP->r); CHKERRQ(ierr); 53 ierr = VecDestroy(&lclP->s); CHKERRQ(ierr); 54 55 ierr = ISDestroy(&tao->state_is); CHKERRQ(ierr); 56 ierr = ISDestroy(&tao->design_is); CHKERRQ(ierr); 57 58 ierr = VecScatterDestroy(&lclP->state_scatter); CHKERRQ(ierr); 59 ierr = VecScatterDestroy(&lclP->design_scatter); CHKERRQ(ierr); 60 } 61 ierr = PetscFree(tao->data); 62 tao->data = PETSC_NULL; 63 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "TaoSetFromOptions_LCL" 69 static PetscErrorCode TaoSetFromOptions_LCL(TaoSolver tao) 70 { 71 TAO_LCL *lclP = (TAO_LCL*)tao->data; 72 PetscBool flg; 73 PetscErrorCode ierr; 74 75 PetscFunctionBegin; 76 ierr = PetscOptionsHead("Linearly-Constrained Augmented Lagrangian Method for PDE-constrained optimization");CHKERRQ(ierr); 77 ierr = PetscOptionsReal("-tao_lcl_eps1","epsilon 1 tolerance","",lclP->eps1,&lclP->eps1,&flg); CHKERRQ(ierr); 78 ierr = PetscOptionsReal("-tao_lcl_eps2","epsilon 2 tolerance","",lclP->eps2,&lclP->eps2,&flg); CHKERRQ(ierr); 79 ierr = PetscOptionsReal("-tao_lcl_rho0","init value for rho","",lclP->rho0,&lclP->rho0,&flg); CHKERRQ(ierr); 80 ierr = PetscOptionsReal("-tao_lcl_rhomax","max value for rho","",lclP->rhomax,&lclP->rhomax,&flg); CHKERRQ(ierr); 81 lclP->phase2_niter = 1; 82 ierr = PetscOptionsInt("-tao_lcl_phase2_niter","Number of phase 2 iterations in LCL algorithm","",lclP->phase2_niter,&lclP->phase2_niter,&flg); CHKERRQ(ierr); 83 lclP->verbose = PETSC_FALSE; 84 ierr = PetscOptionsBool("-tao_lcl_verbose","Print verbose output","",lclP->verbose,&lclP->verbose,&flg); CHKERRQ(ierr); 85 lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4; 86 ierr = PetscOptionsReal("-tao_lcl_tola","Tolerance for first forward solve","",lclP->tau[0],&lclP->tau[0],&flg); CHKERRQ(ierr); 87 ierr = PetscOptionsReal("-tao_lcl_tolb","Tolerance for first adjoint solve","",lclP->tau[1],&lclP->tau[1],&flg); CHKERRQ(ierr); 88 ierr = PetscOptionsReal("-tao_lcl_tolc","Tolerance for second forward solve","",lclP->tau[2],&lclP->tau[2],&flg); CHKERRQ(ierr); 89 ierr = PetscOptionsReal("-tao_lcl_told","Tolerance for second adjoint solve","",lclP->tau[3],&lclP->tau[3],&flg); CHKERRQ(ierr); 90 91 ierr = PetscOptionsTail(); CHKERRQ(ierr); 92 ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr); 93 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "TaoView_LCL" 99 static PetscErrorCode TaoView_LCL(TaoSolver tao, PetscViewer viewer) 100 { 101 /* 102 TAO_LCL *lclP = (TAO_LCL*)tao->data; 103 PetscErrorCode ierr; 104 PetscFunctionBegin; 105 PetscFunctionReturn(0); 106 */ 107 return 0; 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "TaoSetup_LCL" 112 static PetscErrorCode TaoSetup_LCL(TaoSolver tao) 113 { 114 TAO_LCL *lclP = (TAO_LCL*)tao->data; 115 PetscInt lo, hi, nlocalstate, nlocaldesign; 116 PetscErrorCode ierr; 117 IS is_state, is_design; 118 PetscFunctionBegin; 119 /* Check for state/design IS */ 120 if (!tao->state_is) { 121 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONGSTATE,"LCL Solver requires an initial state index set -- use TaoSetStateIS()"); 122 } 123 ierr = VecDuplicate(tao->solution, &tao->gradient); CHKERRQ(ierr); 124 ierr = VecDuplicate(tao->solution, &tao->stepdirection); CHKERRQ(ierr); 125 ierr = VecDuplicate(tao->solution, &lclP->W); CHKERRQ(ierr); 126 ierr = VecDuplicate(tao->solution, &lclP->X0); CHKERRQ(ierr); 127 ierr = VecDuplicate(tao->solution, &lclP->G0); CHKERRQ(ierr); 128 ierr = VecDuplicate(tao->solution, &lclP->GL); CHKERRQ(ierr); 129 ierr = VecDuplicate(tao->solution, &lclP->GAugL); CHKERRQ(ierr); 130 131 ierr = VecDuplicate(tao->constraints, &lclP->lamda); CHKERRQ(ierr); 132 ierr = VecDuplicate(tao->constraints, &lclP->WL); CHKERRQ(ierr); 133 ierr = VecDuplicate(tao->constraints, &lclP->lamda0); CHKERRQ(ierr); 134 ierr = VecDuplicate(tao->constraints, &lclP->con1); CHKERRQ(ierr); 135 136 ierr = VecSet(lclP->lamda,0.0); CHKERRQ(ierr); 137 138 ierr = VecGetSize(tao->solution, &lclP->n); CHKERRQ(ierr); 139 ierr = VecGetSize(tao->constraints, &lclP->m); CHKERRQ(ierr); 140 141 142 ierr = VecCreate(((PetscObject)tao)->comm,&lclP->U); CHKERRQ(ierr); 143 ierr = VecCreate(((PetscObject)tao)->comm,&lclP->V); CHKERRQ(ierr); 144 ierr = ISGetLocalSize(tao->state_is,&nlocalstate); CHKERRQ(ierr); 145 ierr = ISGetLocalSize(tao->design_is,&nlocaldesign); CHKERRQ(ierr); 146 ierr = VecSetSizes(lclP->U,nlocalstate,lclP->m); CHKERRQ(ierr); 147 ierr = VecSetSizes(lclP->V,nlocaldesign,lclP->n-lclP->m); CHKERRQ(ierr); 148 ierr = VecSetType(lclP->U,((PetscObject)(tao->solution))->type_name); CHKERRQ(ierr); 149 ierr = VecSetType(lclP->V,((PetscObject)(tao->solution))->type_name); CHKERRQ(ierr); 150 ierr = VecSetFromOptions(lclP->U); CHKERRQ(ierr); 151 ierr = VecSetFromOptions(lclP->V); CHKERRQ(ierr); 152 ierr = VecDuplicate(lclP->U,&lclP->DU); CHKERRQ(ierr); 153 ierr = VecDuplicate(lclP->U,&lclP->U0); CHKERRQ(ierr); 154 ierr = VecDuplicate(lclP->U,&lclP->GU); CHKERRQ(ierr); 155 ierr = VecDuplicate(lclP->U,&lclP->GU0); CHKERRQ(ierr); 156 ierr = VecDuplicate(lclP->U,&lclP->GAugL_U); CHKERRQ(ierr); 157 ierr = VecDuplicate(lclP->U,&lclP->GL_U); CHKERRQ(ierr); 158 ierr = VecDuplicate(lclP->U,&lclP->GAugL_U0); CHKERRQ(ierr); 159 ierr = VecDuplicate(lclP->U,&lclP->GL_U0); CHKERRQ(ierr); 160 ierr = VecDuplicate(lclP->U,&lclP->WU); CHKERRQ(ierr); 161 ierr = VecDuplicate(lclP->U,&lclP->r); CHKERRQ(ierr); 162 ierr = VecDuplicate(lclP->V,&lclP->V0); CHKERRQ(ierr); 163 ierr = VecDuplicate(lclP->V,&lclP->V1); CHKERRQ(ierr); 164 ierr = VecDuplicate(lclP->V,&lclP->DV); CHKERRQ(ierr); 165 ierr = VecDuplicate(lclP->V,&lclP->s); CHKERRQ(ierr); 166 ierr = VecDuplicate(lclP->V,&lclP->GV); CHKERRQ(ierr); 167 ierr = VecDuplicate(lclP->V,&lclP->GV0); CHKERRQ(ierr); 168 ierr = VecDuplicate(lclP->V,&lclP->dbar); CHKERRQ(ierr); 169 ierr = VecDuplicate(lclP->V,&lclP->GAugL_V); CHKERRQ(ierr); 170 ierr = VecDuplicate(lclP->V,&lclP->GL_V); CHKERRQ(ierr); 171 ierr = VecDuplicate(lclP->V,&lclP->GAugL_V0); CHKERRQ(ierr); 172 ierr = VecDuplicate(lclP->V,&lclP->GL_V0); CHKERRQ(ierr); 173 ierr = VecDuplicate(lclP->V,&lclP->WV); CHKERRQ(ierr); 174 ierr = VecDuplicate(lclP->V,&lclP->g1); CHKERRQ(ierr); 175 ierr = VecDuplicate(lclP->V,&lclP->g2); CHKERRQ(ierr); 176 177 178 179 180 /* create scatters for state, design subvecs */ 181 ierr = VecGetOwnershipRange(lclP->U,&lo,&hi); CHKERRQ(ierr); 182 ierr = ISCreateStride(((PetscObject)lclP->U)->comm,hi-lo,lo,1,&is_state); CHKERRQ(ierr); 183 ierr = VecGetOwnershipRange(lclP->V,&lo,&hi); CHKERRQ(ierr); 184 if (0) { 185 PetscInt sizeU,sizeV; 186 ierr = VecGetSize(lclP->U,&sizeU); 187 ierr = VecGetSize(lclP->V,&sizeV); 188 ierr = PetscPrintf(PETSC_COMM_WORLD,"size(U)=%D, size(V)=%D\n",sizeU,sizeV); 189 } 190 ierr = ISCreateStride(((PetscObject)lclP->V)->comm,hi-lo,lo,1,&is_design); CHKERRQ(ierr); 191 ierr = VecScatterCreate(tao->solution,tao->state_is,lclP->U,is_state,&lclP->state_scatter); CHKERRQ(ierr); 192 ierr = VecScatterCreate(tao->solution,tao->design_is,lclP->V,is_design,&lclP->design_scatter); CHKERRQ(ierr); 193 ierr = ISDestroy(&is_state); CHKERRQ(ierr); 194 ierr = ISDestroy(&is_design); CHKERRQ(ierr); 195 196 197 198 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "TaoSolve_LCL" 204 static PetscErrorCode TaoSolve_LCL(TaoSolver tao) 205 { 206 TAO_LCL *lclP = (TAO_LCL*)tao->data; 207 PetscInt iter=0,phase2_iter,nlocal,its; 208 TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING; 209 TaoLineSearchTerminationReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 210 PetscReal step=1.0,f, descent, aldescent; 211 PetscReal cnorm, mnorm; 212 PetscReal adec,r2,rGL_U,rWU; 213 PetscBool set,pset,flag,pflag,symmetric; 214 PetscErrorCode ierr; 215 PetscFunctionBegin; 216 217 lclP->rho = lclP->rho0; 218 ierr = VecGetLocalSize(lclP->U,&nlocal); CHKERRQ(ierr); 219 ierr = VecGetLocalSize(lclP->V,&nlocal); CHKERRQ(ierr); 220 ierr = MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lclP->n-lclP->m,&lclP->R); CHKERRQ(ierr); 221 ierr = MatLMVMAllocateVectors(lclP->R,lclP->V); CHKERRQ(ierr); 222 lclP->recompute_jacobian_flag = PETSC_TRUE; 223 224 /* Scatter to U,V */ 225 ierr = LCLScatter(lclP,tao->solution,lclP->U,lclP->V); CHKERRQ(ierr); 226 227 /* Evaluate Function, Gradient, Constraints, and Jacobian */ 228 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient); CHKERRQ(ierr); 229 ierr = TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag); CHKERRQ(ierr); 230 ierr = TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design); CHKERRQ(ierr); 231 ierr = TaoComputeConstraints(tao,tao->solution, tao->constraints); CHKERRQ(ierr); 232 233 234 /* Scatter gradient to GU,GV */ 235 ierr = LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV); CHKERRQ(ierr); 236 237 238 /* Evaluate Lagrangian function and gradient */ 239 /* p0 */ 240 ierr = VecSet(lclP->lamda,0.0); CHKERRQ(ierr); /* Initial guess in CG */ 241 ierr = MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); CHKERRQ(ierr); 242 if (tao->jacobian_state_pre) { 243 ierr = MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); CHKERRQ(ierr); 244 } else { 245 pset = pflag = PETSC_TRUE; 246 } 247 if (set && pset && flag && pflag) 248 symmetric = PETSC_TRUE; 249 else 250 symmetric = PETSC_FALSE; 251 252 lclP->solve_type = LCL_ADJOINT2; 253 if (tao->jacobian_state_inv) { 254 if (symmetric) { 255 ierr = MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); CHKERRQ(ierr); } else { 256 ierr = MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda); CHKERRQ(ierr); 257 } 258 } else { 259 ierr = KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre, lclP->statematflag); CHKERRQ(ierr); 260 if (symmetric) { 261 ierr = KSPSolve(tao->ksp, lclP->GU, lclP->lamda); CHKERRQ(ierr); 262 } else { 263 ierr = KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda); CHKERRQ(ierr); 264 } 265 ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr); 266 tao->ksp_its += its; 267 } 268 269 ierr = VecCopy(lclP->lamda,lclP->lamda0); CHKERRQ(ierr); 270 271 ierr = LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao); CHKERRQ(ierr); 272 273 ierr = LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V); CHKERRQ(ierr); 274 ierr = LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V); CHKERRQ(ierr); 275 276 277 /* Evaluate constraint norm */ 278 ierr = VecNorm(tao->constraints, NORM_2, &cnorm); CHKERRQ(ierr); 279 ierr = VecNorm(lclP->GAugL, NORM_2, &mnorm); CHKERRQ(ierr); 280 281 282 /* Monitor convergence */ 283 ierr = TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason); CHKERRQ(ierr); 284 285 while (reason == TAO_CONTINUE_ITERATING) { 286 /* Compute a descent direction for the linearly constrained subproblem 287 minimize f(u+du, v+dv) 288 s.t. A(u0,v0)du + B(u0,v0)dv = -g(u0,v0) */ 289 290 /* Store the points around the linearization */ 291 ierr = VecCopy(lclP->U, lclP->U0); CHKERRQ(ierr); 292 ierr = VecCopy(lclP->V, lclP->V0); CHKERRQ(ierr); 293 ierr = VecCopy(lclP->GU,lclP->GU0); CHKERRQ(ierr); 294 ierr = VecCopy(lclP->GV,lclP->GV0); CHKERRQ(ierr); 295 ierr = VecCopy(lclP->GAugL_U,lclP->GAugL_U0); CHKERRQ(ierr); 296 ierr = VecCopy(lclP->GAugL_V,lclP->GAugL_V0); CHKERRQ(ierr); 297 ierr = VecCopy(lclP->GL_U,lclP->GL_U0); CHKERRQ(ierr); 298 ierr = VecCopy(lclP->GL_V,lclP->GL_V0); CHKERRQ(ierr); 299 300 lclP->aug0 = lclP->aug; 301 lclP->lgn0 = lclP->lgn; 302 303 /* Given the design variables, we need to project the current iterate 304 onto the linearized constraint. We choose to fix the design variables 305 and solve the linear system for the state variables. The resulting 306 point is the Newton direction */ 307 308 /* Solve r = A\con */ 309 lclP->solve_type = LCL_FORWARD1; 310 ierr = VecSet(lclP->r,0.0); CHKERRQ(ierr); /* Initial guess in CG */ 311 312 if (tao->jacobian_state_inv) { 313 ierr = MatMult(tao->jacobian_state_inv, tao->constraints, lclP->r); CHKERRQ(ierr); 314 } else { 315 ierr = KSPSetOperators(tao->ksp, tao->jacobian_state, tao->jacobian_state_pre, lclP->statematflag); CHKERRQ(ierr); 316 ierr = KSPSolve(tao->ksp, tao->constraints, lclP->r); CHKERRQ(ierr); 317 ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr); 318 tao->ksp_its += its; 319 } 320 321 /* Set design step direction dv to zero */ 322 ierr = VecSet(lclP->s, 0.0); CHKERRQ(ierr); 323 324 /* 325 Check sufficient descent for constraint merit function .5*||con||^2 326 con' Ak r >= eps1 ||r||^(2+eps2) 327 */ 328 329 /* Compute WU= Ak' * con */ 330 if (symmetric) { 331 ierr = MatMult(tao->jacobian_state,tao->constraints,lclP->WU); CHKERRQ(ierr); 332 } else { 333 ierr = MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->WU); CHKERRQ(ierr); 334 } 335 /* Compute r * Ak' * con */ 336 ierr = VecDot(lclP->r,lclP->WU,&rWU); CHKERRQ(ierr); 337 338 /* compute ||r||^(2+eps2) */ 339 ierr = VecNorm(lclP->r,NORM_2,&r2); CHKERRQ(ierr); 340 r2 = PetscPowScalar(r2,2.0+lclP->eps2); 341 adec = lclP->eps1 * r2; 342 343 if (rWU < adec) { 344 ierr = PetscInfo(tao,"Newton direction not descent for constraint, feasibility phase required"); CHKERRQ(ierr); 345 if (lclP->verbose) { 346 ierr = PetscPrintf(PETSC_COMM_WORLD,"Newton direction not descent for constraint: %G -- using steepest descent\n",descent); CHKERRQ(ierr); 347 } 348 349 ierr = PetscInfo(tao,"Using steepest descent direction instead.\n"); CHKERRQ(ierr); 350 ierr = VecSet(lclP->r,0.0); CHKERRQ(ierr); 351 ierr = VecAXPY(lclP->r,-1.0,lclP->WU); CHKERRQ(ierr); 352 ierr = VecDot(lclP->r,lclP->r,&rWU); CHKERRQ(ierr); 353 ierr = VecNorm(lclP->r,NORM_2,&r2); CHKERRQ(ierr); 354 r2 = PetscPowScalar(r2,2.0+lclP->eps2); 355 ierr = VecDot(lclP->r,lclP->GAugL_U,&descent); CHKERRQ(ierr); 356 adec = lclP->eps1 * r2; 357 } 358 359 360 /* 361 Check descent for aug. lagrangian 362 r' (GUk - Ak'*yk - rho*Ak'*con) <= -eps1 ||r||^(2+eps2) 363 GL_U = GUk - Ak'*yk 364 WU = Ak'*con 365 adec=eps1||r||^(2+eps2) 366 367 ==> 368 Check r'GL_U - rho*r'WU <= adec 369 */ 370 371 ierr = VecDot(lclP->r,lclP->GL_U,&rGL_U); 372 aldescent = rGL_U - lclP->rho*rWU; 373 if (aldescent > -adec) { 374 if (lclP->verbose) { 375 ierr = PetscPrintf(PETSC_COMM_WORLD," Newton direction not descent for augmented Lagrangian: %G",aldescent); CHKERRQ(ierr); 376 } 377 ierr = PetscInfo1(tao,"Newton direction not descent for augmented Lagrangian: %G",aldescent); CHKERRQ(ierr); 378 lclP->rho = (rGL_U - adec)/rWU; 379 if (lclP->rho > lclP->rhomax) { 380 lclP->rho = lclP->rhomax; 381 SETERRQ1(PETSC_COMM_WORLD,0,"rho=%G > rhomax, case not implemented. Increase rhomax (-tao_lcl_rhomax)",lclP->rho); 382 } 383 if (lclP->verbose) { 384 ierr = PetscPrintf(PETSC_COMM_WORLD," Increasing penalty parameter to %G\n",lclP->rho); CHKERRQ(ierr); 385 } 386 ierr = PetscInfo1(tao," Increasing penalty parameter to %G",lclP->rho); 387 } 388 389 ierr = LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao); CHKERRQ(ierr); 390 ierr = LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V); CHKERRQ(ierr); 391 ierr = LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V); CHKERRQ(ierr); 392 393 /* We now minimize the augmented Lagrangian along the Newton direction */ 394 ierr = VecScale(lclP->r,-1.0); CHKERRQ(ierr); 395 ierr = LCLGather(lclP, lclP->r,lclP->s,tao->stepdirection); CHKERRQ(ierr); 396 ierr = VecScale(lclP->r,-1.0); CHKERRQ(ierr); 397 ierr = LCLGather(lclP, lclP->GAugL_U0, lclP->GAugL_V0, lclP->GAugL); CHKERRQ(ierr); 398 ierr = LCLGather(lclP, lclP->U0,lclP->V0,lclP->X0); CHKERRQ(ierr); 399 400 lclP->recompute_jacobian_flag = PETSC_TRUE; 401 402 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); CHKERRQ(ierr); 403 ierr = TaoLineSearchSetType(tao->linesearch, TAOLINESEARCH_MT); CHKERRQ(ierr); 404 ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr); 405 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection, &step, &ls_reason); CHKERRQ(ierr); 406 if (lclP->verbose) { 407 ierr = PetscPrintf(PETSC_COMM_WORLD,"Steplength = %G\n",step); CHKERRQ(ierr); 408 } 409 410 ierr = LCLScatter(lclP,tao->solution,lclP->U,lclP->V); CHKERRQ(ierr); 411 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient); CHKERRQ(ierr); 412 ierr = LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV); CHKERRQ(ierr); 413 414 ierr = LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V); CHKERRQ(ierr); 415 416 /* Check convergence */ 417 ierr = VecNorm(lclP->GAugL, NORM_2, &mnorm); CHKERRQ(ierr); 418 ierr = VecNorm(tao->constraints, NORM_2, &cnorm); CHKERRQ(ierr); 419 ierr = TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason); CHKERRQ(ierr); 420 if (reason != TAO_CONTINUE_ITERATING){ 421 break; 422 } 423 424 425 /* TODO: use a heuristic to choose how many iterations should be performed within phase 2 */ 426 for (phase2_iter=0; phase2_iter<lclP->phase2_niter; phase2_iter++){ 427 /* We now minimize the objective function starting from the fraction of 428 the Newton point accepted by applying one step of a reduced-space 429 method. The optimization problem is: 430 431 minimize f(u+du, v+dv) 432 s. t. A(u0,v0)du + B(u0,v0)du = -alpha g(u0,v0) 433 434 In particular, we have that 435 du = -inv(A)*(Bdv + alpha g) */ 436 437 ierr = TaoComputeJacobianState(tao,lclP->X0,&tao->jacobian_state,&tao->jacobian_state_pre,&tao->jacobian_state_inv,&lclP->statematflag); CHKERRQ(ierr); 438 ierr = TaoComputeJacobianDesign(tao,lclP->X0,&tao->jacobian_design); CHKERRQ(ierr); 439 440 /* Store V and constraints */ 441 ierr = VecCopy(lclP->V, lclP->V1); CHKERRQ(ierr); 442 ierr = VecCopy(tao->constraints,lclP->con1); CHKERRQ(ierr); 443 444 /* Compute multipliers */ 445 /* p1 */ 446 ierr = VecSet(lclP->lamda,0.0); CHKERRQ(ierr); /* Initial guess in CG */ 447 lclP->solve_type = LCL_ADJOINT1; 448 ierr = MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); CHKERRQ(ierr); 449 if (tao->jacobian_state_pre) { 450 ierr = MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); CHKERRQ(ierr); 451 } else { 452 pset = pflag = PETSC_TRUE; 453 } 454 if (set && pset && flag && pflag) 455 symmetric = PETSC_TRUE; 456 else 457 symmetric = PETSC_FALSE; 458 459 if (tao->jacobian_state_inv) { 460 if (symmetric) { 461 ierr = MatMult(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda); CHKERRQ(ierr); 462 } else { 463 ierr = MatMultTranspose(tao->jacobian_state_inv, lclP->GAugL_U, lclP->lamda); CHKERRQ(ierr); 464 } 465 } else { 466 if (symmetric) { 467 ierr = KSPSolve(tao->ksp, lclP->GAugL_U, lclP->lamda); CHKERRQ(ierr); 468 } else { 469 ierr = KSPSolveTranspose(tao->ksp, lclP->GAugL_U, lclP->lamda); CHKERRQ(ierr); 470 } 471 ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr); 472 tao->ksp_its += its; 473 } 474 ierr = MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g1); CHKERRQ(ierr); 475 ierr = VecAXPY(lclP->g1,-1.0,lclP->GAugL_V); CHKERRQ(ierr); 476 477 /* Compute the limited-memory quasi-newton direction */ 478 if (iter > 0) { 479 ierr = MatLMVMSolve(lclP->R,lclP->g1,lclP->s); CHKERRQ(ierr); 480 ierr = VecDot(lclP->s,lclP->g1,&descent); CHKERRQ(ierr); 481 if (descent < 0e-8) { 482 if (lclP->verbose) { 483 ierr = PetscPrintf(PETSC_COMM_WORLD,"Reduced-space direction not descent: %G\n",descent); CHKERRQ(ierr); 484 } 485 ierr = VecCopy(lclP->g1,lclP->s); CHKERRQ(ierr); 486 } 487 } else { 488 ierr = VecCopy(lclP->g1,lclP->s); CHKERRQ(ierr); 489 } 490 ierr = VecScale(lclP->g1,-1.0); CHKERRQ(ierr); 491 492 493 /* Recover the full space direction */ 494 ierr = MatMult(tao->jacobian_design,lclP->s,lclP->WU); CHKERRQ(ierr); 495 //ierr = VecSet(lclP->r,0.0); CHKERRQ(ierr); /* Initial guess in CG */ 496 lclP->solve_type = LCL_FORWARD2; 497 if (tao->jacobian_state_inv) { 498 ierr = MatMult(tao->jacobian_state_inv,lclP->WU,lclP->r); CHKERRQ(ierr); 499 } else { 500 ierr = KSPSolve(tao->ksp, lclP->WU, lclP->r); CHKERRQ(ierr); 501 ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr); 502 tao->ksp_its += its; 503 } 504 505 /* We now minimize the augmented Lagrangian along the direction -r,s */ 506 ierr = VecScale(lclP->r, -1.0); CHKERRQ(ierr); 507 ierr = LCLGather(lclP,lclP->r,lclP->s,tao->stepdirection); CHKERRQ(ierr); 508 ierr = VecScale(lclP->r, -1.0); CHKERRQ(ierr); 509 lclP->recompute_jacobian_flag = PETSC_TRUE; 510 511 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); CHKERRQ(ierr); 512 ierr = TaoLineSearchSetType(tao->linesearch,TAOLINESEARCH_MT); CHKERRQ(ierr); 513 ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr); 514 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &lclP->aug, lclP->GAugL, tao->stepdirection,&step,&ls_reason); CHKERRQ(ierr); 515 if (lclP->verbose){ 516 ierr = PetscPrintf(PETSC_COMM_WORLD,"Reduced-space steplength = %G\n",step); CHKERRQ(ierr); 517 } 518 519 ierr = LCLScatter(lclP,tao->solution,lclP->U,lclP->V); CHKERRQ(ierr); 520 ierr = LCLScatter(lclP,lclP->GL,lclP->GL_U,lclP->GL_V); CHKERRQ(ierr); 521 ierr = LCLScatter(lclP,lclP->GAugL,lclP->GAugL_U,lclP->GAugL_V); CHKERRQ(ierr); 522 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient); CHKERRQ(ierr); 523 ierr = LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV); CHKERRQ(ierr); 524 525 /* Compute the reduced gradient at the new point */ 526 527 ierr = TaoComputeJacobianState(tao,lclP->X0,&tao->jacobian_state,&tao->jacobian_state_pre,&tao->jacobian_state_inv,&lclP->statematflag); CHKERRQ(ierr); 528 ierr = TaoComputeJacobianDesign(tao,lclP->X0,&tao->jacobian_design); CHKERRQ(ierr); 529 530 /* p2 */ 531 /* Compute multipliers, using lamda-rho*con as an initial guess in PCG */ 532 if (phase2_iter==0){ 533 ierr = VecWAXPY(lclP->lamda,-lclP->rho,lclP->con1,lclP->lamda0); CHKERRQ(ierr); 534 } 535 else { 536 ierr = VecAXPY(lclP->lamda,-lclP->rho,tao->constraints); CHKERRQ(ierr); 537 } 538 539 ierr = MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); CHKERRQ(ierr); 540 if (tao->jacobian_state_pre) { 541 ierr = MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); CHKERRQ(ierr); 542 } else { 543 pset = pflag = PETSC_TRUE; 544 } 545 if (set && pset && flag && pflag) 546 symmetric = PETSC_TRUE; 547 else 548 symmetric = PETSC_FALSE; 549 550 lclP->solve_type = LCL_ADJOINT2; 551 if (tao->jacobian_state_inv) { 552 if (symmetric) { 553 ierr = MatMult(tao->jacobian_state_inv, lclP->GU, lclP->lamda); CHKERRQ(ierr); 554 } else { 555 ierr = MatMultTranspose(tao->jacobian_state_inv, lclP->GU, lclP->lamda); CHKERRQ(ierr); 556 } 557 } else { 558 if (symmetric) { 559 ierr = KSPSolve(tao->ksp, lclP->GU, lclP->lamda); CHKERRQ(ierr); 560 } else { 561 ierr = KSPSolveTranspose(tao->ksp, lclP->GU, lclP->lamda); CHKERRQ(ierr); 562 } 563 ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr); 564 tao->ksp_its += its; 565 } 566 567 568 ierr = MatMultTranspose(tao->jacobian_design,lclP->lamda,lclP->g2); CHKERRQ(ierr); 569 ierr = VecAXPY(lclP->g2,-1.0,lclP->GV); CHKERRQ(ierr); 570 571 ierr = VecScale(lclP->g2,-1.0); CHKERRQ(ierr); 572 573 /* Update the quasi-newton approximation */ 574 if (phase2_iter >= 0){ 575 ierr = MatLMVMSetPrev(lclP->R,lclP->V1,lclP->g1); 576 } 577 ierr = MatLMVMUpdate(lclP->R,lclP->V,lclP->g2); CHKERRQ(ierr); 578 /* 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 */ 579 580 } 581 582 ierr = VecCopy(lclP->lamda,lclP->lamda0); CHKERRQ(ierr); 583 584 /* Evaluate Function, Gradient, Constraints, and Jacobian */ 585 ierr = TaoComputeObjectiveAndGradient(tao,tao->solution,&f,tao->gradient); CHKERRQ(ierr); 586 ierr = LCLScatter(lclP,tao->solution,lclP->U,lclP->V); CHKERRQ(ierr); 587 ierr = LCLScatter(lclP,tao->gradient,lclP->GU,lclP->GV); CHKERRQ(ierr); 588 589 ierr = TaoComputeJacobianState(tao,tao->solution, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag); CHKERRQ(ierr); 590 ierr = TaoComputeJacobianDesign(tao,tao->solution, &tao->jacobian_design); CHKERRQ(ierr); 591 ierr = TaoComputeConstraints(tao,tao->solution, tao->constraints); CHKERRQ(ierr); 592 593 594 ierr = LCLComputeAugmentedLagrangianAndGradient(tao->linesearch,tao->solution,&lclP->aug,lclP->GAugL,tao); CHKERRQ(ierr); 595 596 ierr = VecNorm(lclP->GAugL, NORM_2, &mnorm); CHKERRQ(ierr); 597 598 /* Evaluate constraint norm */ 599 ierr = VecNorm(tao->constraints, NORM_2, &cnorm); CHKERRQ(ierr); 600 601 /* Monitor convergence */ 602 iter++; 603 ierr = TaoMonitor(tao, iter,f,mnorm,cnorm,step,&reason); CHKERRQ(ierr); 604 605 } 606 607 ierr = MatDestroy(&lclP->R); CHKERRQ(ierr); 608 609 PetscFunctionReturn(0); 610 } 611 612 EXTERN_C_BEGIN 613 #undef __FUNCT__ 614 #define __FUNCT__ "TaoCreate_LCL" 615 PetscErrorCode TaoCreate_LCL(TaoSolver tao) 616 { 617 TAO_LCL *lclP; 618 PetscErrorCode ierr; 619 const char *morethuente_type = TAOLINESEARCH_MT; 620 PetscFunctionBegin; 621 622 tao->ops->setup = TaoSetup_LCL; 623 tao->ops->solve = TaoSolve_LCL; 624 tao->ops->view = TaoView_LCL; 625 tao->ops->setfromoptions = TaoSetFromOptions_LCL; 626 tao->ops->destroy = TaoDestroy_LCL; 627 628 629 ierr = PetscNewLog(tao,TAO_LCL,&lclP); CHKERRQ(ierr); 630 tao->data = (void*)lclP; 631 632 633 tao->max_it=200; 634 tao->fatol=1e-8; 635 tao->frtol=1e-8; 636 tao->catol=1e-4; 637 tao->crtol=1e-4; 638 tao->gttol=1e-4; 639 tao->gatol=1e-4; 640 tao->grtol=1e-4; 641 lclP->rho0 = 1.0e-4; 642 lclP->rhomax=1e5; 643 lclP->eps1 = 1.0e-8; 644 lclP->eps2 = 0.0; 645 lclP->solve_type=2; 646 lclP->tau[0] = lclP->tau[1] = lclP->tau[2] = lclP->tau[3] = 1.0e-4; 647 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); CHKERRQ(ierr); 648 ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type); CHKERRQ(ierr); 649 650 ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LCLComputeAugmentedLagrangianAndGradient, tao); CHKERRQ(ierr); 651 ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp); CHKERRQ(ierr); 652 ierr = KSPSetFromOptions(tao->ksp); CHKERRQ(ierr); 653 654 655 656 PetscFunctionReturn(0); 657 658 } 659 EXTERN_C_END 660 661 662 #undef __FUNCT__ 663 #define __FUNCT__ "LCLComputeLagrangianAndGradient" 664 static PetscErrorCode LCLComputeLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr) 665 { 666 TaoSolver tao = (TaoSolver)ptr; 667 TAO_LCL *lclP = (TAO_LCL*)tao->data; 668 PetscBool set,pset,flag,pflag,symmetric; 669 PetscReal cdotl; 670 PetscErrorCode ierr; 671 672 PetscFunctionBegin; 673 ierr = TaoComputeObjectiveAndGradient(tao,X,f,G); CHKERRQ(ierr); 674 ierr = LCLScatter(lclP,G,lclP->GU,lclP->GV); CHKERRQ(ierr); 675 if (lclP->recompute_jacobian_flag) { 676 ierr = TaoComputeJacobianState(tao,X, &tao->jacobian_state, &tao->jacobian_state_pre, &tao->jacobian_state_inv, &lclP->statematflag); CHKERRQ(ierr); 677 ierr = TaoComputeJacobianDesign(tao,X, &tao->jacobian_design); CHKERRQ(ierr); 678 } 679 ierr = TaoComputeConstraints(tao,X, tao->constraints); CHKERRQ(ierr); 680 ierr = MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); CHKERRQ(ierr); 681 if (tao->jacobian_state_pre) { 682 ierr = MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); CHKERRQ(ierr); 683 } else { 684 pset = pflag = PETSC_TRUE; 685 } 686 if (set && pset && flag && pflag) 687 symmetric = PETSC_TRUE; 688 else 689 symmetric = PETSC_FALSE; 690 691 ierr = VecDot(lclP->lamda0, tao->constraints, &cdotl); CHKERRQ(ierr); 692 lclP->lgn = *f - cdotl; 693 694 /* Gradient of Lagrangian GL = G - J' * lamda */ 695 /* WU = A' * WL 696 WV = B' * WL */ 697 if (symmetric) { 698 ierr = MatMult(tao->jacobian_state,lclP->lamda0,lclP->GL_U); CHKERRQ(ierr); 699 } else { 700 ierr = MatMultTranspose(tao->jacobian_state,lclP->lamda0,lclP->GL_U); CHKERRQ(ierr); 701 } 702 ierr = MatMultTranspose(tao->jacobian_design,lclP->lamda0,lclP->GL_V); CHKERRQ(ierr); 703 ierr = VecScale(lclP->GL_U,-1.0); CHKERRQ(ierr); 704 ierr = VecScale(lclP->GL_V,-1.0); CHKERRQ(ierr); 705 ierr = VecAXPY(lclP->GL_U,1.0,lclP->GU); CHKERRQ(ierr); 706 ierr = VecAXPY(lclP->GL_V,1.0,lclP->GV); CHKERRQ(ierr); 707 ierr = LCLGather(lclP,lclP->GL_U,lclP->GL_V,lclP->GL); CHKERRQ(ierr); 708 709 f[0] = lclP->lgn; 710 ierr = VecCopy(lclP->GL,G); 711 712 PetscFunctionReturn(0); 713 } 714 715 #undef __FUNCT__ 716 #define __FUNCT__ "LCLComputeAugmentedLagrangianAndGradient" 717 static PetscErrorCode LCLComputeAugmentedLagrangianAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *ptr) 718 { 719 TaoSolver tao = (TaoSolver)ptr; 720 TAO_LCL *lclP = (TAO_LCL*)tao->data; 721 PetscReal con2; 722 PetscBool flag,pflag,set,pset,symmetric; 723 PetscErrorCode ierr; 724 725 PetscFunctionBegin; 726 ierr = LCLComputeLagrangianAndGradient(tao->linesearch,X,f,G,tao); CHKERRQ(ierr); 727 ierr = LCLScatter(lclP,G,lclP->GL_U,lclP->GL_V); CHKERRQ(ierr); 728 ierr = VecDot(tao->constraints,tao->constraints,&con2); CHKERRQ(ierr); 729 lclP->aug = lclP->lgn + 0.5*lclP->rho*con2; 730 731 /* Gradient of Aug. Lagrangian GAugL = GL + rho * J' c */ 732 /* WU = A' * c 733 WV = B' * c */ 734 ierr = MatIsSymmetricKnown(tao->jacobian_state,&set,&flag); CHKERRQ(ierr); 735 if (tao->jacobian_state_pre) { 736 ierr = MatIsSymmetricKnown(tao->jacobian_state_pre,&pset,&pflag); CHKERRQ(ierr); 737 } else { 738 pset = pflag = PETSC_TRUE; 739 } 740 if (set && pset && flag && pflag) 741 symmetric = PETSC_TRUE; 742 else 743 symmetric = PETSC_FALSE; 744 745 if (symmetric) { 746 ierr = MatMult(tao->jacobian_state,tao->constraints,lclP->GAugL_U); CHKERRQ(ierr); 747 } else { 748 ierr = MatMultTranspose(tao->jacobian_state,tao->constraints,lclP->GAugL_U); CHKERRQ(ierr); 749 } 750 751 ierr = MatMultTranspose(tao->jacobian_design,tao->constraints,lclP->GAugL_V); CHKERRQ(ierr); 752 ierr = VecAYPX(lclP->GAugL_U,lclP->rho,lclP->GL_U); CHKERRQ(ierr); 753 ierr = VecAYPX(lclP->GAugL_V,lclP->rho,lclP->GL_V); CHKERRQ(ierr); 754 ierr = LCLGather(lclP,lclP->GAugL_U,lclP->GAugL_V,lclP->GAugL); CHKERRQ(ierr); 755 756 f[0] = lclP->aug; 757 ierr = VecCopy(lclP->GAugL,G); CHKERRQ(ierr); 758 759 PetscFunctionReturn(0); 760 } 761 762 #undef __FUNCT__ 763 #define __FUNCT__ "LCLGather" 764 PetscErrorCode LCLGather(TAO_LCL *lclP, Vec u, Vec v, Vec x) 765 { 766 PetscErrorCode ierr; 767 PetscFunctionBegin; 768 ierr = VecScatterBegin(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 769 ierr = VecScatterEnd(lclP->state_scatter, u, x, INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 770 ierr = VecScatterBegin(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 771 ierr = VecScatterEnd(lclP->design_scatter, v, x, INSERT_VALUES, SCATTER_REVERSE); CHKERRQ(ierr); 772 PetscFunctionReturn(0); 773 774 } 775 #undef __FUNCT__ 776 #define __FUNCT__ "LCLScatter" 777 PetscErrorCode LCLScatter(TAO_LCL *lclP, Vec x, Vec u, Vec v) 778 { 779 PetscErrorCode ierr; 780 PetscFunctionBegin; 781 ierr = VecScatterBegin(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); 782 ierr = VecScatterEnd(lclP->state_scatter, x, u, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); 783 ierr = VecScatterBegin(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); 784 ierr = VecScatterEnd(lclP->design_scatter, x, v, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr); 785 PetscFunctionReturn(0); 786 787 } 788