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 (ksp->reason) break; 78 if (it < pgmres->max_k + 1) { /* Monitor if we are not done or still iterating, but not before a restart. */ 79 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm)); 80 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm)); 81 } 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 156 if (ksp->reason || ksp->its == ksp->max_it) { 157 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm)); 158 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm)); 159 } 160 PetscFunctionReturn(PETSC_SUCCESS); 161 } 162 163 static PetscErrorCode KSPSolve_PGMRES(KSP ksp) 164 { 165 PetscInt its, itcount; 166 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data; 167 PetscBool guess_zero = ksp->guess_zero; 168 169 PetscFunctionBegin; 170 PetscCheck(!ksp->calc_sings || pgmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called"); 171 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 172 ksp->its = 0; 173 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 174 175 itcount = 0; 176 ksp->reason = KSP_CONVERGED_ITERATING; 177 while (!ksp->reason) { 178 PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs)); 179 PetscCall(KSPPGMRESCycle(&its, ksp)); 180 itcount += its; 181 if (itcount >= ksp->max_it) { 182 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS; 183 break; 184 } 185 ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */ 186 } 187 ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */ 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 static PetscErrorCode KSPDestroy_PGMRES(KSP ksp) 192 { 193 PetscFunctionBegin; 194 PetscCall(KSPDestroy_GMRES(ksp)); 195 PetscFunctionReturn(PETSC_SUCCESS); 196 } 197 198 static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it) 199 { 200 PetscScalar tt; 201 PetscInt k, j; 202 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data; 203 204 PetscFunctionBegin; 205 /* Solve for solution vector that minimizes the residual */ 206 207 if (it < 0) { /* no pgmres steps have been performed */ 208 PetscCall(VecCopy(vguess, vdest)); /* VecCopy() is smart, exits immediately if vguess == vdest */ 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 212 /* solve the upper triangular system - RS is the right side and HH is 213 the upper triangular matrix - put soln in nrs */ 214 if (*HH(it, it) != 0.0) nrs[it] = *RS(it) / *HH(it, it); 215 else nrs[it] = 0.0; 216 217 for (k = it - 1; k >= 0; k--) { 218 tt = *RS(k); 219 for (j = k + 1; j <= it; j++) tt -= *HH(k, j) * nrs[j]; 220 nrs[k] = tt / *HH(k, k); 221 } 222 223 /* Accumulate the correction to the solution of the preconditioned problem in TEMP */ 224 PetscCall(VecMAXPBY(VEC_TEMP, it + 1, nrs, 0, &VEC_VV(0))); 225 PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP)); 226 /* add solution to previous solution */ 227 if (vdest == vguess) { 228 PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP)); 229 } else { 230 PetscCall(VecWAXPY(vdest, 1.0, VEC_TEMP, vguess)); 231 } 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 235 static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool *hapend, PetscReal *res) 236 { 237 PetscScalar *hh, *cc, *ss, *rs; 238 PetscInt j; 239 PetscReal hapbnd; 240 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data; 241 242 PetscFunctionBegin; 243 hh = HH(0, it); /* pointer to beginning of column to update */ 244 cc = CC(0); /* beginning of cosine rotations */ 245 ss = SS(0); /* beginning of sine rotations */ 246 rs = RS(0); /* right-hand side of least squares system */ 247 248 /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */ 249 for (j = 0; j <= it + 1; j++) *HES(j, it) = hh[j]; 250 251 /* check for the happy breakdown */ 252 hapbnd = PetscMin(PetscAbsScalar(hh[it + 1] / rs[it]), pgmres->haptol); 253 if (PetscAbsScalar(hh[it + 1]) < hapbnd) { 254 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)))); 255 *hapend = PETSC_TRUE; 256 } 257 258 /* Apply all the previously computed plane rotations to the new column of the Hessenberg matrix */ 259 /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta), 260 and some refs have [c s ; -conj(s) c] (don't be confused!) */ 261 262 for (j = 0; j < it; j++) { 263 PetscScalar hhj = hh[j]; 264 hh[j] = PetscConj(cc[j]) * hhj + ss[j] * hh[j + 1]; 265 hh[j + 1] = -ss[j] * hhj + cc[j] * hh[j + 1]; 266 } 267 268 /* 269 compute the new plane rotation, and apply it to: 270 1) the right-hand side of the Hessenberg system (RS) 271 note: it affects RS(it) and RS(it+1) 272 2) the new column of the Hessenberg matrix 273 note: it affects HH(it,it) which is currently pointed to 274 by hh and HH(it+1, it) (*(hh+1)) 275 thus obtaining the updated value of the residual... 276 */ 277 278 /* compute new plane rotation */ 279 280 if (!*hapend) { 281 PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it + 1]))); 282 if (delta == 0.0) { 283 ksp->reason = KSP_DIVERGED_NULL; 284 PetscFunctionReturn(PETSC_SUCCESS); 285 } 286 287 cc[it] = hh[it] / delta; /* new cosine value */ 288 ss[it] = hh[it + 1] / delta; /* new sine value */ 289 290 hh[it] = PetscConj(cc[it]) * hh[it] + ss[it] * hh[it + 1]; 291 rs[it + 1] = -ss[it] * rs[it]; 292 rs[it] = PetscConj(cc[it]) * rs[it]; 293 *res = PetscAbsScalar(rs[it + 1]); 294 } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply 295 another rotation matrix (so RH doesn't change). The new residual is 296 always the new sine term times the residual from last time (RS(it)), 297 but now the new sine rotation would be zero...so the residual should 298 be zero...so we will multiply "zero" by the last residual. This might 299 not be exactly what we want to do here -could just return "zero". */ 300 *res = 0.0; 301 } 302 PetscFunctionReturn(PETSC_SUCCESS); 303 } 304 305 static PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp, Vec ptr, Vec *result) 306 { 307 KSP_PGMRES *pgmres = (KSP_PGMRES *)ksp->data; 308 309 PetscFunctionBegin; 310 if (!ptr) { 311 if (!pgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &pgmres->sol_temp)); 312 ptr = pgmres->sol_temp; 313 } 314 if (!pgmres->nrs) { 315 /* allocate the work area */ 316 PetscCall(PetscMalloc1(pgmres->max_k, &pgmres->nrs)); 317 } 318 319 PetscCall(KSPPGMRESBuildSoln(pgmres->nrs, ksp->vec_sol, ptr, ksp, pgmres->it)); 320 if (result) *result = ptr; 321 PetscFunctionReturn(PETSC_SUCCESS); 322 } 323 324 static PetscErrorCode KSPSetFromOptions_PGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject) 325 { 326 PetscFunctionBegin; 327 PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject)); 328 PetscOptionsHeadBegin(PetscOptionsObject, "KSP pipelined GMRES Options"); 329 PetscOptionsHeadEnd(); 330 PetscFunctionReturn(PETSC_SUCCESS); 331 } 332 333 static PetscErrorCode KSPReset_PGMRES(KSP ksp) 334 { 335 PetscFunctionBegin; 336 PetscCall(KSPReset_GMRES(ksp)); 337 PetscFunctionReturn(PETSC_SUCCESS); 338 } 339 340 /*MC 341 KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method {cite}`ghyselsashbymeerbergenvanroose2013`. [](sec_pipelineksp) 342 343 Options Database Keys: 344 + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against 345 . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence) 346 . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially 347 (otherwise groups of vectors are allocated as needed) 348 . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize 349 against the Krylov space (fast) (the default) 350 . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower) 351 . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the 352 stability of the classical Gram-Schmidt orthogonalization. 353 - -ksp_gmres_krylov_monitor - plot the Krylov space generated 354 355 Level: beginner 356 357 Note: 358 MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods. 359 See [](doc_faq_pipelined) 360 361 Developer Note: 362 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 363 364 .seealso: [](ch_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPLGMRES`, `KSPPIPECG`, `KSPPIPECR`, 365 `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`, 366 `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`, 367 `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()` 368 M*/ 369 370 PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp) 371 { 372 KSP_PGMRES *pgmres; 373 374 PetscFunctionBegin; 375 PetscCall(PetscNew(&pgmres)); 376 377 ksp->data = (void *)pgmres; 378 ksp->ops->buildsolution = KSPBuildSolution_PGMRES; 379 ksp->ops->setup = KSPSetUp_PGMRES; 380 ksp->ops->solve = KSPSolve_PGMRES; 381 ksp->ops->reset = KSPReset_PGMRES; 382 ksp->ops->destroy = KSPDestroy_PGMRES; 383 ksp->ops->view = KSPView_GMRES; 384 ksp->ops->setfromoptions = KSPSetFromOptions_PGMRES; 385 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES; 386 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES; 387 388 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3)); 389 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2)); 390 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1)); 391 392 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES)); 393 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES)); 394 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetOrthogonalization_C", KSPGMRESGetOrthogonalization_GMRES)); 395 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES)); 396 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES)); 397 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES)); 398 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetCGSRefinementType_C", KSPGMRESGetCGSRefinementType_GMRES)); 399 400 pgmres->nextra_vecs = 1; 401 pgmres->haptol = 1.0e-30; 402 pgmres->q_preallocate = 0; 403 pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS; 404 pgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization; 405 pgmres->nrs = NULL; 406 pgmres->sol_temp = NULL; 407 pgmres->max_k = PGMRES_DEFAULT_MAXK; 408 pgmres->Rsvd = NULL; 409 pgmres->orthogwork = NULL; 410 pgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER; 411 PetscFunctionReturn(PETSC_SUCCESS); 412 } 413