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