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