Actual source code: cr.c
2: #include <petsc/private/kspimpl.h>
4: static PetscErrorCode KSPSetUp_CR(KSP ksp)
5: {
8: KSPSetWorkVecs(ksp,6);
9: return 0;
10: }
12: static PetscErrorCode KSPSolve_CR(KSP ksp)
13: {
14: PetscInt i = 0;
15: PetscReal dp;
16: PetscScalar ai, bi;
17: PetscScalar apq,btop, bbot;
18: Vec X,B,R,RT,P,AP,ART,Q;
19: Mat Amat, Pmat;
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];
30: /* R is the true residual norm, RT is the preconditioned residual norm */
31: PCGetOperators(ksp->pc,&Amat,&Pmat);
32: if (!ksp->guess_zero) {
33: KSP_MatMult(ksp,Amat,X,R); /* R <- A*X */
34: VecAYPX(R,-1.0,B); /* R <- B-R == B-A*X */
35: } else {
36: 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) {
40: VecSetInf(R);
41: }
42: KSP_PCApply(ksp,R,P); /* P <- B*R */
43: KSP_MatMult(ksp,Amat,P,AP); /* AP <- A*P */
44: VecCopy(P,RT); /* RT <- P */
45: VecCopy(AP,ART); /* ART <- AP */
46: VecDotBegin(RT,ART,&btop); /* (RT,ART) */
48: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
49: VecNormBegin(RT,NORM_2,&dp); /* dp <- RT'*RT */
50: VecDotEnd (RT,ART,&btop); /* (RT,ART) */
51: 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: VecDotEnd (RT,ART,&btop); /* (RT,ART) */
56: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
57: VecNormBegin(R,NORM_2,&dp); /* dp <- R'*R */
58: VecDotEnd (RT,ART,&btop); /* (RT,ART) */
59: VecNormEnd (R,NORM_2,&dp); /* dp <- RT'*RT */
60: KSPCheckNorm(ksp,dp);
61: } else if (ksp->normtype == KSP_NORM_NATURAL) {
62: 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: PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");
68: return 0;
69: }
71: ksp->its = 0;
72: KSPMonitor(ksp,0,dp);
73: PetscObjectSAWsTakeAccess((PetscObject)ksp);
74: ksp->rnorm = dp;
75: PetscObjectSAWsGrantAccess((PetscObject)ksp);
76: KSPLogResidualHistory(ksp,dp);
77: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
78: if (ksp->reason) return 0;
80: i = 0;
81: do {
82: KSP_PCApply(ksp,AP,Q); /* Q <- B* AP */
84: VecDot(AP,Q,&apq);
85: KSPCheckDot(ksp,apq);
86: if (PetscRealPart(apq) <= 0.0) {
87: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
88: 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) */
93: VecAXPY(X,ai,P); /* X <- X + ai*P */
94: VecAXPY(RT,-ai,Q); /* RT <- RT - ai*Q */
95: KSP_MatMult(ksp,Amat,RT,ART); /* ART <- A*RT */
96: bbot = btop;
97: VecDotBegin(RT,ART,&btop);
99: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
100: VecNormBegin(RT,NORM_2,&dp); /* dp <- || RT || */
101: VecDotEnd (RT,ART,&btop);
102: VecNormEnd (RT,NORM_2,&dp); /* dp <- || RT || */
103: KSPCheckNorm(ksp,dp);
104: } else if (ksp->normtype == KSP_NORM_NATURAL) {
105: VecDotEnd(RT,ART,&btop);
106: dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */
107: } else if (ksp->normtype == KSP_NORM_NONE) {
108: 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: VecAXPY(R,ai,AP); /* R <- R - ai*AP */
112: VecNormBegin(R,NORM_2,&dp); /* dp <- R'*R */
113: VecDotEnd (RT,ART,&btop);
114: 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: PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");
120: break;
121: }
123: PetscObjectSAWsTakeAccess((PetscObject)ksp);
124: ksp->its++;
125: ksp->rnorm = dp;
126: PetscObjectSAWsGrantAccess((PetscObject)ksp);
128: KSPLogResidualHistory(ksp,dp);
129: KSPMonitor(ksp,i+1,dp);
130: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
131: if (ksp->reason) break;
133: bi = btop/bbot;
134: VecAYPX(P,bi,RT); /* P <- RT + Bi P */
135: 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: return 0;
140: }
142: /*MC
143: KSPCR - This code implements the (preconditioned) conjugate residuals method
145: Options Database Keys:
146: see KSPSolve()
148: Level: beginner
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.
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
159: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG
160: M*/
161: PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp)
162: {
163: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
164: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
165: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
166: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
168: ksp->ops->setup = KSPSetUp_CR;
169: ksp->ops->solve = KSPSolve_CR;
170: ksp->ops->destroy = KSPDestroyDefault;
171: ksp->ops->buildsolution = KSPBuildSolutionDefault;
172: ksp->ops->buildresidual = KSPBuildResidualDefault;
173: ksp->ops->setfromoptions = NULL;
174: ksp->ops->view = NULL;
175: return 0;
176: }