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 /* scale the initial update */ 359 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 360 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 361 SNESCheckJacobianDomainerror(snes); 362 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 363 } 364 365 for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) { 366 if (qn->scale_type == SNES_QN_SCALE_SHANNO && i_r > 0) { 367 PetscScalar ff,xf; 368 ierr = VecCopy(Dold,Y);CHKERRQ(ierr); 369 ierr = VecCopy(Xold,W);CHKERRQ(ierr); 370 ierr = VecAXPY(Y,-1.0,D);CHKERRQ(ierr); 371 ierr = VecAXPY(W,-1.0,X);CHKERRQ(ierr); 372 ierr = VecDotBegin(Y,Y,&ff);CHKERRQ(ierr); 373 ierr = VecDotBegin(W,Y,&xf);CHKERRQ(ierr); 374 ierr = VecDotEnd(Y,Y,&ff);CHKERRQ(ierr); 375 ierr = VecDotEnd(W,Y,&xf);CHKERRQ(ierr); 376 qn->scaling = PetscRealPart(xf)/PetscRealPart(ff); 377 if (qn->monitor) { 378 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 379 ierr = PetscViewerASCIIPrintf(qn->monitor, "Shanno scaling %D %g\n", i,(double)qn->scaling);CHKERRQ(ierr); 380 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 381 } 382 } 383 switch (qn->type) { 384 case SNES_QN_BADBROYDEN: 385 ierr = SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr); 386 break; 387 case SNES_QN_BROYDEN: 388 ierr = SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D);CHKERRQ(ierr); 389 break; 390 case SNES_QN_LBFGS: 391 SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr); 392 break; 393 } 394 /* line search for lambda */ 395 ynorm = 1; gnorm = fnorm; 396 ierr = VecCopy(D, Dold);CHKERRQ(ierr); 397 ierr = VecCopy(X, Xold);CHKERRQ(ierr); 398 ierr = SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 399 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 400 ierr = SNESLineSearchGetReason(snes->linesearch, &lssucceed);CHKERRQ(ierr); 401 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 402 if (lssucceed) { 403 if (++snes->numFailures >= snes->maxFailures) { 404 snes->reason = SNES_DIVERGED_LINE_SEARCH; 405 break; 406 } 407 } 408 if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) { 409 ierr = SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);CHKERRQ(ierr); 410 } 411 412 /* convergence monitoring */ 413 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); 414 415 if (snes->npc && snes->npcside== PC_RIGHT) { 416 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);CHKERRQ(ierr); 417 ierr = SNESSolve(snes->npc,snes->vec_rhs,X);CHKERRQ(ierr); 418 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);CHKERRQ(ierr); 419 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 420 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 421 snes->reason = SNES_DIVERGED_INNER; 422 PetscFunctionReturn(0); 423 } 424 ierr = SNESGetNPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 425 } 426 427 ierr = SNESSetIterationNumber(snes, i+1);CHKERRQ(ierr); 428 snes->norm = fnorm; 429 snes->xnorm = xnorm; 430 snes->ynorm = ynorm; 431 432 ierr = SNESLogConvergenceHistory(snes,snes->norm,snes->iter);CHKERRQ(ierr); 433 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 434 /* set parameter for default relative tolerance convergence test */ 435 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 436 if (snes->reason) PetscFunctionReturn(0); 437 if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 438 ierr = SNESApplyNPC(snes,X,F,D);CHKERRQ(ierr); 439 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 440 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 441 snes->reason = SNES_DIVERGED_INNER; 442 PetscFunctionReturn(0); 443 } 444 } else { 445 ierr = VecCopy(F, D);CHKERRQ(ierr); 446 } 447 powell = PETSC_FALSE; 448 if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) { 449 /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */ 450 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 451 ierr = MatMult(snes->jacobian_pre,Dold,W);CHKERRQ(ierr); 452 } else { 453 ierr = VecCopy(Dold,W);CHKERRQ(ierr); 454 } 455 ierr = VecDotBegin(W, Dold, &DolddotDold);CHKERRQ(ierr); 456 ierr = VecDotBegin(W, D, &DolddotD);CHKERRQ(ierr); 457 ierr = VecDotEnd(W, Dold, &DolddotDold);CHKERRQ(ierr); 458 ierr = VecDotEnd(W, D, &DolddotD);CHKERRQ(ierr); 459 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE; 460 } 461 periodic = PETSC_FALSE; 462 if (qn->restart_type == SNES_QN_RESTART_PERIODIC) { 463 if (i_r>qn->m-1) periodic = PETSC_TRUE; 464 } 465 /* restart if either powell or periodic restart is satisfied. */ 466 if (powell || periodic) { 467 if (qn->monitor) { 468 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 469 if (powell) { 470 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); 471 } else { 472 ierr = PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %D\n", i_r);CHKERRQ(ierr); 473 } 474 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 475 } 476 i_r = -1; 477 /* general purpose update */ 478 if (snes->ops->update) { 479 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 480 } 481 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 482 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 483 SNESCheckJacobianDomainerror(snes); 484 } 485 } 486 /* general purpose update */ 487 if (snes->ops->update) { 488 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 489 } 490 } 491 if (i == snes->max_its) { 492 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 493 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 494 } 495 PetscFunctionReturn(0); 496 } 497 498 static PetscErrorCode SNESSetUp_QN(SNES snes) 499 { 500 SNES_QN *qn = (SNES_QN*)snes->data; 501 PetscErrorCode ierr; 502 DM dm; 503 504 PetscFunctionBegin; 505 506 if (!snes->vec_sol) { 507 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 508 ierr = DMCreateGlobalVector(dm,&snes->vec_sol);CHKERRQ(ierr); 509 } 510 511 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);CHKERRQ(ierr); 512 if (qn->type != SNES_QN_BROYDEN) ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);CHKERRQ(ierr); 513 ierr = PetscMalloc4(qn->m,&qn->alpha,qn->m,&qn->beta,qn->m,&qn->dXtdF,qn->m,&qn->lambda);CHKERRQ(ierr); 514 515 if (qn->singlereduction) { 516 ierr = PetscMalloc3(qn->m*qn->m,&qn->dXdFmat,qn->m,&qn->dFtdX,qn->m,&qn->YtdX);CHKERRQ(ierr); 517 } 518 ierr = SNESSetWorkVecs(snes,4);CHKERRQ(ierr); 519 /* set method defaults */ 520 if (qn->scale_type == SNES_QN_SCALE_DEFAULT) { 521 if (qn->type == SNES_QN_BADBROYDEN) { 522 qn->scale_type = SNES_QN_SCALE_NONE; 523 } else { 524 qn->scale_type = SNES_QN_SCALE_SHANNO; 525 } 526 } 527 if (qn->restart_type == SNES_QN_RESTART_DEFAULT) { 528 if (qn->type == SNES_QN_LBFGS) { 529 qn->restart_type = SNES_QN_RESTART_POWELL; 530 } else { 531 qn->restart_type = SNES_QN_RESTART_PERIODIC; 532 } 533 } 534 535 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 536 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 537 } 538 if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;} 539 PetscFunctionReturn(0); 540 } 541 542 static PetscErrorCode SNESReset_QN(SNES snes) 543 { 544 PetscErrorCode ierr; 545 SNES_QN *qn; 546 547 PetscFunctionBegin; 548 if (snes->data) { 549 qn = (SNES_QN*)snes->data; 550 if (qn->U) { 551 ierr = VecDestroyVecs(qn->m, &qn->U);CHKERRQ(ierr); 552 } 553 if (qn->V) { 554 ierr = VecDestroyVecs(qn->m, &qn->V);CHKERRQ(ierr); 555 } 556 if (qn->singlereduction) { 557 ierr = PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);CHKERRQ(ierr); 558 } 559 ierr = PetscFree4(qn->alpha,qn->beta,qn->dXtdF,qn->lambda);CHKERRQ(ierr); 560 } 561 PetscFunctionReturn(0); 562 } 563 564 static PetscErrorCode SNESDestroy_QN(SNES snes) 565 { 566 PetscErrorCode ierr; 567 568 PetscFunctionBegin; 569 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 570 ierr = PetscFree(snes->data);CHKERRQ(ierr); 571 ierr = PetscObjectComposeFunction((PetscObject)snes,"",NULL);CHKERRQ(ierr); 572 PetscFunctionReturn(0); 573 } 574 575 static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes) 576 { 577 578 PetscErrorCode ierr; 579 SNES_QN *qn = (SNES_QN*)snes->data; 580 PetscBool monflg = PETSC_FALSE,flg; 581 SNESLineSearch linesearch; 582 SNESQNRestartType rtype = qn->restart_type; 583 SNESQNScaleType stype = qn->scale_type; 584 SNESQNType qtype = qn->type; 585 586 PetscFunctionBegin; 587 ierr = PetscOptionsHead(PetscOptionsObject,"SNES QN options");CHKERRQ(ierr); 588 ierr = PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);CHKERRQ(ierr); 589 ierr = PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);CHKERRQ(ierr); 590 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, NULL);CHKERRQ(ierr); 591 ierr = PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);CHKERRQ(ierr); 592 ierr = PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);CHKERRQ(ierr); 593 if (flg) ierr = SNESQNSetScaleType(snes,stype);CHKERRQ(ierr); 594 595 ierr = PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);CHKERRQ(ierr); 596 if (flg) ierr = SNESQNSetRestartType(snes,rtype);CHKERRQ(ierr); 597 598 ierr = PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,(PetscEnum)qtype,(PetscEnum*)&qtype,&flg);CHKERRQ(ierr); 599 if (flg) {ierr = SNESQNSetType(snes,qtype);CHKERRQ(ierr);} 600 ierr = PetscOptionsTail();CHKERRQ(ierr); 601 if (!snes->linesearch) { 602 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 603 if (qn->type == SNES_QN_LBFGS) { 604 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr); 605 } else if (qn->type == SNES_QN_BROYDEN) { 606 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);CHKERRQ(ierr); 607 } else { 608 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr); 609 } 610 } 611 if (monflg) { 612 qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 613 } 614 PetscFunctionReturn(0); 615 } 616 617 static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer) 618 { 619 SNES_QN *qn = (SNES_QN*)snes->data; 620 PetscBool iascii; 621 PetscErrorCode ierr; 622 623 PetscFunctionBegin; 624 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 625 if (iascii) { 626 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); 627 ierr = PetscViewerASCIIPrintf(viewer," Stored subspace size: %D\n", qn->m);CHKERRQ(ierr); 628 if (qn->singlereduction) { 629 ierr = PetscViewerASCIIPrintf(viewer," Using the single reduction variant.\n");CHKERRQ(ierr); 630 } 631 } 632 PetscFunctionReturn(0); 633 } 634 635 /*@ 636 SNESQNSetRestartType - Sets the restart type for SNESQN. 637 638 Logically Collective on SNES 639 640 Input Parameters: 641 + snes - the iterative context 642 - rtype - restart type 643 644 Options Database: 645 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 646 - -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic 647 648 Level: intermediate 649 650 SNESQNRestartTypes: 651 + SNES_QN_RESTART_NONE - never restart 652 . SNES_QN_RESTART_POWELL - restart based upon descent criteria 653 - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations 654 655 .keywords: SNES, SNESQN, restart, type, set SNESLineSearch, SNESQNRestartType 656 @*/ 657 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) 658 { 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 663 ierr = PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 /*@ 668 SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN. 669 670 Logically Collective on SNES 671 672 Input Parameters: 673 + snes - the iterative context 674 - stype - scale type 675 676 Options Database: 677 . -snes_qn_scale_type <shanno,none,linesearch,jacobian> 678 679 Level: intermediate 680 681 SNESQNScaleTypes: 682 + SNES_QN_SCALE_NONE - don't scale the problem 683 . SNES_QN_SCALE_SHANNO - use shanno scaling 684 . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda 685 - SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration 686 of QN and at ever restart. 687 688 .keywords: scaling type 689 690 .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian() 691 @*/ 692 693 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) 694 { 695 PetscErrorCode ierr; 696 697 PetscFunctionBegin; 698 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 699 ierr = PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));CHKERRQ(ierr); 700 PetscFunctionReturn(0); 701 } 702 703 PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) 704 { 705 SNES_QN *qn = (SNES_QN*)snes->data; 706 707 PetscFunctionBegin; 708 qn->scale_type = stype; 709 PetscFunctionReturn(0); 710 } 711 712 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) 713 { 714 SNES_QN *qn = (SNES_QN*)snes->data; 715 716 PetscFunctionBegin; 717 qn->restart_type = rtype; 718 PetscFunctionReturn(0); 719 } 720 721 /*@ 722 SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN. 723 724 Logically Collective on SNES 725 726 Input Parameters: 727 + snes - the iterative context 728 - qtype - variant type 729 730 Options Database: 731 . -snes_qn_type <lbfgs,broyden,badbroyden> 732 733 Level: beginner 734 735 SNESQNTypes: 736 + SNES_QN_LBFGS - LBFGS variant 737 . SNES_QN_BROYDEN - Broyden variant 738 - SNES_QN_BADBROYDEN - Bad Broyden variant 739 740 .keywords: SNES, SNESQN, type, set, SNESQNType 741 @*/ 742 743 PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype) 744 { 745 PetscErrorCode ierr; 746 747 PetscFunctionBegin; 748 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 749 ierr = PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));CHKERRQ(ierr); 750 PetscFunctionReturn(0); 751 } 752 753 PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype) 754 { 755 SNES_QN *qn = (SNES_QN*)snes->data; 756 757 PetscFunctionBegin; 758 qn->type = qtype; 759 PetscFunctionReturn(0); 760 } 761 762 /* -------------------------------------------------------------------------- */ 763 /*MC 764 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 765 766 Options Database: 767 768 + -snes_qn_m <m> - Number of past states saved for the L-Broyden methods. 769 + -snes_qn_restart_type <powell,periodic,none> - set the restart type 770 . -snes_qn_powell_gamma - Angle condition for restart. 771 . -snes_qn_powell_descent - Descent condition for restart. 772 . -snes_qn_type <lbfgs,broyden,badbroyden> - QN type 773 . -snes_qn_scale_type <shanno,none,linesearch,jacobian> - scaling performed on inner Jacobian 774 . -snes_linesearch_type <cp, l2, basic> - Type of line search. 775 - -snes_qn_monitor - Monitors the quasi-newton Jacobian. 776 777 Notes: 778 This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using 779 previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one 780 updates. 781 782 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 783 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 784 iteration as the current iteration's values when constructing the approximate Jacobian. The second, composed, 785 perturbs the problem the Jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 786 787 Uses left nonlinear preconditioning by default. 788 789 References: 790 + 1. - Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995. 791 . 2. - R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods, 792 Technical Report, Northwestern University, June 1992. 793 . 3. - Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE 794 Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985. 795 - 4. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 796 SIAM Review, 57(4), 2015 797 798 Level: beginner 799 800 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR 801 802 M*/ 803 PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes) 804 { 805 PetscErrorCode ierr; 806 SNES_QN *qn; 807 808 PetscFunctionBegin; 809 snes->ops->setup = SNESSetUp_QN; 810 snes->ops->solve = SNESSolve_QN; 811 snes->ops->destroy = SNESDestroy_QN; 812 snes->ops->setfromoptions = SNESSetFromOptions_QN; 813 snes->ops->view = SNESView_QN; 814 snes->ops->reset = SNESReset_QN; 815 816 snes->npcside= PC_LEFT; 817 818 snes->usesnpc = PETSC_TRUE; 819 snes->usesksp = PETSC_FALSE; 820 821 snes->alwayscomputesfinalresidual = PETSC_TRUE; 822 823 if (!snes->tolerancesset) { 824 snes->max_funcs = 30000; 825 snes->max_its = 10000; 826 } 827 828 ierr = PetscNewLog(snes,&qn);CHKERRQ(ierr); 829 snes->data = (void*) qn; 830 qn->m = 10; 831 qn->scaling = 1.0; 832 qn->U = NULL; 833 qn->V = NULL; 834 qn->dXtdF = NULL; 835 qn->dFtdX = NULL; 836 qn->dXdFmat = NULL; 837 qn->monitor = NULL; 838 qn->singlereduction = PETSC_TRUE; 839 qn->powell_gamma = 0.9999; 840 qn->scale_type = SNES_QN_SCALE_DEFAULT; 841 qn->restart_type = SNES_QN_RESTART_DEFAULT; 842 qn->type = SNES_QN_LBFGS; 843 844 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);CHKERRQ(ierr); 845 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);CHKERRQ(ierr); 846 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);CHKERRQ(ierr); 847 PetscFunctionReturn(0); 848 } 849