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