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 = KSPSolve(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 = KSPSolve(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 = KSPSolve(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 = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 327 snes->iter = 0; 328 snes->norm = 0.; 329 ierr = PetscObjectAMSGrantAccess((PetscObject)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 = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 350 snes->norm = fnorm; 351 ierr = PetscObjectAMSGrantAccess((PetscObject)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 = PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,B,0);CHKERRQ(ierr); 440 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 441 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,B,0);CHKERRQ(ierr); 442 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 443 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 444 snes->reason = SNES_DIVERGED_INNER; 445 PetscFunctionReturn(0); 446 } 447 ierr = SNESGetFunction(snes->pc, &FPC, NULL, NULL);CHKERRQ(ierr); 448 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 449 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 450 ierr = VecCopy(F, D);CHKERRQ(ierr); 451 } else { 452 ierr = VecCopy(F, D);CHKERRQ(ierr); 453 } 454 455 powell = PETSC_FALSE; 456 if (qn->restart_type == SNES_QN_RESTART_POWELL) { 457 /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */ 458 ierr = VecDotBegin(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 459 ierr = VecDotBegin(Dold, D, &DolddotD);CHKERRQ(ierr); 460 ierr = VecDotBegin(D, D, &DdotD);CHKERRQ(ierr); 461 ierr = VecDotBegin(Y, D, &YdotD);CHKERRQ(ierr); 462 ierr = VecDotEnd(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 463 ierr = VecDotEnd(Dold, D, &DolddotD);CHKERRQ(ierr); 464 ierr = VecDotEnd(D, D, &DdotD);CHKERRQ(ierr); 465 ierr = VecDotEnd(Y, D, &YdotD);CHKERRQ(ierr); 466 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE; 467 } 468 periodic = PETSC_FALSE; 469 if (qn->restart_type == SNES_QN_RESTART_PERIODIC) { 470 if (i_r>qn->m-1) periodic = PETSC_TRUE; 471 } 472 /* restart if either powell or periodic restart is satisfied. */ 473 if (powell || periodic) { 474 if (qn->monitor) { 475 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 476 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); 477 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 478 } 479 i_r = -1; 480 /* general purpose update */ 481 if (snes->ops->update) { 482 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 483 } 484 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 485 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 486 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 487 } 488 } 489 /* general purpose update */ 490 if (snes->ops->update) { 491 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 492 } 493 } 494 if (i == snes->max_its) { 495 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 496 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 497 } 498 PetscFunctionReturn(0); 499 } 500 501 #undef __FUNCT__ 502 #define __FUNCT__ "SNESSetUp_QN" 503 static PetscErrorCode SNESSetUp_QN(SNES snes) 504 { 505 SNES_QN *qn = (SNES_QN*)snes->data; 506 PetscErrorCode ierr; 507 508 PetscFunctionBegin; 509 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);CHKERRQ(ierr); 510 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);CHKERRQ(ierr); 511 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, 512 qn->m, PetscScalar, &qn->beta, 513 qn->m, PetscScalar, &qn->dXtdF);CHKERRQ(ierr); 514 515 if (qn->singlereduction) { 516 ierr = PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat, 517 qn->m, PetscScalar, &qn->dFtdX, 518 qn->m, PetscScalar, &qn->YtdX);CHKERRQ(ierr); 519 } 520 ierr = SNESSetWorkVecs(snes,4);CHKERRQ(ierr); 521 522 /* set up the line search */ 523 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 524 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 525 } 526 PetscFunctionReturn(0); 527 } 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "SNESReset_QN" 531 static PetscErrorCode SNESReset_QN(SNES snes) 532 { 533 PetscErrorCode ierr; 534 SNES_QN *qn; 535 536 PetscFunctionBegin; 537 if (snes->data) { 538 qn = (SNES_QN*)snes->data; 539 if (qn->U) { 540 ierr = VecDestroyVecs(qn->m, &qn->U);CHKERRQ(ierr); 541 } 542 if (qn->V) { 543 ierr = VecDestroyVecs(qn->m, &qn->V);CHKERRQ(ierr); 544 } 545 if (qn->singlereduction) { 546 ierr = PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);CHKERRQ(ierr); 547 } 548 ierr = PetscFree3(qn->alpha, qn->beta, qn->dXtdF);CHKERRQ(ierr); 549 } 550 PetscFunctionReturn(0); 551 } 552 553 #undef __FUNCT__ 554 #define __FUNCT__ "SNESDestroy_QN" 555 static PetscErrorCode SNESDestroy_QN(SNES snes) 556 { 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 561 ierr = PetscFree(snes->data);CHKERRQ(ierr); 562 ierr = PetscObjectComposeFunction((PetscObject)snes,"",NULL);CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "SNESSetFromOptions_QN" 568 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 569 { 570 571 PetscErrorCode ierr; 572 SNES_QN *qn = (SNES_QN*)snes->data; 573 PetscBool monflg = PETSC_FALSE,flg; 574 SNESLineSearch linesearch; 575 SNESQNRestartType rtype = qn->restart_type; 576 SNESQNScaleType stype = qn->scale_type; 577 578 PetscFunctionBegin; 579 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 580 ierr = PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);CHKERRQ(ierr); 581 ierr = PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);CHKERRQ(ierr); 582 ierr = PetscOptionsReal("-snes_qn_powell_downhill","Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, NULL);CHKERRQ(ierr); 583 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, NULL);CHKERRQ(ierr); 584 ierr = PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);CHKERRQ(ierr); 585 ierr = PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);CHKERRQ(ierr); 586 if (flg) ierr = SNESQNSetScaleType(snes,stype);CHKERRQ(ierr); 587 588 ierr = PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);CHKERRQ(ierr); 589 if (flg) ierr = SNESQNSetRestartType(snes,rtype);CHKERRQ(ierr); 590 591 ierr = PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes, 592 (PetscEnum)qn->type,(PetscEnum*)&qn->type,NULL);CHKERRQ(ierr); 593 ierr = PetscOptionsTail();CHKERRQ(ierr); 594 if (!snes->linesearch) { 595 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 596 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr); 597 } 598 if (monflg) { 599 qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 600 } 601 PetscFunctionReturn(0); 602 } 603 604 605 #undef __FUNCT__ 606 #define __FUNCT__ "SNESQNSetRestartType" 607 /*@ 608 SNESQNSetRestartType - Sets the restart type for SNESQN. 609 610 Logically Collective on SNES 611 612 Input Parameters: 613 + snes - the iterative context 614 - rtype - restart type 615 616 Options Database: 617 + -snes_qn_restart_type<powell,periodic,none> - set the restart type 618 - -snes_qn_restart[10] - sets the number of iterations before restart for periodic 619 620 Level: intermediate 621 622 SNESQNRestartTypes: 623 + SNES_QN_RESTART_NONE - never restart 624 . SNES_QN_RESTART_POWELL - restart based upon descent criteria 625 - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations 626 627 Notes: 628 The default line search used is the L2 line search and it requires two additional function evaluations. 629 630 .keywords: SNES, SNESQN, restart, type, set SNESLineSearch 631 @*/ 632 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) 633 { 634 PetscErrorCode ierr; 635 636 PetscFunctionBegin; 637 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 638 ierr = PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));CHKERRQ(ierr); 639 PetscFunctionReturn(0); 640 } 641 642 #undef __FUNCT__ 643 #define __FUNCT__ "SNESQNSetScaleType" 644 /*@ 645 SNESQNSetScaleType - Sets the scaling type for the inner inverse jacobian in SNESQN. 646 647 Logically Collective on SNES 648 649 Input Parameters: 650 + snes - the iterative context 651 - stype - scale type 652 653 Options Database: 654 . -snes_qn_scale_type<shanno,none,linesearch,jacobian> 655 656 Level: intermediate 657 658 SNESQNSelectTypes: 659 + SNES_QN_SCALE_NONE - don't scale the problem 660 . SNES_QN_SCALE_SHANNO - use shanno scaling 661 . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda 662 - SNES_QN_SCALE_JACOBIAN - scale by inverting a previously computed Jacobian. 663 664 .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch 665 @*/ 666 667 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) 668 { 669 PetscErrorCode ierr; 670 671 PetscFunctionBegin; 672 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 673 ierr = PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));CHKERRQ(ierr); 674 PetscFunctionReturn(0); 675 } 676 677 #undef __FUNCT__ 678 #define __FUNCT__ "SNESQNSetScaleType_QN" 679 PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) 680 { 681 SNES_QN *qn = (SNES_QN*)snes->data; 682 683 PetscFunctionBegin; 684 qn->scale_type = stype; 685 PetscFunctionReturn(0); 686 } 687 688 #undef __FUNCT__ 689 #define __FUNCT__ "SNESQNSetRestartType_QN" 690 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) 691 { 692 SNES_QN *qn = (SNES_QN*)snes->data; 693 694 PetscFunctionBegin; 695 qn->restart_type = rtype; 696 PetscFunctionReturn(0); 697 } 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 #undef __FUNCT__ 735 #define __FUNCT__ "SNESCreate_QN" 736 PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes) 737 { 738 PetscErrorCode ierr; 739 SNES_QN *qn; 740 741 PetscFunctionBegin; 742 snes->ops->setup = SNESSetUp_QN; 743 snes->ops->solve = SNESSolve_QN; 744 snes->ops->destroy = SNESDestroy_QN; 745 snes->ops->setfromoptions = SNESSetFromOptions_QN; 746 snes->ops->view = 0; 747 snes->ops->reset = SNESReset_QN; 748 749 snes->usespc = PETSC_TRUE; 750 snes->usesksp = PETSC_FALSE; 751 752 if (!snes->tolerancesset) { 753 snes->max_funcs = 30000; 754 snes->max_its = 10000; 755 } 756 757 ierr = PetscNewLog(snes,SNES_QN,&qn);CHKERRQ(ierr); 758 snes->data = (void*) qn; 759 qn->m = 10; 760 qn->scaling = 1.0; 761 qn->U = NULL; 762 qn->V = NULL; 763 qn->dXtdF = NULL; 764 qn->dFtdX = NULL; 765 qn->dXdFmat = NULL; 766 qn->monitor = NULL; 767 qn->singlereduction = PETSC_FALSE; 768 qn->powell_gamma = 0.9999; 769 qn->scale_type = SNES_QN_SCALE_SHANNO; 770 qn->restart_type = SNES_QN_RESTART_POWELL; 771 qn->type = SNES_QN_LBFGS; 772 773 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);CHKERRQ(ierr); 774 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);CHKERRQ(ierr); 775 PetscFunctionReturn(0); 776 } 777