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 PetscScalar scl; 473 ierr = SNESLineSearchGetLambda(snes->linesearch, &scl);CHKERRQ(ierr); 474 qn->scaling = PetscRealPart(scl);CHKERRQ(ierr); 475 } 476 477 /* convergence monitoring */ 478 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); 479 480 ierr = SNESSetIterationNumber(snes, i+1);CHKERRQ(ierr); 481 ierr = SNESSetFunctionNorm(snes, fnorm);CHKERRQ(ierr); 482 483 SNESLogConvHistory(snes,snes->norm,snes->iter); 484 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 485 /* set parameter for default relative tolerance convergence test */ 486 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 487 if (snes->reason) PetscFunctionReturn(0); 488 489 490 if (snes->pc) { 491 if (qn->composition_type == SNES_QN_SEQUENTIAL) { 492 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 493 ierr = SNESSetInitialFunctionNorm(snes->pc, fnorm);CHKERRQ(ierr); 494 ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr); 495 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 496 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 497 snes->reason = SNES_DIVERGED_INNER; 498 PetscFunctionReturn(0); 499 } 500 ierr = SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 501 ierr = VecCopy(FPC, F);CHKERRQ(ierr); 502 ierr = SNESGetFunctionNorm(snes->pc, &fnorm);CHKERRQ(ierr); 503 ierr = VecCopy(F, D);CHKERRQ(ierr); 504 } else { 505 ierr = VecCopy(X, D);CHKERRQ(ierr); 506 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 507 ierr = SNESSetInitialFunctionNorm(snes->pc, fnorm);CHKERRQ(ierr); 508 ierr = SNESSolve(snes->pc, B, D);CHKERRQ(ierr); 509 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 510 if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) { 511 snes->reason = SNES_DIVERGED_INNER; 512 PetscFunctionReturn(0); 513 } 514 ierr = VecAYPX(D,-1.0,X);CHKERRQ(ierr); 515 } 516 } else { 517 ierr = VecCopy(F, D);CHKERRQ(ierr); 518 } 519 520 powell = PETSC_FALSE; 521 if (qn->restart_type == SNES_QN_RESTART_POWELL) { 522 /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */ 523 ierr = VecDotBegin(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 524 ierr = VecDotBegin(Dold, D, &DolddotD);CHKERRQ(ierr); 525 ierr = VecDotBegin(D, D, &DdotD);CHKERRQ(ierr); 526 ierr = VecDotBegin(Y, D, &YdotD);CHKERRQ(ierr); 527 ierr = VecDotEnd(Dold, Dold, &DolddotDold);CHKERRQ(ierr); 528 ierr = VecDotEnd(Dold, D, &DolddotD);CHKERRQ(ierr); 529 ierr = VecDotEnd(D, D, &DdotD);CHKERRQ(ierr); 530 ierr = VecDotEnd(Y, D, &YdotD);CHKERRQ(ierr); 531 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE; 532 } 533 periodic = PETSC_FALSE; 534 if (qn->restart_type == SNES_QN_RESTART_PERIODIC) { 535 if (i_r>qn->m-1) periodic = PETSC_TRUE; 536 } 537 /* restart if either powell or periodic restart is satisfied. */ 538 if (powell || periodic) { 539 if (qn->monitor) { 540 ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 541 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); 542 ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr); 543 } 544 i_r = -1; 545 /* general purpose update */ 546 if (snes->ops->update) { 547 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 548 } 549 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 550 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 551 } 552 } 553 /* general purpose update */ 554 if (snes->ops->update) { 555 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 556 } 557 } 558 if (i == snes->max_its) { 559 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 560 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 561 } 562 PetscFunctionReturn(0); 563 } 564 565 566 #undef __FUNCT__ 567 #define __FUNCT__ "SNESSetUp_QN" 568 static PetscErrorCode SNESSetUp_QN(SNES snes) 569 { 570 SNES_QN *qn = (SNES_QN*)snes->data; 571 PetscErrorCode ierr; 572 573 PetscFunctionBegin; 574 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);CHKERRQ(ierr); 575 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);CHKERRQ(ierr); 576 ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha, 577 qn->m, PetscScalar, &qn->beta, 578 qn->m, PetscScalar, &qn->dXtdF);CHKERRQ(ierr); 579 580 if (qn->singlereduction) { 581 ierr = PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat, 582 qn->m, PetscScalar, &qn->dFtdX, 583 qn->m, PetscScalar, &qn->YtdX);CHKERRQ(ierr); 584 } 585 ierr = SNESDefaultGetWork(snes,4);CHKERRQ(ierr); 586 587 /* set up the line search */ 588 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) { 589 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 590 } 591 PetscFunctionReturn(0); 592 } 593 594 #undef __FUNCT__ 595 #define __FUNCT__ "SNESReset_QN" 596 static PetscErrorCode SNESReset_QN(SNES snes) 597 { 598 PetscErrorCode ierr; 599 SNES_QN *qn; 600 PetscFunctionBegin; 601 if (snes->data) { 602 qn = (SNES_QN*)snes->data; 603 if (qn->U) { 604 ierr = VecDestroyVecs(qn->m, &qn->U);CHKERRQ(ierr); 605 } 606 if (qn->V) { 607 ierr = VecDestroyVecs(qn->m, &qn->V);CHKERRQ(ierr); 608 } 609 if (qn->singlereduction) { 610 ierr = PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);CHKERRQ(ierr); 611 } 612 ierr = PetscFree3(qn->alpha, qn->beta, qn->dXtdF);CHKERRQ(ierr); 613 } 614 PetscFunctionReturn(0); 615 } 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "SNESDestroy_QN" 619 static PetscErrorCode SNESDestroy_QN(SNES snes) 620 { 621 PetscErrorCode ierr; 622 PetscFunctionBegin; 623 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 624 ierr = PetscFree(snes->data);CHKERRQ(ierr); 625 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"","",PETSC_NULL);CHKERRQ(ierr); 626 PetscFunctionReturn(0); 627 } 628 629 #undef __FUNCT__ 630 #define __FUNCT__ "SNESSetFromOptions_QN" 631 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 632 { 633 634 PetscErrorCode ierr; 635 SNES_QN *qn; 636 PetscBool monflg = PETSC_FALSE; 637 SNESLineSearch linesearch; 638 PetscFunctionBegin; 639 640 qn = (SNES_QN*)snes->data; 641 642 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 643 ierr = PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,PETSC_NULL);CHKERRQ(ierr); 644 ierr = PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, PETSC_NULL);CHKERRQ(ierr); 645 ierr = PetscOptionsReal("-snes_qn_powell_downhill","Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, PETSC_NULL);CHKERRQ(ierr); 646 ierr = PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, PETSC_NULL);CHKERRQ(ierr); 647 ierr = PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, PETSC_NULL);CHKERRQ(ierr); 648 ierr = PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)qn->scale_type,(PetscEnum*)&qn->scale_type,PETSC_NULL);CHKERRQ(ierr); 649 ierr = PetscOptionsEnum("-snes_qn_composition_type","Composition type","SNESQNSetCompositionType",SNESQNCompositionTypes, 650 (PetscEnum)qn->composition_type,(PetscEnum*)&qn->composition_type,PETSC_NULL);CHKERRQ(ierr); 651 ierr = PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes, 652 (PetscEnum)qn->restart_type,(PetscEnum*)&qn->restart_type,PETSC_NULL);CHKERRQ(ierr); 653 ierr = PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes, 654 (PetscEnum)qn->type,(PetscEnum*)&qn->type,PETSC_NULL);CHKERRQ(ierr); 655 ierr = PetscOptionsTail();CHKERRQ(ierr); 656 if (!snes->linesearch) { 657 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 658 if (!snes->pc || qn->composition_type == SNES_QN_SEQUENTIAL) { 659 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr); 660 } else { 661 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr); 662 } 663 } 664 if (monflg) { 665 qn->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr); 666 } 667 PetscFunctionReturn(0); 668 } 669 670 671 #undef __FUNCT__ 672 #define __FUNCT__ "SNESQNSetRestartType" 673 /*@ 674 SNESQNSetRestartType - Sets the restart type for SNESQN. 675 676 Logically Collective on SNES 677 678 Input Parameters: 679 + snes - the iterative context 680 - rtype - restart type 681 682 Options Database: 683 + -snes_qn_restart_type<powell,periodic,none> - set the restart type 684 - -snes_qn_restart[10] - sets the number of iterations before restart for periodic 685 686 Level: intermediate 687 688 SNESQNRestartTypes: 689 + SNES_QN_RESTART_NONE - never restart 690 . SNES_QN_RESTART_POWELL - restart based upon descent criteria 691 - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations 692 693 Notes: 694 The default line search used is the L2 line search and it requires two additional function evaluations. 695 696 .keywords: SNES, SNESQN, restart, type, set SNESLineSearch 697 @*/ 698 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) { 699 PetscErrorCode ierr; 700 PetscFunctionBegin; 701 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 702 ierr = PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));CHKERRQ(ierr); 703 PetscFunctionReturn(0); 704 } 705 706 #undef __FUNCT__ 707 #define __FUNCT__ "SNESQNSetScaleType" 708 /*@ 709 SNESQNSetScaleType - Sets the scaling type for the inner inverse jacobian in SNESQN. 710 711 Logically Collective on SNES 712 713 Input Parameters: 714 + snes - the iterative context 715 - stype - scale type 716 717 Options Database: 718 . -snes_qn_scale_type<shanno,none,linesearch,jacobian> 719 720 Level: intermediate 721 722 SNESQNSelectTypes: 723 + SNES_QN_SCALE_NONE - don't scale the problem 724 . SNES_QN_SCALE_SHANNO - use shanno scaling 725 . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda 726 - SNES_QN_SCALE_JACOBIAN - scale by inverting a previously computed Jacobian. 727 728 .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch 729 @*/ 730 731 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) { 732 PetscErrorCode ierr; 733 PetscFunctionBegin; 734 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 735 ierr = PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));CHKERRQ(ierr); 736 PetscFunctionReturn(0); 737 } 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "SNESQNSetCompositionType" 741 /*@ 742 SNESQNSetCompositionType - Sets the composition type 743 744 Logically Collective on SNES 745 746 Input Parameters: 747 + snes - the iterative context 748 - stype - composition type 749 750 Options Database: 751 . -snes_qn_composition_type<sequential, composed> 752 753 Level: intermediate 754 755 SNESQNSelectTypes: 756 + SNES_QN_COMPOSITION_SEQUENTIAL - Solve the system with X = PC(X) and D = F(PC(X)) 757 - SNES_QN_COMPOSITION_COMPOSED - solve the system with X = X and D = PC(X) - X 758 759 .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch 760 @*/ 761 762 PetscErrorCode SNESQNSetCompositionType(SNES snes, SNESQNCompositionType ctype) { 763 PetscErrorCode ierr; 764 PetscFunctionBegin; 765 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 766 ierr = PetscTryMethod(snes,"SNESQNSetCompositionType_C",(SNES,SNESQNCompositionType),(snes,ctype));CHKERRQ(ierr); 767 PetscFunctionReturn(0); 768 } 769 770 EXTERN_C_BEGIN 771 #undef __FUNCT__ 772 #define __FUNCT__ "SNESQNSetScaleType_QN" 773 PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) { 774 SNES_QN *qn = (SNES_QN *)snes->data; 775 PetscFunctionBegin; 776 qn->scale_type = stype; 777 PetscFunctionReturn(0); 778 } 779 780 #undef __FUNCT__ 781 #define __FUNCT__ "SNESQNSetRestartType_QN" 782 PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) { 783 SNES_QN *qn = (SNES_QN *)snes->data; 784 PetscFunctionBegin; 785 qn->restart_type = rtype; 786 PetscFunctionReturn(0); 787 } 788 789 #undef __FUNCT__ 790 #define __FUNCT__ "SNESQNSetCompositionType_QN" 791 792 PetscErrorCode SNESQNSetCompositionType_QN(SNES snes, SNESQNCompositionType ctype) { 793 SNES_QN *qn = (SNES_QN *)snes->data; 794 PetscFunctionBegin; 795 qn->composition_type = ctype; 796 PetscFunctionReturn(0); 797 } 798 EXTERN_C_END 799 800 /* -------------------------------------------------------------------------- */ 801 /*MC 802 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 803 804 Options Database: 805 806 + -snes_qn_m - Number of past states saved for the L-Broyden methods. 807 . -snes_qn_powell_angle - Angle condition for restart. 808 . -snes_qn_powell_descent - Descent condition for restart. 809 . -snes_qn_composition <sequential, composed>- Type of composition. 810 . -snes_linesearch_type <cp, l2, basic> - Type of line search. 811 - -snes_qn_monitor - Monitors the quasi-newton jacobian. 812 813 Notes: This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using 814 previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one 815 updates. 816 817 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of 818 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this 819 iteration as the current iteration's values when constructing the approximate jacobian. The second, composed, 820 perturbs the problem the jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner. 821 822 References: 823 824 L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed, 825 International Journal of Computer Mathematics, vol. 86, 2009. 826 827 Experiments with Quasi-Newton Methods in Solving Stiff ODE Systems, Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker 828 SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985. 829 830 831 Level: beginner 832 833 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR 834 835 M*/ 836 EXTERN_C_BEGIN 837 #undef __FUNCT__ 838 #define __FUNCT__ "SNESCreate_QN" 839 PetscErrorCode SNESCreate_QN(SNES snes) 840 { 841 842 PetscErrorCode ierr; 843 SNES_QN *qn; 844 845 PetscFunctionBegin; 846 snes->ops->setup = SNESSetUp_QN; 847 snes->ops->solve = SNESSolve_QN; 848 snes->ops->destroy = SNESDestroy_QN; 849 snes->ops->setfromoptions = SNESSetFromOptions_QN; 850 snes->ops->view = 0; 851 snes->ops->reset = SNESReset_QN; 852 853 snes->usespc = PETSC_TRUE; 854 snes->usesksp = PETSC_FALSE; 855 856 if (!snes->tolerancesset) { 857 snes->max_funcs = 30000; 858 snes->max_its = 10000; 859 } 860 861 ierr = PetscNewLog(snes,SNES_QN,&qn);CHKERRQ(ierr); 862 snes->data = (void *) qn; 863 qn->m = 10; 864 qn->scaling = 1.0; 865 qn->U = PETSC_NULL; 866 qn->V = PETSC_NULL; 867 qn->dXtdF = PETSC_NULL; 868 qn->dFtdX = PETSC_NULL; 869 qn->dXdFmat = PETSC_NULL; 870 qn->monitor = PETSC_NULL; 871 qn->singlereduction = PETSC_FALSE; 872 qn->powell_gamma = 0.9999; 873 qn->composition_type= SNES_QN_SEQUENTIAL; 874 qn->scale_type = SNES_QN_SCALE_SHANNO; 875 qn->restart_type = SNES_QN_RESTART_POWELL; 876 qn->type = SNES_QN_LBFGS; 877 878 PetscFunctionReturn(0); 879 } 880 881 EXTERN_C_END 882