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