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