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