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