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