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 should be done with a conjugate gradient 9 method, although any method can be used. 10 */ 11 12 static PetscErrorCode TaoSolve_BNTR(Tao tao) 13 { 14 PetscErrorCode ierr; 15 TAO_BNK *bnk = (TAO_BNK *)tao->data; 16 17 PetscReal oldTrust; 18 PetscBool stepAccepted = PETSC_TRUE; 19 PetscInt stepType; 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 raw Newton step */ 67 ierr = TaoBNKComputeStep(tao, PETSC_FALSE, &stepType);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 /* Test the new step for acceptance */ 77 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 78 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 79 ierr = TaoBNKUpdateTrustRadius(tao, bnk->fold, bnk->f, stepType, &stepAccepted);CHKERRQ(ierr); 80 81 if (stepAccepted) { 82 /* Step is good, evaluate the gradient and the hessian */ 83 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 84 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 85 } else { 86 /* Step is bad, revert old solution and re-solve with new radius*/ 87 bnk->f = bnk->fold; 88 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 89 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 90 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 91 if (oldTrust == tao->trust == bnk->min_radius) { 92 /* Can't shrink trust radius any further, so we have to terminate */ 93 tao->reason = TAO_DIVERGED_TR_REDUCTION; 94 } 95 } 96 97 /* Check for termination */ 98 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 99 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 100 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 101 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 102 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 103 } 104 PetscFunctionReturn(0); 105 } 106 107 PETSC_EXTERN PetscErrorCode TaoCreate_BNTR(Tao tao) 108 { 109 TAO_BNK *bnk; 110 PetscErrorCode ierr; 111 112 PetscFunctionBegin; 113 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 114 tao->ops->solve=TaoSolve_BNTR; 115 116 bnk = (TAO_BNK *)tao->data; 117 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 118 bnk->sval = 0.0; /* disable Hessian shifting */ 119 PetscFunctionReturn(0); 120 }