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