1 2 /* 3 This file implements PGMRES (a Pipelined Generalized Minimal Residual method) 4 */ 5 6 #include <../src/ksp/ksp/impls/gmres/pgmres/pgmresimpl.h> /*I "petscksp.h" I*/ 7 #define PGMRES_DELTA_DIRECTIONS 10 8 #define PGMRES_DEFAULT_MAXK 30 9 10 static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool*,PetscReal*); 11 static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt); 12 13 /* 14 15 KSPSetUp_PGMRES - Sets up the workspace needed by pgmres. 16 17 This is called once, usually automatically by KSPSolve() or KSPSetUp(), 18 but can be called directly by KSPSetUp(). 19 20 */ 21 static PetscErrorCode KSPSetUp_PGMRES(KSP ksp) 22 { 23 PetscFunctionBegin; 24 PetscCall(KSPSetUp_GMRES(ksp)); 25 PetscFunctionReturn(0); 26 } 27 28 /* 29 30 KSPPGMRESCycle - Run pgmres, possibly with restart. Return residual 31 history if requested. 32 33 input parameters: 34 . pgmres - structure containing parameters and work areas 35 36 output parameters: 37 . itcount - number of iterations used. If null, ignored. 38 . converged - 0 if not converged 39 40 Notes: 41 On entry, the value in vector VEC_VV(0) should be 42 the initial residual. 43 44 */ 45 static PetscErrorCode KSPPGMRESCycle(PetscInt *itcount,KSP ksp) 46 { 47 KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data); 48 PetscReal res_norm,res,newnorm; 49 PetscInt it = 0,j,k; 50 PetscBool hapend = PETSC_FALSE; 51 52 PetscFunctionBegin; 53 if (itcount) *itcount = 0; 54 PetscCall(VecNormalize(VEC_VV(0),&res_norm)); 55 KSPCheckNorm(ksp,res_norm); 56 res = res_norm; 57 *RS(0) = res_norm; 58 59 /* check for the convergence */ 60 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 61 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res; 62 else ksp->rnorm = 0; 63 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 64 pgmres->it = it-2; 65 PetscCall(KSPLogResidualHistory(ksp,ksp->rnorm)); 66 PetscCall(KSPMonitor(ksp,ksp->its,ksp->rnorm)); 67 if (!res) { 68 ksp->reason = KSP_CONVERGED_ATOL; 69 PetscCall(PetscInfo(ksp,"Converged due to zero residual norm on entry\n")); 70 PetscFunctionReturn(0); 71 } 72 73 PetscCall((*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP)); 74 for (; !ksp->reason; it++) { 75 Vec Zcur,Znext; 76 if (pgmres->vv_allocated <= it + VEC_OFFSET + 1) { 77 PetscCall(KSPGMRESGetNewVectors(ksp,it+1)); 78 } 79 /* VEC_VV(it-1) is orthogonal, it will be normalized once the VecNorm arrives. */ 80 Zcur = VEC_VV(it); /* Zcur is not yet orthogonal, but the VecMDot to orthogonalize it has been started. */ 81 Znext = VEC_VV(it+1); /* This iteration will compute Znext, update with a deferred correction once we know how 82 * Zcur relates to the previous vectors, and start the reduction to orthogonalize it. */ 83 84 if (it < pgmres->max_k+1 && ksp->its+1 < PetscMax(2,ksp->max_it)) { /* We don't know whether what we have computed is enough, so apply the matrix. */ 85 PetscCall(KSP_PCApplyBAorAB(ksp,Zcur,Znext,VEC_TEMP_MATOP)); 86 } 87 88 if (it > 1) { /* Complete the pending reduction */ 89 PetscCall(VecNormEnd(VEC_VV(it-1),NORM_2,&newnorm)); 90 *HH(it-1,it-2) = newnorm; 91 } 92 if (it > 0) { /* Finish the reduction computing the latest column of H */ 93 PetscCall(VecMDotEnd(Zcur,it,&(VEC_VV(0)),HH(0,it-1))); 94 } 95 96 if (it > 1) { 97 /* normalize the base vector from two iterations ago, basis is complete up to here */ 98 PetscCall(VecScale(VEC_VV(it-1),1./ *HH(it-1,it-2))); 99 100 PetscCall(KSPPGMRESUpdateHessenberg(ksp,it-2,&hapend,&res)); 101 pgmres->it = it-2; 102 ksp->its++; 103 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res; 104 else ksp->rnorm = 0; 105 106 PetscCall((*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP)); 107 if (it < pgmres->max_k+1 || ksp->reason || ksp->its == ksp->max_it) { /* Monitor if we are done or still iterating, but not before a restart. */ 108 PetscCall(KSPLogResidualHistory(ksp,ksp->rnorm)); 109 PetscCall(KSPMonitor(ksp,ksp->its,ksp->rnorm)); 110 } 111 if (ksp->reason) break; 112 /* Catch error in happy breakdown and signal convergence and break from loop */ 113 if (hapend) { 114 PetscCheck(!ksp->errorifnotconverged,PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the happy break down, but convergence was not indicated. Residual norm = %g",(double)res); 115 else { 116 ksp->reason = KSP_DIVERGED_BREAKDOWN; 117 break; 118 } 119 } 120 121 if (!(it < pgmres->max_k+1 && ksp->its < ksp->max_it)) break; 122 123 /* The it-2 column of H was not scaled when we computed Zcur, apply correction */ 124 PetscCall(VecScale(Zcur,1./ *HH(it-1,it-2))); 125 /* And Znext computed in this iteration was computed using the under-scaled Zcur */ 126 PetscCall(VecScale(Znext,1./ *HH(it-1,it-2))); 127 128 /* In the previous iteration, we projected an unnormalized Zcur against the Krylov basis, so we need to fix the column of H resulting from that projection. */ 129 for (k=0; k<it; k++) *HH(k,it-1) /= *HH(it-1,it-2); 130 /* When Zcur was projected against the Krylov basis, VV(it-1) was still not normalized, so fix that too. This 131 * column is complete except for HH(it,it-1) which we won't know until the next iteration. */ 132 *HH(it-1,it-1) /= *HH(it-1,it-2); 133 } 134 135 if (it > 0) { 136 PetscScalar *work; 137 if (!pgmres->orthogwork) PetscCall(PetscMalloc1(pgmres->max_k + 2,&pgmres->orthogwork)); 138 work = pgmres->orthogwork; 139 /* Apply correction computed by the VecMDot in the last iteration to Znext. The original form is 140 * 141 * Znext -= sum_{j=0}^{i-1} Z[j+1] * H[j,i-1] 142 * 143 * where 144 * 145 * Z[j] = sum_{k=0}^j V[k] * H[k,j-1] 146 * 147 * substituting 148 * 149 * Znext -= sum_{j=0}^{i-1} sum_{k=0}^{j+1} V[k] * H[k,j] * H[j,i-1] 150 * 151 * rearranging the iteration space from row-column to column-row 152 * 153 * Znext -= sum_{k=0}^i sum_{j=k-1}^{i-1} V[k] * H[k,j] * H[j,i-1] 154 * 155 * Note that column it-1 of HH is correct. For all previous columns, we must look at HES because HH has already 156 * been transformed to upper triangular form. 157 */ 158 for (k=0; k<it+1; k++) { 159 work[k] = 0; 160 for (j=PetscMax(0,k-1); j<it-1; j++) work[k] -= *HES(k,j) * *HH(j,it-1); 161 } 162 PetscCall(VecMAXPY(Znext,it+1,work,&VEC_VV(0))); 163 PetscCall(VecAXPY(Znext,-*HH(it-1,it-1),Zcur)); 164 165 /* Orthogonalize Zcur against existing basis vectors. */ 166 for (k=0; k<it; k++) work[k] = -*HH(k,it-1); 167 PetscCall(VecMAXPY(Zcur,it,work,&VEC_VV(0))); 168 /* Zcur is now orthogonal, and will be referred to as VEC_VV(it) again, though it is still not normalized. */ 169 /* Begin computing the norm of the new vector, will be normalized after the MatMult in the next iteration. */ 170 PetscCall(VecNormBegin(VEC_VV(it),NORM_2,&newnorm)); 171 } 172 173 /* Compute column of H (to the diagonal, but not the subdiagonal) to be able to orthogonalize the newest vector. */ 174 PetscCall(VecMDotBegin(Znext,it+1,&VEC_VV(0),HH(0,it))); 175 176 /* Start an asynchronous split-mode reduction, the result of the MDot and Norm will be collected on the next iteration. */ 177 PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Znext))); 178 } 179 180 if (itcount) *itcount = it-1; /* Number of iterations actually completed. */ 181 182 /* 183 Down here we have to solve for the "best" coefficients of the Krylov 184 columns, add the solution values together, and possibly unwind the 185 preconditioning from the solution 186 */ 187 /* Form the solution (or the solution so far) */ 188 PetscCall(KSPPGMRESBuildSoln(RS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-2)); 189 PetscFunctionReturn(0); 190 } 191 192 /* 193 KSPSolve_PGMRES - This routine applies the PGMRES method. 194 195 Input Parameter: 196 . ksp - the Krylov space object that was set to use pgmres 197 198 Output Parameter: 199 . outits - number of iterations used 200 201 */ 202 static PetscErrorCode KSPSolve_PGMRES(KSP ksp) 203 { 204 PetscInt its,itcount; 205 KSP_PGMRES *pgmres = (KSP_PGMRES*)ksp->data; 206 PetscBool guess_zero = ksp->guess_zero; 207 208 PetscFunctionBegin; 209 PetscCheck(!ksp->calc_sings || pgmres->Rsvd,PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called"); 210 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 211 ksp->its = 0; 212 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 213 214 itcount = 0; 215 ksp->reason = KSP_CONVERGED_ITERATING; 216 while (!ksp->reason) { 217 PetscCall(KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs)); 218 PetscCall(KSPPGMRESCycle(&its,ksp)); 219 itcount += its; 220 if (itcount >= ksp->max_it) { 221 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS; 222 break; 223 } 224 ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */ 225 } 226 ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */ 227 PetscFunctionReturn(0); 228 } 229 230 static PetscErrorCode KSPDestroy_PGMRES(KSP ksp) 231 { 232 PetscFunctionBegin; 233 PetscCall(KSPDestroy_GMRES(ksp)); 234 PetscFunctionReturn(0); 235 } 236 237 /* 238 KSPPGMRESBuildSoln - create the solution from the starting vector and the 239 current iterates. 240 241 Input parameters: 242 nrs - work area of size it + 1. 243 vguess - index of initial guess 244 vdest - index of result. Note that vguess may == vdest (replace 245 guess with the solution). 246 it - HH upper triangular part is a block of size (it+1) x (it+1) 247 248 This is an internal routine that knows about the PGMRES internals. 249 */ 250 static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it) 251 { 252 PetscScalar tt; 253 PetscInt k,j; 254 KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data); 255 256 PetscFunctionBegin; 257 /* Solve for solution vector that minimizes the residual */ 258 259 if (it < 0) { /* no pgmres steps have been performed */ 260 PetscCall(VecCopy(vguess,vdest)); /* VecCopy() is smart, exits immediately if vguess == vdest */ 261 PetscFunctionReturn(0); 262 } 263 264 /* solve the upper triangular system - RS is the right side and HH is 265 the upper triangular matrix - put soln in nrs */ 266 if (*HH(it,it) != 0.0) nrs[it] = *RS(it) / *HH(it,it); 267 else nrs[it] = 0.0; 268 269 for (k=it-1; k>=0; k--) { 270 tt = *RS(k); 271 for (j=k+1; j<=it; j++) tt -= *HH(k,j) * nrs[j]; 272 nrs[k] = tt / *HH(k,k); 273 } 274 275 /* Accumulate the correction to the solution of the preconditioned problem in TEMP */ 276 PetscCall(VecZeroEntries(VEC_TEMP)); 277 PetscCall(VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0))); 278 PetscCall(KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP)); 279 /* add solution to previous solution */ 280 if (vdest == vguess) { 281 PetscCall(VecAXPY(vdest,1.0,VEC_TEMP)); 282 } else { 283 PetscCall(VecWAXPY(vdest,1.0,VEC_TEMP,vguess)); 284 } 285 PetscFunctionReturn(0); 286 } 287 288 /* 289 290 KSPPGMRESUpdateHessenberg - Do the scalar work for the orthogonalization. 291 Return new residual. 292 293 input parameters: 294 295 . ksp - Krylov space object 296 . it - plane rotations are applied to the (it+1)th column of the 297 modified hessenberg (i.e. HH(:,it)) 298 . hapend - PETSC_FALSE not happy breakdown ending. 299 300 output parameters: 301 . res - the new residual 302 303 */ 304 /* 305 . it - column of the Hessenberg that is complete, PGMRES is actually computing two columns ahead of this 306 */ 307 static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool *hapend,PetscReal *res) 308 { 309 PetscScalar *hh,*cc,*ss,*rs; 310 PetscInt j; 311 PetscReal hapbnd; 312 KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data); 313 314 PetscFunctionBegin; 315 hh = HH(0,it); /* pointer to beginning of column to update */ 316 cc = CC(0); /* beginning of cosine rotations */ 317 ss = SS(0); /* beginning of sine rotations */ 318 rs = RS(0); /* right hand side of least squares system */ 319 320 /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */ 321 for (j=0; j<=it+1; j++) *HES(j,it) = hh[j]; 322 323 /* check for the happy breakdown */ 324 hapbnd = PetscMin(PetscAbsScalar(hh[it+1] / rs[it]),pgmres->haptol); 325 if (PetscAbsScalar(hh[it+1]) < hapbnd) { 326 PetscCall(PetscInfo(ksp,"Detected happy breakdown, current hapbnd = %14.12e H(%D,%D) = %14.12e\n",(double)hapbnd,it+1,it,(double)PetscAbsScalar(*HH(it+1,it)))); 327 *hapend = PETSC_TRUE; 328 } 329 330 /* Apply all the previously computed plane rotations to the new column 331 of the Hessenberg matrix */ 332 /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta), 333 and some refs have [c s ; -conj(s) c] (don't be confused!) */ 334 335 for (j=0; j<it; j++) { 336 PetscScalar hhj = hh[j]; 337 hh[j] = PetscConj(cc[j])*hhj + ss[j]*hh[j+1]; 338 hh[j+1] = -ss[j] *hhj + cc[j]*hh[j+1]; 339 } 340 341 /* 342 compute the new plane rotation, and apply it to: 343 1) the right-hand-side of the Hessenberg system (RS) 344 note: it affects RS(it) and RS(it+1) 345 2) the new column of the Hessenberg matrix 346 note: it affects HH(it,it) which is currently pointed to 347 by hh and HH(it+1, it) (*(hh+1)) 348 thus obtaining the updated value of the residual... 349 */ 350 351 /* compute new plane rotation */ 352 353 if (!*hapend) { 354 PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it+1]))); 355 if (delta == 0.0) { 356 ksp->reason = KSP_DIVERGED_NULL; 357 PetscFunctionReturn(0); 358 } 359 360 cc[it] = hh[it] / delta; /* new cosine value */ 361 ss[it] = hh[it+1] / delta; /* new sine value */ 362 363 hh[it] = PetscConj(cc[it])*hh[it] + ss[it]*hh[it+1]; 364 rs[it+1] = -ss[it]*rs[it]; 365 rs[it] = PetscConj(cc[it])*rs[it]; 366 *res = PetscAbsScalar(rs[it+1]); 367 } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply 368 another rotation matrix (so RH doesn't change). The new residual is 369 always the new sine term times the residual from last time (RS(it)), 370 but now the new sine rotation would be zero...so the residual should 371 be zero...so we will multiply "zero" by the last residual. This might 372 not be exactly what we want to do here -could just return "zero". */ 373 374 *res = 0.0; 375 } 376 PetscFunctionReturn(0); 377 } 378 379 /* 380 KSPBuildSolution_PGMRES 381 382 Input Parameter: 383 . ksp - the Krylov space object 384 . ptr- 385 386 Output Parameter: 387 . result - the solution 388 389 Note: this calls KSPPGMRESBuildSoln - the same function that KSPPGMRESCycle 390 calls directly. 391 392 */ 393 PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp,Vec ptr,Vec *result) 394 { 395 KSP_PGMRES *pgmres = (KSP_PGMRES*)ksp->data; 396 397 PetscFunctionBegin; 398 if (!ptr) { 399 if (!pgmres->sol_temp) { 400 PetscCall(VecDuplicate(ksp->vec_sol,&pgmres->sol_temp)); 401 PetscCall(PetscLogObjectParent((PetscObject)ksp,(PetscObject)pgmres->sol_temp)); 402 } 403 ptr = pgmres->sol_temp; 404 } 405 if (!pgmres->nrs) { 406 /* allocate the work area */ 407 PetscCall(PetscMalloc1(pgmres->max_k,&pgmres->nrs)); 408 PetscCall(PetscLogObjectMemory((PetscObject)ksp,pgmres->max_k*sizeof(PetscScalar))); 409 } 410 411 PetscCall(KSPPGMRESBuildSoln(pgmres->nrs,ksp->vec_sol,ptr,ksp,pgmres->it)); 412 if (result) *result = ptr; 413 PetscFunctionReturn(0); 414 } 415 416 PetscErrorCode KSPSetFromOptions_PGMRES(PetscOptionItems *PetscOptionsObject,KSP ksp) 417 { 418 PetscFunctionBegin; 419 PetscCall(KSPSetFromOptions_GMRES(PetscOptionsObject,ksp)); 420 PetscCall(PetscOptionsHead(PetscOptionsObject,"KSP pipelined GMRES Options")); 421 PetscCall(PetscOptionsTail()); 422 PetscFunctionReturn(0); 423 } 424 425 PetscErrorCode KSPReset_PGMRES(KSP ksp) 426 { 427 PetscFunctionBegin; 428 PetscCall(KSPReset_GMRES(ksp)); 429 PetscFunctionReturn(0); 430 } 431 432 /*MC 433 KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method. 434 435 Options Database Keys: 436 + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against 437 . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence) 438 . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of 439 vectors are allocated as needed) 440 . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default) 441 . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower) 442 . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the 443 stability of the classical Gram-Schmidt orthogonalization. 444 - -ksp_gmres_krylov_monitor - plot the Krylov space generated 445 446 Level: beginner 447 448 Notes: 449 MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods. 450 See the FAQ on the PETSc website for details. 451 452 Reference: 453 Ghysels, Ashby, Meerbergen, Vanroose, Hiding global communication latencies in the GMRES algorithm on massively parallel machines, 2012. 454 455 Developer Notes: 456 This object is subclassed off of KSPGMRES 457 458 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPLGMRES, KSPPIPECG, KSPPIPECR, 459 KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(), 460 KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(), 461 KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov() 462 M*/ 463 464 PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp) 465 { 466 KSP_PGMRES *pgmres; 467 468 PetscFunctionBegin; 469 PetscCall(PetscNewLog(ksp,&pgmres)); 470 471 ksp->data = (void*)pgmres; 472 ksp->ops->buildsolution = KSPBuildSolution_PGMRES; 473 ksp->ops->setup = KSPSetUp_PGMRES; 474 ksp->ops->solve = KSPSolve_PGMRES; 475 ksp->ops->reset = KSPReset_PGMRES; 476 ksp->ops->destroy = KSPDestroy_PGMRES; 477 ksp->ops->view = KSPView_GMRES; 478 ksp->ops->setfromoptions = KSPSetFromOptions_PGMRES; 479 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES; 480 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES; 481 482 PetscCall(KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3)); 483 PetscCall(KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2)); 484 PetscCall(KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1)); 485 486 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES)); 487 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES)); 488 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES)); 489 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES)); 490 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES)); 491 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES)); 492 PetscCall(PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES)); 493 494 pgmres->nextra_vecs = 1; 495 pgmres->haptol = 1.0e-30; 496 pgmres->q_preallocate = 0; 497 pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS; 498 pgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization; 499 pgmres->nrs = NULL; 500 pgmres->sol_temp = NULL; 501 pgmres->max_k = PGMRES_DEFAULT_MAXK; 502 pgmres->Rsvd = NULL; 503 pgmres->orthogwork = NULL; 504 pgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER; 505 PetscFunctionReturn(0); 506 } 507