1 #include <../src/tao/bound/impls/bnk/bnk.h> 2 #include <petscksp.h> 3 4 /* 5 Implements Newton's Method with a trust region approach for solving 6 bound constrained minimization problems. 7 8 The linear system solve has to be done with a conjugate gradient method. 9 */ 10 11 static PetscErrorCode TaoSolve_BNTR(Tao tao) 12 { 13 PetscErrorCode ierr; 14 TAO_BNK *bnk = (TAO_BNK *)tao->data; 15 KSPConvergedReason ksp_reason; 16 17 PetscReal oldTrust, prered, actred, stepNorm, steplen; 18 PetscBool stepAccepted = PETSC_TRUE; 19 PetscInt stepType = BNK_NEWTON; 20 21 PetscFunctionBegin; 22 /* Project the current point onto the feasible set */ 23 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 24 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 25 26 /* Project the initial point onto the feasible region */ 27 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 28 29 /* Check convergence criteria */ 30 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient);CHKERRQ(ierr); 31 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 32 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 33 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 34 35 tao->reason = TAO_CONTINUE_ITERATING; 36 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 37 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 38 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 39 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 40 41 /* Initialize the preconditioner and trust radius */ 42 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 43 44 /* Have not converged; continue with Newton method */ 45 while (tao->reason == TAO_CONTINUE_ITERATING) { 46 47 if (stepAccepted) { 48 tao->niter++; 49 tao->ksp_its=0; 50 /* Compute the Hessian */ 51 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 52 /* Update the BFGS preconditioner */ 53 if (BNK_PC_BFGS == bnk->pc_type) { 54 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 55 /* Obtain diagonal for the bfgs preconditioner */ 56 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 57 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 58 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 59 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 60 } 61 /* Update the limited memory preconditioner and get existing # of updates */ 62 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 63 } 64 } 65 66 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 67 ierr = TaoBNKComputeStep(tao, &ksp_reason);CHKERRQ(ierr); 68 69 /* Store current solution before it changes */ 70 oldTrust = tao->trust; 71 bnk->fold = bnk->f; 72 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 73 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 74 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 75 76 /* Temporarily accept the step and project it into the bounds */ 77 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 78 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 79 80 /* Check if the projection changed the step direction */ 81 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 82 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 83 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 84 if (stepNorm != bnk->dnorm) { 85 /* Projection changed the step, so we have to recompute predicted reduction. 86 However, we deliberately do not change the step norm and the trust radius 87 in order for the safeguard to more closely mimic a piece-wise linesearch 88 along the bounds. */ 89 ierr = MatMult(tao->hessian, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 90 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 91 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 92 } else { 93 /* Step did not change, so we can just recover the pre-computed prediction */ 94 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 95 } 96 prered = -prered; 97 98 /* Compute the actual reduction and update the trust radius */ 99 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 100 actred = bnk->fold - bnk->f; 101 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 102 103 if (stepAccepted) { 104 /* Step is good, evaluate the gradient and the hessian */ 105 steplen = 1.0; 106 ++bnk->newt; 107 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 108 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 109 } else { 110 /* Step is bad, revert old solution and re-solve with new radius*/ 111 steplen = 0.0; 112 bnk->f = bnk->fold; 113 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 114 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 115 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 116 if (oldTrust == tao->trust) { 117 /* Can't change the radius anymore so just terminate */ 118 tao->reason = TAO_DIVERGED_TR_REDUCTION; 119 } 120 } 121 122 /* Check for termination */ 123 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 124 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 125 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 126 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 127 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 128 } 129 PetscFunctionReturn(0); 130 } 131 132 /*------------------------------------------------------------*/ 133 134 PETSC_EXTERN PetscErrorCode TaoCreate_BNTR(Tao tao) 135 { 136 TAO_BNK *bnk; 137 PetscErrorCode ierr; 138 139 PetscFunctionBegin; 140 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 141 tao->ops->solve=TaoSolve_BNTR; 142 143 bnk = (TAO_BNK *)tao->data; 144 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 145 bnk->sval = 0.0; /* disable Hessian shifting */ 146 PetscFunctionReturn(0); 147 }