1 #include <../src/tao/bound/impls/bnk/bnk.h> 2 #include <petscksp.h> 3 4 /* 5 Implements Newton's Method with a line search approach for solving 6 bound constrained minimization problems. A projected More'-Thuente line 7 search is used to guarantee that the BFGS preconditioner remains positive 8 definite. 9 10 The method can shift the Hessian matrix. The shifting procedure is 11 adapted from the PATH algorithm for solving complementarity 12 problems. 13 14 The linear system solve should be done with a conjugate gradient 15 method, although any method can be used. 16 */ 17 18 static PetscErrorCode TaoSolve_BNLS(Tao tao) 19 { 20 PetscErrorCode ierr; 21 TAO_BNK *bnk = (TAO_BNK *)tao->data; 22 KSPConvergedReason ksp_reason; 23 TaoLineSearchConvergedReason ls_reason; 24 25 PetscReal f_full, prered, actred; 26 PetscReal steplen = 1.0; 27 28 PetscBool trustAccept; 29 PetscInt stepType; 30 PetscInt bfgsUpdates = 0; 31 32 PetscFunctionBegin; 33 /* Project the current point onto the feasible set */ 34 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 35 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 36 37 /* Project the initial point onto the feasible region */ 38 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 39 40 /* Check convergence criteria */ 41 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient);CHKERRQ(ierr); 42 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 43 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 44 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 45 46 tao->reason = TAO_CONTINUE_ITERATING; 47 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 48 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 49 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 50 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 51 52 /* Initialize the preconditioner and trust radius */ 53 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 54 55 /* Have not converged; continue with Newton method */ 56 while (tao->reason == TAO_CONTINUE_ITERATING) { 57 ++tao->niter; 58 tao->ksp_its=0; 59 60 /* Compute the Hessian */ 61 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 62 /* Shift the Hessian matrix */ 63 if (bnk->pert > 0) { 64 ierr = MatShift(tao->hessian, bnk->pert);CHKERRQ(ierr); 65 if (tao->hessian != tao->hessian_pre) { 66 ierr = MatShift(tao->hessian_pre, bnk->pert);CHKERRQ(ierr); 67 } 68 } 69 /* Update the BFGS preconditioner */ 70 if (BNK_PC_BFGS == bnk->pc_type) { 71 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 72 /* Obtain diagonal for the bfgs preconditioner */ 73 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 74 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 75 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 76 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 77 } 78 /* Update the limited memory preconditioner and get existing # of updates */ 79 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 80 ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr); 81 } 82 83 /* Use the common BNK kernel to compute the safeguarded Newton step (for inactive variables only) */ 84 ierr = TaoBNKComputeStep(tao, &ksp_reason);CHKERRQ(ierr); 85 ierr = TaoBNKSafeguardStep(tao, ksp_reason, &stepType);CHKERRQ(ierr); 86 87 /* Store current solution before it changes */ 88 bnk->fold = bnk->f; 89 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 90 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 91 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 92 93 /* Trigger the line search */ 94 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 95 96 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 97 /* Failed to find an improving point */ 98 bnk->f = bnk->fold; 99 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 100 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 101 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 102 steplen = 0.0; 103 tao->reason = TAO_DIVERGED_LS_FAILURE; 104 break; 105 } else { 106 /* count the accepted step types */ 107 switch (stepType) { 108 case BNK_NEWTON: 109 ++bnk->newt; 110 break; 111 case BNK_BFGS: 112 ++bnk->bfgs; 113 break; 114 case BNK_SCALED_GRADIENT: 115 ++bnk->sgrad; 116 break; 117 case BNK_GRADIENT: 118 ++bnk->grad; 119 break; 120 default: 121 break; 122 } 123 } 124 125 /* update trust radius */ 126 ierr = TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full);CHKERRQ(ierr); 127 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 128 prered = -prered; 129 actred = bnk->fold - f_full; 130 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &trustAccept);CHKERRQ(ierr); 131 132 /* Check for termination */ 133 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 134 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 135 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 136 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 137 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 138 } 139 PetscFunctionReturn(0); 140 } 141 142 /*------------------------------------------------------------*/ 143 144 PETSC_EXTERN PetscErrorCode TaoCreate_BNLS(Tao tao) 145 { 146 TAO_BNK *bnk; 147 PetscErrorCode ierr; 148 149 PetscFunctionBegin; 150 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 151 tao->ops->solve = TaoSolve_BNLS; 152 153 bnk = (TAO_BNK *)tao->data; 154 bnk->update_type = BNK_UPDATE_STEP; /* trust region updates based on line search step length */ 155 PetscFunctionReturn(0); 156 }