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