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