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