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 steplen = 1.0; 26 PetscBool shift = PETSC_TRUE; 27 PetscInt stepType; 28 29 PetscFunctionBegin; 30 /* Initialize the preconditioner, KSP solver and trust radius/line search */ 31 tao->reason = TAO_CONTINUE_ITERATING; 32 ierr = TaoBNKInitialize(tao, BNK_INIT_CONSTANT);CHKERRQ(ierr); 33 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 34 35 /* Have not converged; continue with Newton method */ 36 while (tao->reason == TAO_CONTINUE_ITERATING) { 37 ++tao->niter; 38 tao->ksp_its=0; 39 40 /* Compute the hessian and update the BFGS preconditioner at the new iterate*/ 41 ierr = TaoBNKComputeHessian(tao);CHKERRQ(ierr); 42 43 /* Use the common BNK kernel to compute the safeguarded Newton step (for inactive variables only) */ 44 tao->trust = bnk->max_radius; 45 ierr = TaoBNKComputeStep(tao, shift, &ksp_reason);CHKERRQ(ierr); 46 ierr = TaoBNKSafeguardStep(tao, ksp_reason, &stepType);CHKERRQ(ierr); 47 48 /* Store current solution before it changes */ 49 bnk->fold = bnk->f; 50 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 51 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 52 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 53 54 /* Trigger the line search */ 55 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 56 57 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 58 /* Failed to find an improving point */ 59 bnk->f = bnk->fold; 60 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 61 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 62 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 63 steplen = 0.0; 64 tao->reason = TAO_DIVERGED_LS_FAILURE; 65 break; 66 } else { 67 /* count the accepted step type */ 68 ierr = TaoBNKAddStepCounts(tao, stepType);CHKERRQ(ierr); 69 } 70 71 /* Check for termination */ 72 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 73 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 74 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 75 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 76 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 77 } 78 PetscFunctionReturn(0); 79 } 80 81 /*------------------------------------------------------------*/ 82 83 PETSC_EXTERN PetscErrorCode TaoCreate_BNLS(Tao tao) 84 { 85 TAO_BNK *bnk; 86 PetscErrorCode ierr; 87 88 PetscFunctionBegin; 89 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 90 tao->ops->solve = TaoSolve_BNLS; 91 92 bnk = (TAO_BNK *)tao->data; 93 bnk->init_type = BNK_INIT_CONSTANT; 94 bnk->update_type = BNK_UPDATE_STEP; /* trust region updates based on line search step length */ 95 PetscFunctionReturn(0); 96 }