1 #include <private/snesimpl.h> 2 3 typedef struct { 4 Vec * dX; /* The change in X */ 5 Vec * dF; /* The change in F */ 6 PetscInt m; /* the number of kept previous steps */ 7 PetscScalar * alpha; 8 PetscScalar * beta; 9 PetscScalar * rho; 10 PetscViewer monitor; 11 } QNContext; 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "LBGFSApplyJinv_Private" 15 PetscErrorCode LBGFSApplyJinv_Private(SNES snes, PetscInt it, Vec g, Vec z) { 16 17 PetscErrorCode ierr; 18 19 QNContext * qn = (QNContext *)snes->data; 20 21 Vec * dX = qn->dX; 22 Vec * dF = qn->dF; 23 24 PetscScalar * alpha = qn->alpha; 25 PetscScalar * beta = qn->beta; 26 PetscScalar * rho = qn->rho; 27 28 PetscInt k, i; 29 PetscInt m = qn->m; 30 PetscScalar t; 31 PetscInt l = m; 32 33 PetscFunctionBegin; 34 35 if (it < m) l = it; 36 37 ierr = VecCopy(g, z);CHKERRQ(ierr); 38 39 /* outward recursion starting at iteration k's update and working back */ 40 for (i = 0; i < l; i++) { 41 k = (it - i - 1) % l; 42 ierr = VecDot(dX[k], z, &t);CHKERRQ(ierr); 43 alpha[k] = t*rho[k]; 44 if (qn->monitor) { 45 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 46 ierr = PetscViewerASCIIPrintf(qn->monitor, "k: %d alpha: %14.12e\n", k, PetscRealPart(alpha[k]));CHKERRQ(ierr); 47 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 48 } 49 ierr = VecAXPY(z, -alpha[k], dF[k]);CHKERRQ(ierr); 50 } 51 52 /* inner application of the initial inverse jacobian approximation */ 53 /* right now it's just the identity. Nothing needs to go here. */ 54 55 /* inward recursion starting at the first update and working forward*/ 56 for (i = 0; i < l; i++) { 57 k = (it + i - l) % l; 58 ierr = VecDot(dF[k], z, &t);CHKERRQ(ierr); 59 beta[k] = rho[k]*t; 60 ierr = VecAXPY(z, (alpha[k] - beta[k]), dX[k]); 61 if (qn->monitor) { 62 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 63 ierr = PetscViewerASCIIPrintf(qn->monitor, "k: %d alpha - beta: %14.12e\n", k, PetscRealPart(alpha[k] - beta[k]));CHKERRQ(ierr); 64 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 65 } 66 } 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "SNESSolve_QN" 72 static PetscErrorCode SNESSolve_QN(SNES snes) 73 { 74 75 PetscErrorCode ierr; 76 QNContext * qn = (QNContext*) snes->data; 77 78 Vec X, Xold; 79 Vec F, Fold, G; 80 Vec W, Y; 81 82 PetscInt i, k; 83 84 PetscReal fnorm, xnorm = 0, ynorm, gnorm; 85 PetscInt m = qn->m; 86 PetscBool lssucceed; 87 88 PetscScalar rhosc; 89 90 Vec * dX = qn->dX; 91 Vec * dF = qn->dF; 92 PetscScalar * rho = qn->rho; 93 94 /* basically just a regular newton's method except for the application of the jacobian */ 95 PetscFunctionBegin; 96 97 X = snes->vec_sol; /* solution vector */ 98 F = snes->vec_func; /* residual vector */ 99 Y = snes->vec_sol_update; /* search direction */ 100 G = snes->work[0]; 101 W = snes->work[1]; 102 Xold = snes->work[2]; 103 Fold = snes->work[3]; 104 105 snes->reason = SNES_CONVERGED_ITERATING; 106 107 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 108 snes->iter = 0; 109 snes->norm = 0.; 110 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 111 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 112 if (snes->domainerror) { 113 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 114 PetscFunctionReturn(0); 115 } 116 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 117 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 118 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 119 snes->norm = fnorm; 120 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 121 SNESLogConvHistory(snes,fnorm,0); 122 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 123 124 /* set parameter for default relative tolerance convergence test */ 125 snes->ttol = fnorm*snes->rtol; 126 /* test convergence */ 127 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 128 if (snes->reason) PetscFunctionReturn(0); 129 130 /* initialize the search direction as steepest descent */ 131 ierr = VecCopy(F, Y);CHKERRQ(ierr); 132 for(i = 0; i < snes->max_its; i++) { 133 /* line search for lambda */ 134 ynorm = 1; gnorm = fnorm; 135 ierr = VecCopy(F, Fold);CHKERRQ(ierr); 136 ierr = VecCopy(X, Xold);CHKERRQ(ierr); 137 ierr = VecScale(Y, -1.0);CHKERRQ(ierr); 138 ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,G,W,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 139 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);CHKERRQ(ierr); 140 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 141 if (snes->domainerror) { 142 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 143 PetscFunctionReturn(0); 144 } 145 if (!lssucceed) { 146 if (++snes->numFailures >= snes->maxFailures) { 147 snes->reason = SNES_DIVERGED_LINE_SEARCH; 148 break; 149 } 150 } 151 /* Update function and solution vectors */ 152 fnorm = gnorm; 153 ierr = VecCopy(G,F);CHKERRQ(ierr); 154 ierr = VecCopy(W,X);CHKERRQ(ierr); 155 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 156 snes->iter = i+1; 157 snes->norm = fnorm; 158 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 159 SNESLogConvHistory(snes,snes->norm,snes->iter); 160 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 161 /* set parameter for default relative tolerance convergence test */ 162 ierr = (*snes->ops->converged)(snes,i+1,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 163 if (snes->reason) PetscFunctionReturn(0); 164 165 /* set the differences */ 166 k = i % m; 167 ierr = VecCopy(F, dF[k]);CHKERRQ(ierr); 168 ierr = VecAXPY(dF[k], -1.0, Fold);CHKERRQ(ierr); 169 ierr = VecCopy(X, dX[k]);CHKERRQ(ierr); 170 ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr); 171 ierr = VecDot(dX[k], dF[k], &rhosc);CHKERRQ(ierr); 172 rho[k] = 1. / rhosc; 173 174 /* general purpose update */ 175 if (snes->ops->update) { 176 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 177 } 178 /* apply the current iteration of the approximate jacobian in order to get the next search direction*/ 179 ierr = LBGFSApplyJinv_Private(snes, i, F, Y);CHKERRQ(ierr); 180 } 181 if (i == snes->max_its) { 182 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 183 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 184 } 185 PetscFunctionReturn(0); 186 } 187 188 189 #undef __FUNCT__ 190 #define __FUNCT__ "SNESSetUp_QN" 191 static PetscErrorCode SNESSetUp_QN(SNES snes) 192 { 193 QNContext * qn = (QNContext *)snes->data; 194 PetscErrorCode ierr; 195 PetscFunctionBegin; 196 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dX);CHKERRQ(ierr); 197 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dF);CHKERRQ(ierr); 198 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, qn->m, PetscScalar, &qn->beta, qn->m, PetscScalar, &qn->rho);CHKERRQ(ierr); 199 ierr = SNESDefaultGetWork(snes,4);CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "SNESReset_QN" 205 static PetscErrorCode SNESReset_QN(SNES snes) 206 { 207 PetscErrorCode ierr; 208 QNContext * qn; 209 PetscFunctionBegin; 210 if (snes->data) { 211 qn = (QNContext *)snes->data; 212 if (qn->dX) { 213 ierr = VecDestroyVecs(qn->m, &qn->dX);CHKERRQ(ierr); 214 } 215 if (qn->dF) { 216 ierr = VecDestroyVecs(qn->m, &qn->dF);CHKERRQ(ierr); 217 } 218 ierr = PetscFree3(qn->alpha, qn->beta, qn->rho);CHKERRQ(ierr); 219 } 220 if (snes->work) {ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr);} 221 PetscFunctionReturn(0); 222 } 223 224 #undef __FUNCT__ 225 #define __FUNCT__ "SNESDestroy_QN" 226 static PetscErrorCode SNESDestroy_QN(SNES snes) 227 { 228 PetscErrorCode ierr; 229 PetscFunctionBegin; 230 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 231 ierr = PetscFree(snes->data);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "SNESSetFromOptions_QN" 237 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 238 { 239 240 PetscErrorCode ierr; 241 QNContext * qn; 242 PetscBool monflg = PETSC_FALSE; 243 PetscFunctionBegin; 244 245 qn = (QNContext *)snes->data; 246 247 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 248 ierr = PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-Broyden methods", "SNES", qn->m, &qn->m, PETSC_NULL);CHKERRQ(ierr); 249 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNES", monflg, &monflg, PETSC_NULL);CHKERRQ(ierr); 250 ierr = PetscOptionsTail();CHKERRQ(ierr); 251 if (monflg) { 252 qn->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 253 } 254 PetscFunctionReturn(0); 255 } 256 257 EXTERN_C_BEGIN 258 #undef __FUNCT__ 259 #define __FUNCT__ "SNESLineSearchSetType_QN" 260 PetscErrorCode SNESLineSearchSetType_QN(SNES snes, SNESLineSearchType type) 261 { 262 PetscErrorCode ierr; 263 PetscFunctionBegin; 264 265 switch (type) { 266 case SNES_LS_BASIC: 267 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 268 break; 269 case SNES_LS_BASIC_NONORMS: 270 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 271 break; 272 case SNES_LS_QUADRATIC: 273 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 274 break; 275 case SNES_LS_SECANT: 276 ierr = SNESLineSearchSet(snes,SNESLineSearchSecant,PETSC_NULL);CHKERRQ(ierr); 277 break; 278 default: 279 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type"); 280 break; 281 } 282 snes->ls_type = type; 283 PetscFunctionReturn(0); 284 } 285 EXTERN_C_END 286 287 288 /* -------------------------------------------------------------------------- */ 289 /*MC 290 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 291 292 Options Database: 293 294 + -snes_qn_m - Number of past states saved for the L-Broyden methods. 295 - -snes_qn_monitor - Monitors the quasi-newton jacobian. 296 297 Notes: This implements the L-BFGS algorithm for the solution of F(x) = b using previous change in F(x) and x to 298 form the approximate inverse Jacobian using a series of multiplicative rank-one updates. This will eventually be 299 generalized to implement several limited-memory Broyden methods. 300 301 References: 302 303 L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed, 304 International Journal of Computer Mathematics, vol. 86, 2009. 305 306 307 Level: beginner 308 309 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR 310 311 M*/ 312 EXTERN_C_BEGIN 313 #undef __FUNCT__ 314 #define __FUNCT__ "SNESCreate_QN" 315 PetscErrorCode SNESCreate_QN(SNES snes) 316 { 317 318 PetscErrorCode ierr; 319 QNContext * qn; 320 321 PetscFunctionBegin; 322 snes->ops->setup = SNESSetUp_QN; 323 snes->ops->solve = SNESSolve_QN; 324 snes->ops->destroy = SNESDestroy_QN; 325 snes->ops->setfromoptions = SNESSetFromOptions_QN; 326 snes->ops->view = 0; 327 snes->ops->reset = SNESReset_QN; 328 329 snes->usespc = PETSC_TRUE; 330 snes->usesksp = PETSC_FALSE; 331 332 ierr = PetscNewLog(snes, QNContext, &qn);CHKERRQ(ierr); 333 snes->data = (void *) qn; 334 qn->m = 10; 335 qn->dX = PETSC_NULL; 336 qn->dF = PETSC_NULL; 337 qn->monitor = PETSC_NULL; 338 339 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_QN",SNESLineSearchSetType_QN);CHKERRQ(ierr); 340 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 341 342 PetscFunctionReturn(0); 343 } 344 EXTERN_C_END 345