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