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