1 #include <petsc/private/kspimpl.h> 2 3 typedef struct { 4 PetscReal haptol; 5 } KSP_SYMMLQ; 6 7 static PetscErrorCode KSPSetUp_SYMMLQ(KSP ksp) 8 { 9 PetscFunctionBegin; 10 PetscCall(KSPSetWorkVecs(ksp, 9)); 11 PetscFunctionReturn(PETSC_SUCCESS); 12 } 13 14 static PetscErrorCode KSPSolve_SYMMLQ(KSP ksp) 15 { 16 PetscInt i; 17 PetscScalar alpha, beta, ibeta, betaold, beta1, ceta = 0, ceta_oold = 0.0, ceta_old = 0.0, ceta_bar; 18 PetscScalar c = 1.0, cold = 1.0, s = 0.0, sold = 0.0, coold, soold, rho0, rho1, rho2, rho3; 19 PetscScalar dp = 0.0; 20 PetscReal np = 0.0, s_prod; 21 Vec X, B, R, Z, U, V, W, UOLD, VOLD, Wbar; 22 Mat Amat, Pmat; 23 KSP_SYMMLQ *symmlq = (KSP_SYMMLQ *)ksp->data; 24 PetscBool diagonalscale; 25 26 PetscFunctionBegin; 27 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale)); 28 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name); 29 30 X = ksp->vec_sol; 31 B = ksp->vec_rhs; 32 R = ksp->work[0]; 33 Z = ksp->work[1]; 34 U = ksp->work[2]; 35 V = ksp->work[3]; 36 W = ksp->work[4]; 37 UOLD = ksp->work[5]; 38 VOLD = ksp->work[6]; 39 Wbar = ksp->work[7]; 40 41 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat)); 42 43 ksp->its = 0; 44 45 PetscCall(VecSet(UOLD, 0.0)); /* u_old <- zeros; */ 46 PetscCall(VecCopy(UOLD, VOLD)); /* v_old <- u_old; */ 47 PetscCall(VecCopy(UOLD, W)); /* w <- u_old; */ 48 PetscCall(VecCopy(UOLD, Wbar)); /* w_bar <- u_old; */ 49 if (!ksp->guess_zero) { 50 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */ 51 PetscCall(VecAYPX(R, -1.0, B)); 52 } else { 53 PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */ 54 } 55 56 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */ 57 PetscCall(VecDot(R, Z, &dp)); /* dp = r'*z; */ 58 KSPCheckDot(ksp, dp); 59 if (PetscAbsScalar(dp) < symmlq->haptol) { 60 PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol)); 61 ksp->rnorm = 0.0; /* what should we really put here? */ 62 ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN; /* bugfix proposed by Lourens (lourens.vanzanen@shell.com) */ 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 #if !defined(PETSC_USE_COMPLEX) 67 if (dp < 0.0) { 68 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 69 PetscFunctionReturn(PETSC_SUCCESS); 70 } 71 #endif 72 dp = PetscSqrtScalar(dp); 73 beta = dp; /* beta <- sqrt(r'*z) */ 74 beta1 = beta; 75 s_prod = PetscAbsScalar(beta1); 76 77 PetscCall(VecCopy(R, V)); /* v <- r; */ 78 PetscCall(VecCopy(Z, U)); /* u <- z; */ 79 ibeta = 1.0 / beta; 80 PetscCall(VecScale(V, ibeta)); /* v <- ibeta*v; */ 81 PetscCall(VecScale(U, ibeta)); /* u <- ibeta*u; */ 82 PetscCall(VecCopy(U, Wbar)); /* w_bar <- u; */ 83 if (ksp->normtype != KSP_NORM_NONE) { 84 PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */ 85 KSPCheckNorm(ksp, np); 86 } 87 PetscCall(KSPLogResidualHistory(ksp, np)); 88 PetscCall(KSPMonitor(ksp, 0, np)); 89 ksp->rnorm = np; 90 PetscCall((*ksp->converged)(ksp, 0, np, &ksp->reason, ksp->cnvP)); /* test for convergence */ 91 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS); 92 93 i = 0; 94 ceta = 0.; 95 do { 96 ksp->its = i + 1; 97 98 /* Update */ 99 if (ksp->its > 1) { 100 PetscCall(VecCopy(V, VOLD)); /* v_old <- v; */ 101 PetscCall(VecCopy(U, UOLD)); /* u_old <- u; */ 102 103 PetscCall(VecCopy(R, V)); 104 PetscCall(VecScale(V, 1.0 / beta)); /* v <- ibeta*r; */ 105 PetscCall(VecCopy(Z, U)); 106 PetscCall(VecScale(U, 1.0 / beta)); /* u <- ibeta*z; */ 107 108 PetscCall(VecCopy(Wbar, W)); 109 PetscCall(VecScale(W, c)); 110 PetscCall(VecAXPY(W, s, U)); /* w <- c*w_bar + s*u; (w_k) */ 111 PetscCall(VecScale(Wbar, -s)); 112 PetscCall(VecAXPY(Wbar, c, U)); /* w_bar <- -s*w_bar + c*u; (w_bar_(k+1)) */ 113 PetscCall(VecAXPY(X, ceta, W)); /* x <- x + ceta * w; (xL_k) */ 114 115 ceta_oold = ceta_old; 116 ceta_old = ceta; 117 } 118 119 /* Lanczos */ 120 PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- Amat*u; */ 121 PetscCall(VecDot(U, R, &alpha)); /* alpha <- u'*r; */ 122 PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r; */ 123 124 PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha* v; */ 125 PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha* u; */ 126 PetscCall(VecAXPY(R, -beta, VOLD)); /* r <- r - beta * v_old; */ 127 PetscCall(VecAXPY(Z, -beta, UOLD)); /* z <- z - beta * u_old; */ 128 betaold = beta; /* beta_k */ 129 PetscCall(VecDot(R, Z, &dp)); /* dp <- r'*z; */ 130 KSPCheckDot(ksp, dp); 131 if (PetscAbsScalar(dp) < symmlq->haptol) { 132 PetscCall(PetscInfo(ksp, "Detected happy breakdown %g tolerance %g\n", (double)PetscAbsScalar(dp), (double)symmlq->haptol)); 133 dp = 0.0; 134 } 135 136 #if !defined(PETSC_USE_COMPLEX) 137 if (dp < 0.0) { 138 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 139 break; 140 } 141 #endif 142 beta = PetscSqrtScalar(dp); /* beta = sqrt(dp); */ 143 144 /* QR factorization */ 145 coold = cold; 146 cold = c; 147 soold = sold; 148 sold = s; 149 rho0 = cold * alpha - coold * sold * betaold; /* gamma_bar */ 150 rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta); /* gamma */ 151 rho2 = sold * alpha + coold * cold * betaold; /* delta */ 152 rho3 = soold * betaold; /* epsilon */ 153 154 /* Givens rotation: [c -s; s c] (different from the Reference!) */ 155 c = rho0 / rho1; 156 s = beta / rho1; 157 158 if (ksp->its == 1) ceta = beta1 / rho1; 159 else ceta = -(rho2 * ceta_old + rho3 * ceta_oold) / rho1; 160 161 s_prod = s_prod * PetscAbsScalar(s); 162 if (c == 0.0) np = s_prod * 1.e16; 163 else np = s_prod / PetscAbsScalar(c); /* residual norm for xc_k (CGNORM) */ 164 165 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np; 166 else ksp->rnorm = 0.0; 167 PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm)); 168 PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm)); 169 PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */ 170 if (ksp->reason) break; 171 i++; 172 } while (i < ksp->max_it); 173 174 /* move to the CG point: xc_(k+1) */ 175 if (c == 0.0) ceta_bar = ceta * 1.e15; 176 else ceta_bar = ceta / c; 177 178 PetscCall(VecAXPY(X, ceta_bar, Wbar)); /* x <- x + ceta_bar*w_bar */ 179 180 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 181 PetscFunctionReturn(PETSC_SUCCESS); 182 } 183 184 /*MC 185 KSPSYMMLQ - Implements the SYMMLQ method {cite}`paige.saunders:solution`. 186 187 Level: beginner 188 189 Notes: 190 The operator and the preconditioner must be symmetric for this method. 191 192 The preconditioner must be POSITIVE-DEFINITE. 193 194 Supports only left preconditioning. 195 196 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP` 197 M*/ 198 PETSC_EXTERN PetscErrorCode KSPCreate_SYMMLQ(KSP ksp) 199 { 200 KSP_SYMMLQ *symmlq; 201 202 PetscFunctionBegin; 203 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3)); 204 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1)); 205 206 PetscCall(PetscNew(&symmlq)); 207 symmlq->haptol = 1.e-18; 208 ksp->data = (void *)symmlq; 209 210 /* 211 Sets the functions that are associated with this data structure 212 (in C++ this is the same as defining virtual functions) 213 */ 214 ksp->ops->setup = KSPSetUp_SYMMLQ; 215 ksp->ops->solve = KSPSolve_SYMMLQ; 216 ksp->ops->destroy = KSPDestroyDefault; 217 ksp->ops->setfromoptions = NULL; 218 ksp->ops->buildsolution = KSPBuildSolutionDefault; 219 ksp->ops->buildresidual = KSPBuildResidualDefault; 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222