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