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