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