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