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