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