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