1 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 2 3 #define H(i,j) qn->dXdFmat[i*qn->m + j] 4 5 const char *const SNESQNScaleTypes[] = {"NONE","SHANNO","LINESEARCH","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",0}; 6 const char *const SNESQNRestartTypes[] = {"NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0}; 7 const char *const SNESQNTypes[] = {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_RESTART_",0}; 8 9 typedef enum {SNES_QN_LBFGS = 0, 10 SNES_QN_BROYDEN = 1, 11 SNES_QN_BADBROYDEN = 2 12 } SNESQNType; 13 14 typedef struct { 15 Vec *U; /* Stored past states (vary from method to method) */ 16 Vec *V; /* Stored past states (vary from method to method) */ 17 PetscInt m; /* The number of kept previous steps */ 18 PetscScalar *alpha, *beta; 19 PetscScalar *dXtdF, *dFtdX, *YtdX; 20 PetscBool singlereduction; /* Aggregated reduction implementation */ 21 PetscScalar *dXdFmat; /* A matrix of values for dX_i dot dF_j */ 22 PetscViewer monitor; 23 PetscReal powell_gamma; /* Powell angle restart condition */ 24 PetscReal powell_downhill; /* Powell descent restart condition */ 25 PetscReal scaling; /* scaling of H0 */ 26 SNESQNType type; /* the type of quasi-newton method used */ 27 SNESQNScaleType scale_type; /* the type of scaling used */ 28 SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */ 29 } SNES_QN; 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "SNESQNApply_Broyden" 33 PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold, Vec D, Vec Dold) 34 { 35 PetscErrorCode ierr; 36 SNES_QN *qn = (SNES_QN*)snes->data; 37 Vec W = snes->work[3]; 38 Vec *U = qn->U; 39 Vec *V = qn->V; 40 KSPConvergedReason kspreason; 41 PetscInt k,i,lits; 42 PetscInt m = qn->m; 43 PetscScalar gdot; 44 PetscInt l = m; 45 46 PetscFunctionBegin; 47 if (it < m) l = it; 48 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 49 ierr = SNES_KSPSolve(snes,snes->ksp,D,W);CHKERRQ(ierr); 50 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 51 if (kspreason < 0) { 52 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 53 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 54 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 55 PetscFunctionReturn(0); 56 } 57 } 58 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 59 snes->linear_its += lits; 60 ierr = VecCopy(W,Y);CHKERRQ(ierr); 61 } else { 62 ierr = VecCopy(D,Y);CHKERRQ(ierr); 63 ierr = VecScale(Y,qn->scaling);CHKERRQ(ierr); 64 } 65 66 /* inward recursion starting at the first update and working forward */ 67 if (it > 1) { 68 for (i = 0; i < l-1; i++) { 69 k = (it+i-l)%l; 70 ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr); 71 ierr = VecAXPY(Y,gdot,V[k]);CHKERRQ(ierr); 72 if (qn->monitor) { 73 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 74 ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr); 75 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 76 } 77 } 78 } 79 if (it < m) l = it; 80 81 /* set W to be the last step, post-linesearch */ 82 ierr = VecCopy(Xold,W);CHKERRQ(ierr); 83 ierr = VecAXPY(W,-1.0,X);CHKERRQ(ierr); 84 85 if (l > 0) { 86 k = (it-1)%l; 87 ierr = VecCopy(W,U[k]);CHKERRQ(ierr); 88 ierr = VecAXPY(W,-1.0,Y);CHKERRQ(ierr); 89 ierr = VecDot(U[k],W,&gdot);CHKERRQ(ierr); 90 ierr = VecCopy(Y,V[k]);CHKERRQ(ierr); 91 ierr = VecScale(V[k],1.0/gdot);CHKERRQ(ierr); 92 ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr); 93 ierr = VecAXPY(Y,gdot,V[k]);CHKERRQ(ierr); 94 if (qn->monitor) { 95 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 96 ierr = PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr); 97 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 98 } 99 } 100 PetscFunctionReturn(0); 101 } 102 103 #undef __FUNCT__ 104 #define __FUNCT__ "SNESQNApply_BadBroyden" 105 PetscErrorCode SNESQNApply_BadBroyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold) 106 { 107 PetscErrorCode ierr; 108 SNES_QN *qn = (SNES_QN*)snes->data; 109 Vec W = snes->work[3]; 110 Vec *U = qn->U; 111 Vec *T = qn->V; 112 113 /* ksp thing for jacobian scaling */ 114 KSPConvergedReason kspreason; 115 PetscInt k, i, lits; 116 PetscInt m = qn->m; 117 PetscScalar gdot; 118 PetscInt l = m; 119 120 PetscFunctionBegin; 121 if (it < m) l = it; 122 ierr = VecCopy(D,Y);CHKERRQ(ierr); 123 if (l > 0) { 124 k = (it-1)%l; 125 ierr = VecCopy(Dold,U[k]);CHKERRQ(ierr); 126 ierr = VecAXPY(U[k],-1.0,D);CHKERRQ(ierr); 127 ierr = VecDot(U[k],U[k],&gdot);CHKERRQ(ierr); 128 ierr = VecCopy(D,T[k]);CHKERRQ(ierr); 129 ierr = VecScale(T[k],1.0/gdot);CHKERRQ(ierr); 130 ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr); 131 ierr = VecAXPY(Y,gdot,T[k]);CHKERRQ(ierr); 132 if (qn->monitor) { 133 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr); 135 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 136 } 137 } 138 139 /* inward recursion starting at the first update and working forward */ 140 if (it > 1) { 141 for (i = 0; i < l-1; i++) { 142 k = (it+i-l)%l; 143 ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr); 144 ierr = VecAXPY(Y,gdot,T[k]);CHKERRQ(ierr); 145 if (qn->monitor) { 146 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 147 ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr); 148 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 149 } 150 } 151 } 152 153 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 154 ierr = SNES_KSPSolve(snes,snes->ksp,Y,W);CHKERRQ(ierr); 155 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 156 if (kspreason < 0) { 157 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 158 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 159 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 160 PetscFunctionReturn(0); 161 } 162 } 163 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 164 snes->linear_its += lits; 165 ierr = VecCopy(W,Y);CHKERRQ(ierr); 166 } else { 167 ierr = VecScale(Y,qn->scaling);CHKERRQ(ierr); 168 } 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "SNESQNApply_LBFGS" 174 PetscErrorCode SNESQNApply_LBFGS(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold) 175 { 176 PetscErrorCode ierr; 177 SNES_QN *qn = (SNES_QN*)snes->data; 178 Vec W = snes->work[3]; 179 Vec *dX = qn->U; 180 Vec *dF = qn->V; 181 PetscScalar *alpha = qn->alpha; 182 PetscScalar *beta = qn->beta; 183 PetscScalar *dXtdF = qn->dXtdF; 184 PetscScalar *dFtdX = qn->dFtdX; 185 PetscScalar *YtdX = qn->YtdX; 186 187 /* ksp thing for jacobian scaling */ 188 KSPConvergedReason kspreason; 189 PetscInt k,i,j,g,lits; 190 PetscInt m = qn->m; 191 PetscScalar t; 192 PetscInt l = m; 193 194 PetscFunctionBegin; 195 if (it < m) l = it; 196 if (it > 0) { 197 k = (it - 1) % l; 198 ierr = VecCopy(D, dF[k]);CHKERRQ(ierr); 199 ierr = VecAXPY(dF[k], -1.0, Dold);CHKERRQ(ierr); 200 ierr = VecCopy(X, dX[k]);CHKERRQ(ierr); 201 ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr); 202 if (qn->singlereduction) { 203 ierr = VecMDot(dF[k], l, dX, dXtdF);CHKERRQ(ierr); 204 ierr = VecMDot(dX[k], l, dF, dFtdX);CHKERRQ(ierr); 205 for (j = 0; j < l; j++) { 206 H(k, j) = dFtdX[j]; 207 H(j, k) = dXtdF[j]; 208 } 209 /* copy back over to make the computation of alpha and beta easier */ 210 for (j = 0; j < l; j++) dXtdF[j] = H(j, j); 211 } else { 212 ierr = VecDot(dX[k], dF[k], &dXtdF[k]);CHKERRQ(ierr); 213 } 214 if (qn->scale_type == SNES_QN_SCALE_SHANNO) { 215 PetscReal dFtdF; 216 ierr = VecDotRealPart(dF[k],dF[k],&dFtdF);CHKERRQ(ierr); 217 qn->scaling = PetscRealPart(dXtdF[k])/dFtdF; 218 } else if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) { 219 ierr = SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);CHKERRQ(ierr); 220 } 221 } 222 223 ierr = VecCopy(D,Y);CHKERRQ(ierr); 224 225 if (qn->singlereduction) { 226 ierr = VecMDot(Y,l,dX,YtdX);CHKERRQ(ierr); 227 } 228 /* outward recursion starting at iteration k's update and working back */ 229 for (i=0; i<l; i++) { 230 k = (it-i-1)%l; 231 if (qn->singlereduction) { 232 /* construct t = dX[k] dot Y as Y_0 dot dX[k] + sum(-alpha[j]dX[k]dF[j]) */ 233 t = YtdX[k]; 234 for (j=0; j<i; j++) { 235 g = (it-j-1)%l; 236 t += -alpha[g]*H(g, k); 237 } 238 alpha[k] = t/H(k,k); 239 } else { 240 ierr = VecDot(dX[k],Y,&t);CHKERRQ(ierr); 241 alpha[k] = t/dXtdF[k]; 242 } 243 if (qn->monitor) { 244 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha: %14.12e\n", it, k, PetscRealPart(alpha[k]));CHKERRQ(ierr); 246 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 247 } 248 ierr = VecAXPY(Y,-alpha[k],dF[k]);CHKERRQ(ierr); 249 } 250 251 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 252 ierr = SNES_KSPSolve(snes,snes->ksp,Y,W);CHKERRQ(ierr); 253 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 254 if (kspreason < 0) { 255 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 256 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 257 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 258 PetscFunctionReturn(0); 259 } 260 } 261 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 262 snes->linear_its += lits; 263 ierr = VecCopy(W, Y);CHKERRQ(ierr); 264 } else { 265 ierr = VecScale(Y, qn->scaling);CHKERRQ(ierr); 266 } 267 if (qn->singlereduction) { 268 ierr = VecMDot(Y,l,dF,YtdX);CHKERRQ(ierr); 269 } 270 /* inward recursion starting at the first update and working forward */ 271 for (i = 0; i < l; i++) { 272 k = (it + i - l) % l; 273 if (qn->singlereduction) { 274 t = YtdX[k]; 275 for (j = 0; j < i; j++) { 276 g = (it + j - l) % l; 277 t += (alpha[g] - beta[g])*H(k, g); 278 } 279 beta[k] = t / H(k, k); 280 } else { 281 ierr = VecDot(dF[k], Y, &t);CHKERRQ(ierr); 282 beta[k] = t / dXtdF[k]; 283 } 284 ierr = VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);CHKERRQ(ierr); 285 if (qn->monitor) { 286 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 287 ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha - beta: %14.12e\n", it, k, PetscRealPart(alpha[k] - beta[k]));CHKERRQ(ierr); 288 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 289 } 290 } 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "SNESSolve_QN" 296 static PetscErrorCode SNESSolve_QN(SNES snes) 297 { 298 PetscErrorCode ierr; 299 SNES_QN *qn = (SNES_QN*) snes->data; 300 Vec X,Xold; 301 Vec F,B; 302 Vec Y,FPC,D,Dold; 303 SNESConvergedReason reason; 304 PetscInt i, i_r; 305 PetscReal fnorm,xnorm,ynorm,gnorm; 306 PetscBool lssucceed,powell,periodic; 307 PetscScalar DolddotD,DolddotDold,DdotD,YdotD; 308 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 309 310 /* basically just a regular newton's method except for the application of the jacobian */ 311 312 PetscFunctionBegin; 313 F = snes->vec_func; /* residual vector */ 314 Y = snes->vec_sol_update; /* search direction generated by J^-1D*/ 315 B = snes->vec_rhs; 316 317 X = snes->vec_sol; /* solution vector */ 318 Xold = snes->work[0]; 319 320 /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */ 321 D = snes->work[1]; 322 Dold = snes->work[2]; 323 324 snes->reason = SNES_CONVERGED_ITERATING; 325 326 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 327 snes->iter = 0; 328 snes->norm = 0.; 329 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 330 if (!snes->vec_func_init_set) { 331 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 332 if (snes->domainerror) { 333 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 334 PetscFunctionReturn(0); 335 } 336 } else snes->vec_func_init_set = PETSC_FALSE; 337 338 if (!snes->norm_init_set) { 339 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 340 if (PetscIsInfOrNanReal(fnorm)) { 341 snes->reason = SNES_DIVERGED_FNORM_NAN; 342 PetscFunctionReturn(0); 343 } 344 } else { 345 fnorm = snes->norm_init; 346 snes->norm_init_set = PETSC_FALSE; 347 } 348 349 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 350 snes->norm = fnorm; 351 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 352 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 353 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 354 355 /* set parameter for default relative tolerance convergence test */ 356 snes->ttol = fnorm*snes->rtol; 357 /* test convergence */ 358 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 359 if (snes->reason) PetscFunctionReturn(0); 360 361 /* composed solve */ 362 if (snes->pc && snes->pcside == PC_RIGHT) { 363 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 364 ierr = SNESSetInitialFunctionNorm(snes->pc, fnorm);CHKERRQ(ierr); 365 366 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,B,0);CHKERRQ(ierr); 367 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 368 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,B,0);CHKERRQ(ierr); 369 370 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 371 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 372 snes->reason = SNES_DIVERGED_INNER; 373 PetscFunctionReturn(0); 374 } 375 ierr = SNESGetFunction(snes->pc, &FPC, NULL, NULL);CHKERRQ(ierr); 376 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 377 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 378 ierr = VecCopy(F, Y);CHKERRQ(ierr); 379 } else { 380 ierr = VecCopy(F, Y);CHKERRQ(ierr); 381 } 382 ierr = VecCopy(Y, D);CHKERRQ(ierr); 383 384 /* scale the initial update */ 385 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 386 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 387 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 388 } 389 390 for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) { 391 switch (qn->type) { 392 case SNES_QN_BADBROYDEN: 393 ierr = SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr); 394 break; 395 case SNES_QN_BROYDEN: 396 ierr = SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr); 397 break; 398 case SNES_QN_LBFGS: 399 SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr); 400 break; 401 } 402 /* line search for lambda */ 403 ynorm = 1; gnorm = fnorm; 404 ierr = VecCopy(D, Dold);CHKERRQ(ierr); 405 ierr = VecCopy(X, Xold);CHKERRQ(ierr); 406 ierr = SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 407 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 408 if (snes->domainerror) { 409 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 410 PetscFunctionReturn(0); 411 } 412 ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);CHKERRQ(ierr); 413 if (!lssucceed) { 414 if (++snes->numFailures >= snes->maxFailures) { 415 snes->reason = SNES_DIVERGED_LINE_SEARCH; 416 break; 417 } 418 } 419 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 420 if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) { 421 ierr = SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);CHKERRQ(ierr); 422 } 423 424 /* convergence monitoring */ 425 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); 426 427 ierr = SNESSetIterationNumber(snes, i+1);CHKERRQ(ierr); 428 ierr = SNESSetFunctionNorm(snes, fnorm);CHKERRQ(ierr); 429 430 ierr = SNESLogConvergenceHistory(snes,snes->norm,snes->iter);CHKERRQ(ierr); 431 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 432 /* set parameter for default relative tolerance convergence test */ 433 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 434 if (snes->reason) PetscFunctionReturn(0); 435 436 if (snes->pc && snes->pcside == PC_RIGHT) { 437 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 438 ierr = SNESSetInitialFunctionNorm(snes->pc, fnorm);CHKERRQ(ierr); 439 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 440 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 441 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 442 snes->reason = SNES_DIVERGED_INNER; 443 PetscFunctionReturn(0); 444 } 445 ierr = SNESGetFunction(snes->pc, &FPC, NULL, NULL);CHKERRQ(ierr); 446 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 447 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 448 ierr = VecCopy(F, D);CHKERRQ(ierr); 449 } else { 450 ierr = VecCopy(F, D);CHKERRQ(ierr); 451 } 452 453 powell = PETSC_FALSE; 454 if (qn->restart_type == SNES_QN_RESTART_POWELL) { 455 /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */ 456 ierr = VecDotBegin(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 457 ierr = VecDotBegin(Dold, D, &DolddotD);CHKERRQ(ierr); 458 ierr = VecDotBegin(D, D, &DdotD);CHKERRQ(ierr); 459 ierr = VecDotBegin(Y, D, &YdotD);CHKERRQ(ierr); 460 ierr = VecDotEnd(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 461 ierr = VecDotEnd(Dold, D, &DolddotD);CHKERRQ(ierr); 462 ierr = VecDotEnd(D, D, &DdotD);CHKERRQ(ierr); 463 ierr = VecDotEnd(Y, D, &YdotD);CHKERRQ(ierr); 464 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE; 465 } 466 periodic = PETSC_FALSE; 467 if (qn->restart_type == SNES_QN_RESTART_PERIODIC) { 468 if (i_r>qn->m-1) periodic = PETSC_TRUE; 469 } 470 /* restart if either powell or periodic restart is satisfied. */ 471 if (powell || periodic) { 472 if (qn->monitor) { 473 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 474 ierr = PetscViewerASCIIPrintf(qn->monitor, "restart! |%14.12e| > %4.2f*|%14.12e| or i_r = %d\n", PetscRealPart(DolddotD), qn->powell_gamma, PetscRealPart(DolddotDold), i_r);CHKERRQ(ierr); 475 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 476 } 477 i_r = -1; 478 /* general purpose update */ 479 if (snes->ops->update) { 480 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 481 } 482 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 483 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 484 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 485 } 486 } 487 /* general purpose update */ 488 if (snes->ops->update) { 489 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 490 } 491 } 492 if (i == snes->max_its) { 493 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 494 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 495 } 496 PetscFunctionReturn(0); 497 } 498 499 #undef __FUNCT__ 500 #define __FUNCT__ "SNESSetUp_QN" 501 static PetscErrorCode SNESSetUp_QN(SNES snes) 502 { 503 SNES_QN *qn = (SNES_QN*)snes->data; 504 PetscErrorCode ierr; 505 506 PetscFunctionBegin; 507 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);CHKERRQ(ierr); 508 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);CHKERRQ(ierr); 509 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, 510 qn->m, PetscScalar, &qn->beta, 511 qn->m, PetscScalar, &qn->dXtdF);CHKERRQ(ierr); 512 513 if (qn->singlereduction) { 514 ierr = PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat, 515 qn->m, PetscScalar, &qn->dFtdX, 516 qn->m, PetscScalar, &qn->YtdX);CHKERRQ(ierr); 517 } 518 ierr = SNESDefaultGetWork(snes,4);CHKERRQ(ierr); 519 520 /* set up the line search */ 521 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 522 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 523 } 524 PetscFunctionReturn(0); 525 } 526 527 #undef __FUNCT__ 528 #define __FUNCT__ "SNESReset_QN" 529 static PetscErrorCode SNESReset_QN(SNES snes) 530 { 531 PetscErrorCode ierr; 532 SNES_QN *qn; 533 534 PetscFunctionBegin; 535 if (snes->data) { 536 qn = (SNES_QN*)snes->data; 537 if (qn->U) { 538 ierr = VecDestroyVecs(qn->m, &qn->U);CHKERRQ(ierr); 539 } 540 if (qn->V) { 541 ierr = VecDestroyVecs(qn->m, &qn->V);CHKERRQ(ierr); 542 } 543 if (qn->singlereduction) { 544 ierr = PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);CHKERRQ(ierr); 545 } 546 ierr = PetscFree3(qn->alpha, qn->beta, qn->dXtdF);CHKERRQ(ierr); 547 } 548 PetscFunctionReturn(0); 549 } 550 551 #undef __FUNCT__ 552 #define __FUNCT__ "SNESDestroy_QN" 553 static PetscErrorCode SNESDestroy_QN(SNES snes) 554 { 555 PetscErrorCode ierr; 556 557 PetscFunctionBegin; 558 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 559 ierr = PetscFree(snes->data);CHKERRQ(ierr); 560 ierr = PetscObjectComposeFunction((PetscObject)snes,"","",NULL);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563 564 #undef __FUNCT__ 565 #define __FUNCT__ "SNESSetFromOptions_QN" 566 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 567 { 568 569 PetscErrorCode ierr; 570 SNES_QN *qn = (SNES_QN*)snes->data; 571 PetscBool monflg = PETSC_FALSE,flg; 572 SNESLineSearch linesearch; 573 SNESQNRestartType rtype = qn->restart_type; 574 SNESQNScaleType stype = qn->scale_type; 575 576 PetscFunctionBegin; 577 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 578 ierr = PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);CHKERRQ(ierr); 579 ierr = PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);CHKERRQ(ierr); 580 ierr = PetscOptionsReal("-snes_qn_powell_downhill","Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, NULL);CHKERRQ(ierr); 581 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, NULL);CHKERRQ(ierr); 582 ierr = PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);CHKERRQ(ierr); 583 ierr = PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);CHKERRQ(ierr); 584 if (flg) ierr = SNESQNSetScaleType(snes,stype);CHKERRQ(ierr); 585 586 ierr = PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);CHKERRQ(ierr); 587 if (flg) ierr = SNESQNSetRestartType(snes,rtype);CHKERRQ(ierr); 588 589 ierr = PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes, 590 (PetscEnum)qn->type,(PetscEnum*)&qn->type,NULL);CHKERRQ(ierr); 591 ierr = PetscOptionsTail();CHKERRQ(ierr); 592 if (!snes->linesearch) { 593 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 594 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr); 595 } 596 if (monflg) { 597 qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 598 } 599 PetscFunctionReturn(0); 600 } 601 602 603 #undef __FUNCT__ 604 #define __FUNCT__ "SNESQNSetRestartType" 605 /*@ 606 SNESQNSetRestartType - Sets the restart type for SNESQN. 607 608 Logically Collective on SNES 609 610 Input Parameters: 611 + snes - the iterative context 612 - rtype - restart type 613 614 Options Database: 615 + -snes_qn_restart_type<powell,periodic,none> - set the restart type 616 - -snes_qn_restart[10] - sets the number of iterations before restart for periodic 617 618 Level: intermediate 619 620 SNESQNRestartTypes: 621 + SNES_QN_RESTART_NONE - never restart 622 . SNES_QN_RESTART_POWELL - restart based upon descent criteria 623 - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations 624 625 Notes: 626 The default line search used is the L2 line search and it requires two additional function evaluations. 627 628 .keywords: SNES, SNESQN, restart, type, set SNESLineSearch 629 @*/ 630 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) 631 { 632 PetscErrorCode ierr; 633 634 PetscFunctionBegin; 635 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 636 ierr = PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));CHKERRQ(ierr); 637 PetscFunctionReturn(0); 638 } 639 640 #undef __FUNCT__ 641 #define __FUNCT__ "SNESQNSetScaleType" 642 /*@ 643 SNESQNSetScaleType - Sets the scaling type for the inner inverse jacobian in SNESQN. 644 645 Logically Collective on SNES 646 647 Input Parameters: 648 + snes - the iterative context 649 - stype - scale type 650 651 Options Database: 652 . -snes_qn_scale_type<shanno,none,linesearch,jacobian> 653 654 Level: intermediate 655 656 SNESQNSelectTypes: 657 + SNES_QN_SCALE_NONE - don't scale the problem 658 . SNES_QN_SCALE_SHANNO - use shanno scaling 659 . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda 660 - SNES_QN_SCALE_JACOBIAN - scale by inverting a previously computed Jacobian. 661 662 .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch 663 @*/ 664 665 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) 666 { 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 671 ierr = PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 EXTERN_C_BEGIN 676 #undef __FUNCT__ 677 #define __FUNCT__ "SNESQNSetScaleType_QN" 678 PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) 679 { 680 SNES_QN *qn = (SNES_QN*)snes->data; 681 682 PetscFunctionBegin; 683 qn->scale_type = stype; 684 PetscFunctionReturn(0); 685 } 686 687 #undef __FUNCT__ 688 #define __FUNCT__ "SNESQNSetRestartType_QN" 689 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) 690 { 691 SNES_QN *qn = (SNES_QN*)snes->data; 692 693 PetscFunctionBegin; 694 qn->restart_type = rtype; 695 PetscFunctionReturn(0); 696 } 697 EXTERN_C_END 698 699 /* -------------------------------------------------------------------------- */ 700 /*MC 701 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 702 703 Options Database: 704 705 + -snes_qn_m - Number of past states saved for the L-Broyden methods. 706 . -snes_qn_powell_angle - Angle condition for restart. 707 . -snes_qn_powell_descent - Descent condition for restart. 708 . -snes_linesearch_type <cp, l2, basic> - Type of line search. 709 - -snes_qn_monitor - Monitors the quasi-newton jacobian. 710 711 Notes: This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using 712 previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one 713 updates. 714 715 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 716 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 717 iteration as the current iteration's values when constructing the approximate jacobian. The second, composed, 718 perturbs the problem the jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 719 720 References: 721 722 L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed, 723 International Journal of Computer Mathematics, vol. 86, 2009. 724 725 Experiments with Quasi-Newton Methods in Solving Stiff ODE Systems, Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker 726 SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985. 727 728 729 Level: beginner 730 731 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR 732 733 M*/ 734 EXTERN_C_BEGIN 735 #undef __FUNCT__ 736 #define __FUNCT__ "SNESCreate_QN" 737 PetscErrorCode SNESCreate_QN(SNES snes) 738 { 739 PetscErrorCode ierr; 740 SNES_QN *qn; 741 742 PetscFunctionBegin; 743 snes->ops->setup = SNESSetUp_QN; 744 snes->ops->solve = SNESSolve_QN; 745 snes->ops->destroy = SNESDestroy_QN; 746 snes->ops->setfromoptions = SNESSetFromOptions_QN; 747 snes->ops->view = 0; 748 snes->ops->reset = SNESReset_QN; 749 750 snes->usespc = PETSC_TRUE; 751 snes->usesksp = PETSC_FALSE; 752 753 if (!snes->tolerancesset) { 754 snes->max_funcs = 30000; 755 snes->max_its = 10000; 756 } 757 758 ierr = PetscNewLog(snes,SNES_QN,&qn);CHKERRQ(ierr); 759 snes->data = (void*) qn; 760 qn->m = 10; 761 qn->scaling = 1.0; 762 qn->U = NULL; 763 qn->V = NULL; 764 qn->dXtdF = NULL; 765 qn->dFtdX = NULL; 766 qn->dXdFmat = NULL; 767 qn->monitor = NULL; 768 qn->singlereduction = PETSC_FALSE; 769 qn->powell_gamma = 0.9999; 770 qn->scale_type = SNES_QN_SCALE_SHANNO; 771 qn->restart_type = SNES_QN_RESTART_POWELL; 772 qn->type = SNES_QN_LBFGS; 773 774 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C","SNESQNSetScaleType_QN",SNESQNSetScaleType_QN);CHKERRQ(ierr); 775 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C","SNESQNSetRestartType_QN",SNESQNSetRestartType_QN);CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 EXTERN_C_END 780