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 = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,G,W,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 138 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); 139 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 140 if (snes->domainerror) { 141 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 142 PetscFunctionReturn(0); 143 } 144 if (!lssucceed) { 145 if (++snes->numFailures >= snes->maxFailures) { 146 snes->reason = SNES_DIVERGED_LINE_SEARCH; 147 break; 148 } 149 } 150 /* Update function and solution vectors */ 151 fnorm = gnorm; 152 ierr = VecCopy(G,F);CHKERRQ(ierr); 153 ierr = VecCopy(W,X);CHKERRQ(ierr); 154 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 155 snes->iter = i+1; 156 snes->norm = fnorm; 157 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 158 SNESLogConvHistory(snes,snes->norm,snes->iter); 159 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 160 /* set parameter for default relative tolerance convergence test */ 161 ierr = (*snes->ops->converged)(snes,i+1,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 162 if (snes->reason) PetscFunctionReturn(0); 163 164 /* set the differences */ 165 k = i % m; 166 ierr = VecCopy(F, dF[k]);CHKERRQ(ierr); 167 ierr = VecAXPY(dF[k], -1.0, Fold);CHKERRQ(ierr); 168 ierr = VecCopy(X, dX[k]);CHKERRQ(ierr); 169 ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr); 170 ierr = VecDot(dX[k], dF[k], &rhosc);CHKERRQ(ierr); 171 rho[k] = 1. / rhosc; 172 173 /* general purpose update */ 174 if (snes->ops->update) { 175 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 176 } 177 /* apply the current iteration of the approximate jacobian in order to get the next search direction*/ 178 ierr = LBGFSApplyJinv_Private(snes, i, F, Y);CHKERRQ(ierr); 179 } 180 if (i == snes->max_its) { 181 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 182 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 183 } 184 PetscFunctionReturn(0); 185 } 186 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "SNESSetUp_QN" 190 static PetscErrorCode SNESSetUp_QN(SNES snes) 191 { 192 QNContext * qn = (QNContext *)snes->data; 193 PetscErrorCode ierr; 194 PetscFunctionBegin; 195 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dX);CHKERRQ(ierr); 196 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dF);CHKERRQ(ierr); 197 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, qn->m, PetscScalar, &qn->beta, qn->m, PetscScalar, &qn->rho);CHKERRQ(ierr); 198 ierr = SNESDefaultGetWork(snes,4);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "SNESReset_QN" 204 static PetscErrorCode SNESReset_QN(SNES snes) 205 { 206 PetscErrorCode ierr; 207 QNContext * qn; 208 PetscFunctionBegin; 209 if (snes->data) { 210 qn = (QNContext *)snes->data; 211 if (qn->dX) { 212 ierr = VecDestroyVecs(qn->m, &qn->dX);CHKERRQ(ierr); 213 } 214 if (qn->dF) { 215 ierr = VecDestroyVecs(qn->m, &qn->dF);CHKERRQ(ierr); 216 } 217 ierr = PetscFree3(qn->alpha, qn->beta, qn->rho);CHKERRQ(ierr); 218 } 219 if (snes->work) {ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr);} 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "SNESDestroy_QN" 225 static PetscErrorCode SNESDestroy_QN(SNES snes) 226 { 227 PetscErrorCode ierr; 228 PetscFunctionBegin; 229 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 230 ierr = PetscFree(snes->data);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "SNESSetFromOptions_QN" 236 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 237 { 238 239 PetscErrorCode ierr; 240 QNContext * qn; 241 PetscBool monflg = PETSC_FALSE; 242 PetscFunctionBegin; 243 244 qn = (QNContext *)snes->data; 245 246 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 247 ierr = PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-Broyden methods", "SNES", qn->m, &qn->m, PETSC_NULL);CHKERRQ(ierr); 248 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNES", monflg, &monflg, PETSC_NULL);CHKERRQ(ierr); 249 ierr = PetscOptionsTail();CHKERRQ(ierr); 250 if (monflg) { 251 qn->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 252 } 253 PetscFunctionReturn(0); 254 } 255 256 EXTERN_C_BEGIN 257 #undef __FUNCT__ 258 #define __FUNCT__ "SNESLineSearchSetType_QN" 259 PetscErrorCode SNESLineSearchSetType_QN(SNES snes, SNESLineSearchType type) 260 { 261 PetscErrorCode ierr; 262 PetscFunctionBegin; 263 264 switch (type) { 265 case SNES_LS_BASIC: 266 ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr); 267 break; 268 case SNES_LS_BASIC_NONORMS: 269 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr); 270 break; 271 case SNES_LS_QUADRATIC: 272 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr); 273 break; 274 case SNES_LS_SECANT: 275 ierr = SNESLineSearchSet(snes,SNESLineSearchSecant,PETSC_NULL);CHKERRQ(ierr); 276 break; 277 default: 278 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type"); 279 break; 280 } 281 snes->ls_type = type; 282 PetscFunctionReturn(0); 283 } 284 EXTERN_C_END 285 286 287 /* -------------------------------------------------------------------------- */ 288 /*MC 289 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 290 291 Options Database: 292 293 + -snes_qn_m - Number of past states saved for the L-Broyden methods. 294 - -snes_qn_monitor - Monitors the quasi-newton jacobian. 295 296 Notes: This implements the L-BFGS algorithm for the solution of F(x) = b using previous change in F(x) and x to 297 form the approximate inverse Jacobian using a series of multiplicative rank-one updates. This will eventually be 298 generalized to implement several limited-memory Broyden methods. 299 300 References: 301 302 L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed, 303 International Journal of Computer Mathematics, vol. 86, 2009. 304 305 306 Level: beginner 307 308 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR 309 310 M*/ 311 EXTERN_C_BEGIN 312 #undef __FUNCT__ 313 #define __FUNCT__ "SNESCreate_QN" 314 PetscErrorCode SNESCreate_QN(SNES snes) 315 { 316 317 PetscErrorCode ierr; 318 QNContext * qn; 319 320 PetscFunctionBegin; 321 snes->ops->setup = SNESSetUp_QN; 322 snes->ops->solve = SNESSolve_QN; 323 snes->ops->destroy = SNESDestroy_QN; 324 snes->ops->setfromoptions = SNESSetFromOptions_QN; 325 snes->ops->view = 0; 326 snes->ops->reset = SNESReset_QN; 327 328 snes->usespc = PETSC_TRUE; 329 snes->usesksp = PETSC_FALSE; 330 331 ierr = PetscNewLog(snes, QNContext, &qn);CHKERRQ(ierr); 332 snes->data = (void *) qn; 333 qn->m = 10; 334 qn->dX = PETSC_NULL; 335 qn->dF = PETSC_NULL; 336 qn->monitor = PETSC_NULL; 337 338 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_QN",SNESLineSearchSetType_QN);CHKERRQ(ierr); 339 ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr); 340 341 PetscFunctionReturn(0); 342 } 343 EXTERN_C_END 344