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 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 References: 152 . * - Magnus R. Hestenes and Eduard Stiefel, Methods of Conjugate Gradients for Solving Linear Systems, 153 Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379 154 155 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG` 156 M*/ 157 PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp) 158 { 159 PetscFunctionBegin; 160 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3)); 161 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2)); 162 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2)); 163 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1)); 164 165 ksp->ops->setup = KSPSetUp_CR; 166 ksp->ops->solve = KSPSolve_CR; 167 ksp->ops->destroy = KSPDestroyDefault; 168 ksp->ops->buildsolution = KSPBuildSolutionDefault; 169 ksp->ops->buildresidual = KSPBuildResidualDefault; 170 ksp->ops->setfromoptions = NULL; 171 ksp->ops->view = NULL; 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174