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