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 /* Initialize the preconditioner, KSP solver and trust radius/line search */ 34 tao->reason = TAO_CONTINUE_ITERATING; 35 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 36 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 37 38 /* Have not converged; continue with Newton method */ 39 while (tao->reason == TAO_CONTINUE_ITERATING) { 40 ++tao->niter; 41 tao->ksp_its=0; 42 43 /* Compute the Hessian */ 44 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 45 /* Shift the Hessian matrix */ 46 if (bnk->pert > 0) { 47 ierr = MatShift(tao->hessian, bnk->pert);CHKERRQ(ierr); 48 if (tao->hessian != tao->hessian_pre) { 49 ierr = MatShift(tao->hessian_pre, bnk->pert);CHKERRQ(ierr); 50 } 51 } 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 ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr); 64 } 65 66 /* Use the common BNK kernel to compute the safeguarded Newton step (for inactive variables only) */ 67 tao->trust = bnk->max_radius; 68 ierr = TaoBNKComputeStep(tao, &ksp_reason);CHKERRQ(ierr); 69 ierr = TaoBNKSafeguardStep(tao, ksp_reason, &stepType);CHKERRQ(ierr); 70 71 /* Store current solution before it changes */ 72 bnk->fold = bnk->f; 73 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 74 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 75 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 76 77 /* Trigger the line search */ 78 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 79 80 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 81 /* Failed to find an improving point */ 82 bnk->f = bnk->fold; 83 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 84 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 85 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 86 steplen = 0.0; 87 tao->reason = TAO_DIVERGED_LS_FAILURE; 88 break; 89 } else { 90 /* count the accepted step types */ 91 switch (stepType) { 92 case BNK_NEWTON: 93 ++bnk->newt; 94 break; 95 case BNK_BFGS: 96 ++bnk->bfgs; 97 break; 98 case BNK_SCALED_GRADIENT: 99 ++bnk->sgrad; 100 break; 101 case BNK_GRADIENT: 102 ++bnk->grad; 103 break; 104 default: 105 break; 106 } 107 } 108 109 /* Check for termination */ 110 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 111 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 112 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 113 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 114 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 115 } 116 PetscFunctionReturn(0); 117 } 118 119 /*------------------------------------------------------------*/ 120 121 PETSC_EXTERN PetscErrorCode TaoCreate_BNLS(Tao tao) 122 { 123 TAO_BNK *bnk; 124 PetscErrorCode ierr; 125 126 PetscFunctionBegin; 127 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 128 tao->ops->solve = TaoSolve_BNLS; 129 130 bnk = (TAO_BNK *)tao->data; 131 bnk->init_type = BNK_INIT_CONSTANT; 132 bnk->update_type = BNK_UPDATE_STEP; /* trust region updates based on line search step length */ 133 PetscFunctionReturn(0); 134 }