1 #include "tao-private/taolinesearch_impl.h" 2 #include "gpcglinesearch.h" 3 4 /* ---------------------------------------------------------- */ 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "TaoLineSearchDestroy_GPCG" 8 static PetscErrorCode TaoLineSearchDestroy_GPCG(TaoLineSearch ls) 9 { 10 PetscErrorCode ierr; 11 TAOLINESEARCH_GPCG_CTX *ctx = (TAOLINESEARCH_GPCG_CTX *)ls->data; 12 13 PetscFunctionBegin; 14 if (ctx->W1) { 15 ierr = VecDestroy(&ctx->W1);CHKERRQ(ierr); 16 } 17 if (ctx->W2) { 18 ierr = VecDestroy(&ctx->W2);CHKERRQ(ierr); 19 } 20 if (ctx->Gold) { 21 ierr = VecDestroy(&ctx->Gold);CHKERRQ(ierr); 22 } 23 if (ctx->x) { 24 ierr = PetscObjectDereference((PetscObject)ctx->x); CHKERRQ(ierr); 25 } 26 27 ierr = PetscFree(ctx);CHKERRQ(ierr); 28 ls->data = 0; 29 PetscFunctionReturn(0); 30 } 31 32 33 /*------------------------------------------------------------*/ 34 #undef __FUNCT__ 35 #define __FUNCT__ "TaoLineSearchView_GPCG" 36 static PetscErrorCode TaoLineSearchView_GPCG(TaoLineSearch ls, PetscViewer viewer) 37 { 38 PetscBool isascii; 39 PetscErrorCode ierr; 40 PetscFunctionBegin; 41 42 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 43 44 if (isascii) { 45 ierr = PetscViewerASCIIPrintf(viewer," GPCG Line search"); CHKERRQ(ierr); 46 } else { 47 SETERRQ1(((PetscObject)ls)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for GPCG LineSearch",((PetscObject)viewer)->type_name); 48 } 49 PetscFunctionReturn(0); 50 } 51 52 /*------------------------------------------------------------*/ 53 #undef __FUNCT__ 54 #define __FUNCT__ "TaoLineSearchApply_GPCG" 55 static PetscErrorCode TaoLineSearchApply_GPCG(TaoLineSearch ls, Vec x, 56 PetscReal *f, Vec g, Vec s) 57 58 { 59 TAOLINESEARCH_GPCG_CTX *neP = (TAOLINESEARCH_GPCG_CTX *)ls->data; 60 PetscErrorCode ierr; 61 PetscInt i; 62 PetscBool g_computed=PETSC_FALSE; /* to prevent extra gradient computation */ 63 PetscReal d1,finit,actred,prered,rho, gdx; 64 65 PetscFunctionBegin; 66 /* ls->stepmin - lower bound for step */ 67 /* ls->stepmax - upper bound for step */ 68 /* ls->rtol - relative tolerance for an acceptable step */ 69 /* ls->ftol - tolerance for sufficient decrease condition */ 70 /* ls->gtol - tolerance for curvature condition */ 71 /* ls->nfeval - number of function evaluations */ 72 /* ls->nfeval - number of function/gradient evaluations */ 73 /* ls->max_funcs - maximum number of function evaluations */ 74 75 ls->reason = TAOLINESEARCH_CONTINUE_ITERATING; 76 ls->step = ls->initstep; 77 if (!neP->W2) { 78 ierr = VecDuplicate(x,&neP->W2); CHKERRQ(ierr); 79 ierr = VecDuplicate(x,&neP->W1); CHKERRQ(ierr); 80 ierr = VecDuplicate(x,&neP->Gold); CHKERRQ(ierr); 81 neP->x = x; 82 ierr = PetscObjectReference((PetscObject)neP->x); CHKERRQ(ierr); 83 } 84 85 /* If X has changed, remake work vectors */ 86 else if (x != neP->x) { 87 ierr = VecDestroy(&neP->x); CHKERRQ(ierr); 88 ierr = VecDestroy(&neP->W1); CHKERRQ(ierr); 89 ierr = VecDestroy(&neP->W2); CHKERRQ(ierr); 90 ierr = VecDestroy(&neP->Gold); CHKERRQ(ierr); 91 ierr = VecDuplicate(x,&neP->W1); CHKERRQ(ierr); 92 ierr = VecDuplicate(x,&neP->W2); CHKERRQ(ierr); 93 ierr = VecDuplicate(x,&neP->Gold); CHKERRQ(ierr); 94 ierr = PetscObjectDereference((PetscObject)neP->x); CHKERRQ(ierr); 95 neP->x = x; 96 ierr = PetscObjectReference((PetscObject)neP->x); CHKERRQ(ierr); 97 } 98 99 ierr = VecDot(g,s,&gdx); CHKERRQ(ierr); 100 if (gdx > 0) { 101 ierr = PetscInfo1(ls,"Line search error: search direction is not descent direction. dot(g,s) = %G",gdx); CHKERRQ(ierr); 102 ls->reason = TAOLINESEARCH_FAILED_ASCENT; 103 PetscFunctionReturn(0); 104 } 105 ierr = VecCopy(x,neP->W2); CHKERRQ(ierr); 106 ierr = VecCopy(g,neP->Gold); CHKERRQ(ierr); 107 if (ls->bounded) { 108 /* Compute the smallest steplength that will make one nonbinding variable 109 equal the bound */ 110 ierr = VecStepBoundInfo(x,ls->lower,ls->upper,s,&rho,&actred,&d1); CHKERRQ(ierr); 111 ls->step = PetscMin(ls->step,d1); 112 } 113 rho=0; actred=0; 114 115 if (ls->step < 0) { 116 ierr = PetscInfo1(ls,"Line search error: initial step parameter %G < 0\n",ls->step); CHKERRQ(ierr); 117 ls->reason = TAOLINESEARCH_HALTED_OTHER; 118 PetscFunctionReturn(0); 119 } 120 121 /* Initialization */ 122 finit = *f; 123 for (i=0; i< ls->max_funcs; i++) { 124 125 /* Force the step to be within the bounds */ 126 ls->step = PetscMax(ls->step,ls->stepmin); 127 ls->step = PetscMin(ls->step,ls->stepmax); 128 129 ierr = VecCopy(x,neP->W2); CHKERRQ(ierr); 130 ierr = VecAXPY(neP->W2,ls->step,s); CHKERRQ(ierr); 131 if (ls->bounded) { 132 /* Make sure new vector is numerically within bounds */ 133 ierr = VecMedian(neP->W2,ls->lower,ls->upper,neP->W2);CHKERRQ(ierr); 134 135 } 136 137 /* Gradient is not needed here. Unless there is a separate 138 gradient routine, compute it here anyway to prevent recomputing at 139 the end of the line search */ 140 if (ls->hasobjective) { 141 ierr = TaoLineSearchComputeObjective(ls,neP->W2,f); CHKERRQ(ierr); 142 g_computed=PETSC_FALSE; 143 } else if (ls->usegts){ 144 ierr = TaoLineSearchComputeObjectiveAndGTS(ls,neP->W2,f,&gdx); CHKERRQ(ierr); 145 g_computed=PETSC_FALSE; 146 } else { 147 ierr = TaoLineSearchComputeObjectiveAndGradient(ls,neP->W2,f,g); CHKERRQ(ierr); 148 g_computed=PETSC_TRUE; 149 } 150 151 if (0 == i) { 152 ls->f_fullstep = *f; 153 } 154 155 actred = *f - finit; 156 ierr = VecCopy(neP->W2,neP->W1); CHKERRQ(ierr); 157 ierr = VecAXPY(neP->W1,-1.0,x); CHKERRQ(ierr); /* W1 = W2 - X */ 158 ierr = VecDot(neP->W1,neP->Gold,&prered); CHKERRQ(ierr); 159 160 if (fabs(prered)<1.0e-100) prered=1.0e-12; 161 rho = actred/prered; 162 163 /* 164 If sufficient progress has been obtained, accept the 165 point. Otherwise, backtrack. 166 */ 167 168 if (actred > 0) { 169 ierr = PetscInfo(ls,"Step resulted in ascent, rejecting.\n"); CHKERRQ(ierr); 170 ls->step = (ls->step)/2; 171 } else if (rho > ls->ftol){ 172 break; 173 } else{ 174 ls->step = (ls->step)/2; 175 } 176 177 /* Convergence testing */ 178 179 if (ls->step <= ls->stepmin || ls->step >= ls->stepmax) { 180 ls->reason = TAOLINESEARCH_HALTED_OTHER; 181 ierr = PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n"); CHKERRQ(ierr); 182 ierr = PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n"); CHKERRQ(ierr); 183 break; 184 } 185 if (ls->step == ls->stepmax) { 186 ierr = PetscInfo1(ls,"Step is at the upper bound, stepmax (%G)\n",ls->stepmax); CHKERRQ(ierr); 187 ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND; 188 break; 189 } 190 if (ls->step == ls->stepmin) { 191 ierr = PetscInfo1(ls,"Step is at the lower bound, stepmin (%G)\n",ls->stepmin); CHKERRQ(ierr); 192 ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND; 193 break; 194 } 195 if ((ls->nfeval+ls->nfgeval) >= ls->max_funcs) { 196 ierr = PetscInfo2(ls,"Number of line search function evals (%D) > maximum (%D)\n",ls->nfeval+ls->nfgeval,ls->max_funcs); CHKERRQ(ierr); 197 ls->reason = TAOLINESEARCH_HALTED_MAXFCN; 198 break; 199 } 200 if ((neP->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)){ 201 ierr = PetscInfo1(ls,"Relative width of interval of uncertainty is at most rtol (%G)\n",ls->rtol); CHKERRQ(ierr); 202 ls->reason = TAOLINESEARCH_HALTED_RTOL; 203 break; 204 } 205 } 206 ierr = PetscInfo2(ls,"%D function evals in line search, step = %G\n",ls->nfeval+ls->nfgeval,ls->step); CHKERRQ(ierr); 207 /* set new solution vector and compute gradient if necessary */ 208 ierr = VecCopy(neP->W2, x); CHKERRQ(ierr); 209 if (!g_computed) { 210 ierr = TaoLineSearchComputeGradient(ls,x,g); CHKERRQ(ierr); 211 } 212 PetscFunctionReturn(0); 213 } 214 215 /* ---------------------------------------------------------- */ 216 EXTERN_C_BEGIN 217 #undef __FUNCT__ 218 #define __FUNCT__ "TaoLineSearchCreate_GPCG" 219 PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls) 220 { 221 PetscErrorCode ierr; 222 TAOLINESEARCH_GPCG_CTX *neP; 223 224 PetscFunctionBegin; 225 226 227 228 ls->ftol = 0.05; 229 ls->rtol = 0.0; 230 ls->gtol = 0.0; 231 ls->stepmin = 1.0e-20; 232 ls->stepmax = 1.0e+20; 233 ls->nfeval = 0; 234 ls->max_funcs = 30; 235 ls->step = 1.0; 236 237 ierr = PetscNewLog(ls,&neP);CHKERRQ(ierr); 238 neP->bracket = 0; 239 neP->infoc = 1; 240 ls->data = (void*)neP; 241 242 ls->ops->setup = 0; 243 ls->ops->apply=TaoLineSearchApply_GPCG; 244 ls->ops->view =TaoLineSearchView_GPCG; 245 ls->ops->destroy=TaoLineSearchDestroy_GPCG; 246 ls->ops->setfromoptions=0; 247 248 PetscFunctionReturn(0); 249 } 250 EXTERN_C_END 251