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 method can shift the Hessian matrix. The shifting procedure is 9 adapted from the PATH algorithm for solving complementarity 10 problems. 11 12 The linear system solve should be done with a conjugate gradient 13 method, although any method can be used. 14 */ 15 16 static PetscErrorCode TaoSolve_BNTR(Tao tao) 17 { 18 PetscErrorCode ierr; 19 TAO_BNK *bnk = (TAO_BNK *)tao->data; 20 21 PetscReal oldTrust; 22 PetscBool stepAccepted = PETSC_TRUE; 23 PetscInt stepType; 24 25 PetscFunctionBegin; 26 /* Project the current point onto the feasible set */ 27 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 28 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 29 30 /* Project the initial point onto the feasible region */ 31 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 32 33 /* Check convergence criteria */ 34 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient);CHKERRQ(ierr); 35 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 36 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 37 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 38 39 tao->reason = TAO_CONTINUE_ITERATING; 40 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 41 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 42 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 43 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 44 45 /* Initialize the preconditioner and trust radius */ 46 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 47 48 /* Have not converged; continue with Newton method */ 49 while (tao->reason == TAO_CONTINUE_ITERATING) { 50 if (stepAccepted) { 51 tao->niter++; 52 tao->ksp_its=0; 53 /* Compute the Hessian */ 54 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 55 /* Update the BFGS preconditioner */ 56 if (BNK_PC_BFGS == bnk->pc_type) { 57 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 58 /* Obtain diagonal for the bfgs preconditioner */ 59 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 60 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 61 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 62 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 63 } 64 /* Update the limited memory preconditioner and get existing # of updates */ 65 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 66 } 67 } 68 69 /* Use the common BNK kernel to compute the raw Newton step */ 70 ierr = TaoBNKComputeStep(tao, PETSC_FALSE, &stepType);CHKERRQ(ierr); 71 72 /* Store current solution before it changes */ 73 oldTrust = tao->trust; 74 bnk->fold = bnk->f; 75 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 76 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 77 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 78 79 /* Test the new step for acceptance */ 80 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 81 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 82 ierr = TaoBNKUpdateTrustRadius(tao, bnk->fold, bnk->f, stepType, &stepAccepted);CHKERRQ(ierr); 83 84 if (stepAccepted) { 85 /* Step is good, evaluate gradient and project it */ 86 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 87 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 88 } else { 89 /* Step is bad, revert old solution and re-solve with new radius*/ 90 bnk->f = bnk->fold; 91 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 92 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 93 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 94 if (oldTrust == tao->trust) { 95 /* Can't shrink trust radius any further, so we have to terminate */ 96 tao->reason = TAO_DIVERGED_TR_REDUCTION; 97 } 98 } 99 100 /* Check for termination */ 101 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 102 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 103 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 104 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 105 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 PETSC_EXTERN PetscErrorCode TaoCreate_BNTR(Tao tao) 111 { 112 TAO_BNK *bnk; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 117 tao->ops->solve=TaoSolve_BNTR; 118 119 bnk = (TAO_BNK *)tao->data; 120 bnk->update_type = BNK_UPDATE_REDUCTION; 121 bnk->sval = 0.0; 122 PetscFunctionReturn(0); 123 }