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