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