1 #include <../src/tao/complementarity/impls/ssls/ssls.h> 2 3 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(0); 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 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 52 if (tao->reason!=TAO_CONTINUE_ITERATING) break; 53 54 /* Call general purpose update function */ 55 if (tao->ops->update) { 56 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 57 } 58 tao->niter++; 59 60 /* Calculate direction. (Really negative of newton direction. Therefore, 61 rest of the code uses -d.) */ 62 PetscCall(KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre)); 63 PetscCall(KSPSolve(tao->ksp,ssls->ff,tao->stepdirection)); 64 PetscCall(KSPGetIterationNumber(tao->ksp,&tao->ksp_its)); 65 tao->ksp_tot_its+=tao->ksp_its; 66 67 PetscCall(VecCopy(tao->stepdirection,ssls->w)); 68 PetscCall(VecScale(ssls->w,-1.0)); 69 PetscCall(VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w)); 70 71 PetscCall(VecNorm(ssls->w,NORM_2,&normd)); 72 PetscCall(VecDot(ssls->w,ssls->dpsi,&innerd)); 73 74 /* Make sure that we have a descent direction */ 75 if (innerd >= -delta*PetscPowReal(normd, rho)) { 76 PetscCall(PetscInfo(tao, "newton direction not descent\n")); 77 PetscCall(VecCopy(ssls->dpsi,tao->stepdirection)); 78 PetscCall(VecDot(ssls->w,ssls->dpsi,&innerd)); 79 } 80 81 PetscCall(VecScale(tao->stepdirection, -1.0)); 82 innerd = -innerd; 83 84 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 85 PetscCall(TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason)); 86 PetscCall(VecNorm(ssls->dpsi,NORM_2,&ndpsi)); 87 } 88 PetscFunctionReturn(0); 89 } 90 91 PetscErrorCode TaoDestroy_SSFLS(Tao tao) 92 { 93 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 94 95 PetscFunctionBegin; 96 PetscCall(VecDestroy(&ssls->ff)); 97 PetscCall(VecDestroy(&ssls->w)); 98 PetscCall(VecDestroy(&ssls->dpsi)); 99 PetscCall(VecDestroy(&ssls->da)); 100 PetscCall(VecDestroy(&ssls->db)); 101 PetscCall(VecDestroy(&ssls->t1)); 102 PetscCall(VecDestroy(&ssls->t2)); 103 PetscCall(KSPDestroy(&tao->ksp)); 104 PetscCall(PetscFree(tao->data)); 105 PetscFunctionReturn(0); 106 } 107 108 /* ---------------------------------------------------------- */ 109 /*MC 110 TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving 111 complementarity constraints 112 113 Options Database Keys: 114 + -tao_ssls_delta - descent test fraction 115 - -tao_ssls_rho - descent test power 116 117 Level: beginner 118 M*/ 119 120 PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao) 121 { 122 TAO_SSLS *ssls; 123 const char *armijo_type = TAOLINESEARCHARMIJO; 124 125 PetscFunctionBegin; 126 PetscCall(PetscNewLog(tao,&ssls)); 127 tao->data = (void*)ssls; 128 tao->ops->solve=TaoSolve_SSFLS; 129 tao->ops->setup=TaoSetUp_SSFLS; 130 tao->ops->view=TaoView_SSLS; 131 tao->ops->setfromoptions = TaoSetFromOptions_SSLS; 132 tao->ops->destroy = TaoDestroy_SSFLS; 133 134 ssls->delta = 1e-10; 135 ssls->rho = 2.1; 136 137 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch)); 138 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 139 PetscCall(TaoLineSearchSetType(tao->linesearch,armijo_type)); 140 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 141 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 142 /* Linesearch objective and objectivegradient routines are set in solve routine */ 143 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 144 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 145 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 146 147 /* Override default settings (unless already changed) */ 148 if (!tao->max_it_changed) tao->max_it = 2000; 149 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 150 if (!tao->gttol_changed) tao->gttol = 0; 151 if (!tao->grtol_changed) tao->grtol = 0; 152 #if defined(PETSC_USE_REAL_SINGLE) 153 if (!tao->gatol_changed) tao->gatol = 1.0e-6; 154 if (!tao->fmin_changed) tao->fmin = 1.0e-4; 155 #else 156 if (!tao->gatol_changed) tao->gatol = 1.0e-16; 157 if (!tao->fmin_changed) tao->fmin = 1.0e-8; 158 #endif 159 PetscFunctionReturn(0); 160 } 161