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