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