1 #include <../src/tao/complementarity/impls/ssls/ssls.h> 2 3 PetscErrorCode TaoSetUp_SSILS(Tao tao) 4 { 5 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 6 7 PetscFunctionBegin; 8 PetscCall(VecDuplicate(tao->solution,&tao->gradient)); 9 PetscCall(VecDuplicate(tao->solution,&tao->stepdirection)); 10 PetscCall(VecDuplicate(tao->solution,&ssls->ff)); 11 PetscCall(VecDuplicate(tao->solution,&ssls->dpsi)); 12 PetscCall(VecDuplicate(tao->solution,&ssls->da)); 13 PetscCall(VecDuplicate(tao->solution,&ssls->db)); 14 PetscCall(VecDuplicate(tao->solution,&ssls->t1)); 15 PetscCall(VecDuplicate(tao->solution,&ssls->t2)); 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode TaoDestroy_SSILS(Tao tao) 20 { 21 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 22 23 PetscFunctionBegin; 24 PetscCall(VecDestroy(&ssls->ff)); 25 PetscCall(VecDestroy(&ssls->dpsi)); 26 PetscCall(VecDestroy(&ssls->da)); 27 PetscCall(VecDestroy(&ssls->db)); 28 PetscCall(VecDestroy(&ssls->t1)); 29 PetscCall(VecDestroy(&ssls->t2)); 30 PetscCall(PetscFree(tao->data)); 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode TaoSolve_SSILS(Tao tao) 35 { 36 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 37 PetscReal psi, ndpsi, normd, innerd, t=0; 38 PetscReal delta, rho; 39 TaoLineSearchConvergedReason ls_reason; 40 41 PetscFunctionBegin; 42 /* Assume that Setup has been called! 43 Set the structure for the Jacobian and create a linear solver. */ 44 delta = ssls->delta; 45 rho = ssls->rho; 46 47 PetscCall(TaoComputeVariableBounds(tao)); 48 PetscCall(VecMedian(tao->XL,tao->solution,tao->XU,tao->solution)); 49 PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao)); 50 PetscCall(TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao)); 51 52 /* Calculate the function value and fischer function value at the 53 current iterate */ 54 PetscCall(TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi)); 55 PetscCall(VecNorm(ssls->dpsi,NORM_2,&ndpsi)); 56 57 tao->reason = TAO_CONTINUE_ITERATING; 58 while (PETSC_TRUE) { 59 PetscCall(PetscInfo(tao, "iter: %D, merit: %g, ndpsi: %g\n",tao->niter, (double)ssls->merit, (double)ndpsi)); 60 /* Check the termination criteria */ 61 PetscCall(TaoLogConvergenceHistory(tao,ssls->merit,ndpsi,0.0,tao->ksp_its)); 62 PetscCall(TaoMonitor(tao,tao->niter,ssls->merit,ndpsi,0.0,t)); 63 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 64 if (tao->reason!=TAO_CONTINUE_ITERATING) break; 65 66 /* Call general purpose update function */ 67 if (tao->ops->update) { 68 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 69 } 70 tao->niter++; 71 72 /* Calculate direction. (Really negative of newton direction. Therefore, 73 rest of the code uses -d.) */ 74 PetscCall(KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre)); 75 PetscCall(KSPSolve(tao->ksp,ssls->ff,tao->stepdirection)); 76 PetscCall(KSPGetIterationNumber(tao->ksp,&tao->ksp_its)); 77 tao->ksp_tot_its+=tao->ksp_its; 78 PetscCall(VecNorm(tao->stepdirection,NORM_2,&normd)); 79 PetscCall(VecDot(tao->stepdirection,ssls->dpsi,&innerd)); 80 81 /* Make sure that we have a descent direction */ 82 if (innerd <= delta*PetscPowReal(normd, rho)) { 83 PetscCall(PetscInfo(tao, "newton direction not descent\n")); 84 PetscCall(VecCopy(ssls->dpsi,tao->stepdirection)); 85 PetscCall(VecDot(tao->stepdirection,ssls->dpsi,&innerd)); 86 } 87 88 PetscCall(VecScale(tao->stepdirection, -1.0)); 89 innerd = -innerd; 90 91 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 92 PetscCall(TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason)); 93 PetscCall(VecNorm(ssls->dpsi,NORM_2,&ndpsi)); 94 } 95 PetscFunctionReturn(0); 96 } 97 98 /* ---------------------------------------------------------- */ 99 /*MC 100 TAOSSILS - semi-smooth infeasible linesearch algorithm for solving 101 complementarity constraints 102 103 Options Database Keys: 104 + -tao_ssls_delta - descent test fraction 105 - -tao_ssls_rho - descent test power 106 107 Level: beginner 108 M*/ 109 PETSC_EXTERN PetscErrorCode TaoCreate_SSILS(Tao tao) 110 { 111 TAO_SSLS *ssls; 112 const char *armijo_type = TAOLINESEARCHARMIJO; 113 114 PetscFunctionBegin; 115 PetscCall(PetscNewLog(tao,&ssls)); 116 tao->data = (void*)ssls; 117 tao->ops->solve=TaoSolve_SSILS; 118 tao->ops->setup=TaoSetUp_SSILS; 119 tao->ops->view=TaoView_SSLS; 120 tao->ops->setfromoptions=TaoSetFromOptions_SSLS; 121 tao->ops->destroy = TaoDestroy_SSILS; 122 123 ssls->delta = 1e-10; 124 ssls->rho = 2.1; 125 126 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch)); 127 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 128 PetscCall(TaoLineSearchSetType(tao->linesearch,armijo_type)); 129 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 130 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 131 /* Note: linesearch objective and objectivegradient routines are set in solve routine */ 132 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 133 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 134 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 135 136 /* Override default settings (unless already changed) */ 137 if (!tao->max_it_changed) tao->max_it = 2000; 138 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 139 if (!tao->gttol_changed) tao->gttol = 0; 140 if (!tao->grtol_changed) tao->grtol = 0; 141 #if defined(PETSC_USE_REAL_SINGLE) 142 if (!tao->gatol_changed) tao->gatol = 1.0e-6; 143 if (!tao->fmin_changed) tao->fmin = 1.0e-4; 144 #else 145 if (!tao->gatol_changed) tao->gatol = 1.0e-16; 146 if (!tao->fmin_changed) tao->fmin = 1.0e-8; 147 #endif 148 PetscFunctionReturn(0); 149 } 150