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