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