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