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