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