1 2 /* 3 Note that for the complex numbers version, the VecDot() arguments 4 within the code MUST remain in the order given for correct computation 5 of inner products. 6 */ 7 #include <petsc/private/kspimpl.h> 8 9 static PetscErrorCode KSPSetUp_CGS(KSP ksp) 10 { 11 PetscFunctionBegin; 12 PetscCall(KSPSetWorkVecs(ksp, 7)); 13 PetscFunctionReturn(PETSC_SUCCESS); 14 } 15 16 static PetscErrorCode KSPSolve_CGS(KSP ksp) 17 { 18 PetscInt i; 19 PetscScalar rho, rhoold, a, s, b; 20 Vec X, B, V, P, R, RP, T, Q, U, AUQ; 21 PetscReal dp = 0.0; 22 PetscBool diagonalscale; 23 24 PetscFunctionBegin; 25 /* not sure what residual norm it does use, should use for right preconditioning */ 26 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 RP = ksp->work[1]; 34 V = ksp->work[2]; 35 T = ksp->work[3]; 36 Q = ksp->work[4]; 37 P = ksp->work[5]; 38 U = ksp->work[6]; 39 AUQ = V; 40 41 /* Compute initial preconditioned residual */ 42 PetscCall(KSPInitialResidual(ksp, X, V, T, R, B)); 43 44 /* Test for nothing to do */ 45 if (ksp->normtype != KSP_NORM_NONE) { 46 PetscCall(VecNorm(R, NORM_2, &dp)); 47 KSPCheckNorm(ksp, dp); 48 if (ksp->normtype == KSP_NORM_NATURAL) dp *= dp; 49 } else dp = 0.0; 50 51 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 52 ksp->its = 0; 53 ksp->rnorm = dp; 54 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 55 PetscCall(KSPLogResidualHistory(ksp, dp)); 56 PetscCall(KSPMonitor(ksp, 0, dp)); 57 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); 58 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS); 59 60 /* Make the initial Rp == R */ 61 PetscCall(VecCopy(R, RP)); 62 /* added for Fidap */ 63 /* Penalize Startup - Isaac Hasbani Trick for CGS 64 Since most initial conditions result in a mostly 0 residual, 65 we change all the 0 values in the vector RP to the maximum. 66 */ 67 if (ksp->normtype == KSP_NORM_NATURAL) { 68 PetscReal vr0max; 69 PetscScalar *tmp_RP = NULL; 70 PetscInt numnp = 0, *max_pos = NULL; 71 PetscCall(VecMax(RP, max_pos, &vr0max)); 72 PetscCall(VecGetArray(RP, &tmp_RP)); 73 PetscCall(VecGetLocalSize(RP, &numnp)); 74 for (i = 0; i < numnp; i++) { 75 if (tmp_RP[i] == 0.0) tmp_RP[i] = vr0max; 76 } 77 PetscCall(VecRestoreArray(RP, &tmp_RP)); 78 } 79 /* end of addition for Fidap */ 80 81 /* Set the initial conditions */ 82 PetscCall(VecDot(R, RP, &rhoold)); /* rhoold = (r,rp) */ 83 PetscCall(VecCopy(R, U)); 84 PetscCall(VecCopy(R, P)); 85 PetscCall(KSP_PCApplyBAorAB(ksp, P, V, T)); 86 87 i = 0; 88 do { 89 PetscCall(VecDot(V, RP, &s)); /* s <- (v,rp) */ 90 KSPCheckDot(ksp, s); 91 a = rhoold / s; /* a <- rho / s */ 92 PetscCall(VecWAXPY(Q, -a, V, U)); /* q <- u - a v */ 93 PetscCall(VecWAXPY(T, 1.0, U, Q)); /* t <- u + q */ 94 PetscCall(VecAXPY(X, a, T)); /* x <- x + a (u + q) */ 95 PetscCall(KSP_PCApplyBAorAB(ksp, T, AUQ, U)); 96 PetscCall(VecAXPY(R, -a, AUQ)); /* r <- r - a K (u + q) */ 97 PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp) */ 98 KSPCheckDot(ksp, rho); 99 if (ksp->normtype == KSP_NORM_NATURAL) { 100 dp = PetscAbsScalar(rho); 101 } else if (ksp->normtype != KSP_NORM_NONE) { 102 PetscCall(VecNorm(R, NORM_2, &dp)); 103 KSPCheckNorm(ksp, dp); 104 } else dp = 0.0; 105 106 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 107 ksp->its++; 108 ksp->rnorm = dp; 109 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 110 PetscCall(KSPLogResidualHistory(ksp, dp)); 111 PetscCall(KSPMonitor(ksp, i + 1, dp)); 112 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP)); 113 if (ksp->reason) break; 114 115 b = rho / rhoold; /* b <- rho / rhoold */ 116 PetscCall(VecWAXPY(U, b, Q, R)); /* u <- r + b q */ 117 PetscCall(VecAXPY(Q, b, P)); 118 PetscCall(VecWAXPY(P, b, Q, U)); /* p <- u + b(q + b p) */ 119 PetscCall(KSP_PCApplyBAorAB(ksp, P, V, Q)); /* v <- K p */ 120 rhoold = rho; 121 i++; 122 } while (i < ksp->max_it); 123 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 124 125 PetscCall(KSPUnwindPreconditioner(ksp, X, T)); 126 PetscFunctionReturn(PETSC_SUCCESS); 127 } 128 129 /*MC 130 KSPCGS - This code implements the CGS (Conjugate Gradient Squared) method. 131 132 Level: beginner 133 134 Notes: 135 Does not require a symmetric matrix. Does not apply transpose of the matrix. 136 137 Supports left and right preconditioning, but not symmetric. 138 139 Developer Note: 140 Has this weird support for doing the convergence test with the natural norm, I assume this works only with 141 no preconditioning and symmetric positive definite operator. 142 143 References: 144 . * - Sonneveld, 1989. 145 146 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBCGS`, `KSPSetPCSide()` 147 M*/ 148 PETSC_EXTERN PetscErrorCode KSPCreate_CGS(KSP ksp) 149 { 150 PetscFunctionBegin; 151 ksp->data = (void *)0; 152 153 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3)); 154 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2)); 155 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2)); 156 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_RIGHT, 2)); 157 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1)); 158 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1)); 159 160 ksp->ops->setup = KSPSetUp_CGS; 161 ksp->ops->solve = KSPSolve_CGS; 162 ksp->ops->destroy = KSPDestroyDefault; 163 ksp->ops->buildsolution = KSPBuildSolutionDefault; 164 ksp->ops->buildresidual = KSPBuildResidualDefault; 165 ksp->ops->setfromoptions = NULL; 166 ksp->ops->view = NULL; 167 PetscFunctionReturn(PETSC_SUCCESS); 168 } 169