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