1 #include <petsc/private/kspimpl.h> 2 3 static PetscErrorCode KSPSetUp_CR(KSP ksp) 4 { 5 PetscFunctionBegin; 6 PetscCall(KSPSetWorkVecs(ksp, 6)); 7 PetscFunctionReturn(PETSC_SUCCESS); 8 } 9 10 static PetscErrorCode KSPSolve_CR(KSP ksp) 11 { 12 PetscInt i = 0; 13 PetscReal dp; 14 PetscScalar ai, bi; 15 PetscScalar apq, btop, bbot; 16 Vec X, B, R, RT, P, AP, ART, Q; 17 Mat Amat, Pmat; 18 19 PetscFunctionBegin; 20 X = ksp->vec_sol; 21 B = ksp->vec_rhs; 22 R = ksp->work[0]; 23 RT = ksp->work[1]; 24 P = ksp->work[2]; 25 AP = ksp->work[3]; 26 ART = ksp->work[4]; 27 Q = ksp->work[5]; 28 29 /* R is the true residual norm, RT is the preconditioned residual norm */ 30 PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat)); 31 if (!ksp->guess_zero) { 32 PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* R <- A*X */ 33 PetscCall(VecAYPX(R, -1.0, B)); /* R <- B-R == B-A*X */ 34 } else { 35 PetscCall(VecCopy(B, R)); /* R <- B (X is 0) */ 36 } 37 /* This may be true only on a subset of MPI ranks; setting it here so it will be detected by the first norm computation below */ 38 if (ksp->reason == KSP_DIVERGED_PC_FAILED) PetscCall(VecSetInf(R)); 39 PetscCall(KSP_PCApply(ksp, R, P)); /* P <- B*R */ 40 PetscCall(KSP_MatMult(ksp, Amat, P, AP)); /* AP <- A*P */ 41 PetscCall(VecCopy(P, RT)); /* RT <- P */ 42 PetscCall(VecCopy(AP, ART)); /* ART <- AP */ 43 PetscCall(VecDotBegin(RT, ART, &btop)); /* (RT,ART) */ 44 45 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 46 PetscCall(VecNormBegin(RT, NORM_2, &dp)); /* dp <- RT'*RT */ 47 PetscCall(VecDotEnd(RT, ART, &btop)); /* (RT,ART) */ 48 PetscCall(VecNormEnd(RT, NORM_2, &dp)); /* dp <- RT'*RT */ 49 KSPCheckNorm(ksp, dp); 50 } else if (ksp->normtype == KSP_NORM_NONE) { 51 dp = 0.0; /* meaningless value that is passed to monitor and convergence test */ 52 PetscCall(VecDotEnd(RT, ART, &btop)); /* (RT,ART) */ 53 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 54 PetscCall(VecNormBegin(R, NORM_2, &dp)); /* dp <- R'*R */ 55 PetscCall(VecDotEnd(RT, ART, &btop)); /* (RT,ART) */ 56 PetscCall(VecNormEnd(R, NORM_2, &dp)); /* dp <- RT'*RT */ 57 KSPCheckNorm(ksp, dp); 58 } else if (ksp->normtype == KSP_NORM_NATURAL) { 59 PetscCall(VecDotEnd(RT, ART, &btop)); /* (RT,ART) */ 60 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 61 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSPNormType of %d not supported", (int)ksp->normtype); 62 if (PetscAbsScalar(btop) < 0.0) { 63 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 64 PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite matrix\n")); 65 PetscFunctionReturn(PETSC_SUCCESS); 66 } 67 68 ksp->its = 0; 69 PetscCall(KSPMonitor(ksp, 0, dp)); 70 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 71 ksp->rnorm = dp; 72 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 73 PetscCall(KSPLogResidualHistory(ksp, dp)); 74 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP)); 75 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS); 76 77 i = 0; 78 do { 79 PetscCall(KSP_PCApply(ksp, AP, Q)); /* Q <- B* AP */ 80 81 PetscCall(VecDot(AP, Q, &apq)); 82 KSPCheckDot(ksp, apq); 83 if (PetscRealPart(apq) <= 0.0) { 84 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 85 PetscCall(PetscInfo(ksp, "KSPSolve_CR:diverging due to indefinite or negative definite PC\n")); 86 break; 87 } 88 ai = btop / apq; /* ai = (RT,ART)/(AP,Q) */ 89 90 PetscCall(VecAXPY(X, ai, P)); /* X <- X + ai*P */ 91 PetscCall(VecAXPY(RT, -ai, Q)); /* RT <- RT - ai*Q */ 92 PetscCall(KSP_MatMult(ksp, Amat, RT, ART)); /* ART <- A*RT */ 93 bbot = btop; 94 PetscCall(VecDotBegin(RT, ART, &btop)); 95 96 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 97 PetscCall(VecNormBegin(RT, NORM_2, &dp)); /* dp <- || RT || */ 98 PetscCall(VecDotEnd(RT, ART, &btop)); 99 PetscCall(VecNormEnd(RT, NORM_2, &dp)); /* dp <- || RT || */ 100 KSPCheckNorm(ksp, dp); 101 } else if (ksp->normtype == KSP_NORM_NATURAL) { 102 PetscCall(VecDotEnd(RT, ART, &btop)); 103 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 104 } else if (ksp->normtype == KSP_NORM_NONE) { 105 PetscCall(VecDotEnd(RT, ART, &btop)); 106 dp = 0.0; /* meaningless value that is passed to monitor and convergence test */ 107 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 108 PetscCall(VecAXPY(R, ai, AP)); /* R <- R - ai*AP */ 109 PetscCall(VecNormBegin(R, NORM_2, &dp)); /* dp <- R'*R */ 110 PetscCall(VecDotEnd(RT, ART, &btop)); 111 PetscCall(VecNormEnd(R, NORM_2, &dp)); /* dp <- R'*R */ 112 KSPCheckNorm(ksp, dp); 113 } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSPNormType of %d not supported", (int)ksp->normtype); 114 if (PetscAbsScalar(btop) < 0.0) { 115 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 116 PetscCall(PetscInfo(ksp, "diverging due to indefinite or negative definite PC\n")); 117 break; 118 } 119 120 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 121 ksp->its++; 122 ksp->rnorm = dp; 123 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 124 125 PetscCall(KSPLogResidualHistory(ksp, dp)); 126 PetscCall(KSPMonitor(ksp, i + 1, dp)); 127 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP)); 128 if (ksp->reason) break; 129 130 bi = btop / bbot; 131 PetscCall(VecAYPX(P, bi, RT)); /* P <- RT + Bi P */ 132 PetscCall(VecAYPX(AP, bi, ART)); /* AP <- ART + Bi AP */ 133 i++; 134 } while (i < ksp->max_it); 135 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 /*MC 140 KSPCR - This code implements the (preconditioned) conjugate residuals method {cite}`hs:52` 141 142 Level: beginner 143 144 Notes: 145 The operator and the preconditioner must be symmetric for this method. 146 147 The preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE. 148 149 Support only for left preconditioning. 150 151 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG` 152 M*/ 153 PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp) 154 { 155 PetscFunctionBegin; 156 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3)); 157 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2)); 158 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2)); 159 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1)); 160 161 ksp->ops->setup = KSPSetUp_CR; 162 ksp->ops->solve = KSPSolve_CR; 163 ksp->ops->destroy = KSPDestroyDefault; 164 ksp->ops->buildsolution = KSPBuildSolutionDefault; 165 ksp->ops->buildresidual = KSPBuildResidualDefault; 166 ksp->ops->setfromoptions = NULL; 167 ksp->ops->view = NULL; 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170