Actual source code: cr.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/kspimpl.h>
6: static PetscErrorCode KSPSetUp_CR(KSP ksp)
7: {
11: if (ksp->pc_side == PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no right preconditioning for KSPCR");
12: else if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"no symmetric preconditioning for KSPCR");
13: KSPSetWorkVecs(ksp,6);
14: return(0);
15: }
19: static PetscErrorCode KSPSolve_CR(KSP ksp)
20: {
22: PetscInt i = 0;
23: PetscReal dp;
24: PetscScalar ai, bi;
25: PetscScalar apq,btop, bbot;
26: Vec X,B,R,RT,P,AP,ART,Q;
27: Mat Amat, Pmat;
30: X = ksp->vec_sol;
31: B = ksp->vec_rhs;
32: R = ksp->work[0];
33: RT = ksp->work[1];
34: P = ksp->work[2];
35: AP = ksp->work[3];
36: ART = ksp->work[4];
37: Q = ksp->work[5];
39: /* R is the true residual norm, RT is the preconditioned residual norm */
40: PCGetOperators(ksp->pc,&Amat,&Pmat);
41: if (!ksp->guess_zero) {
42: KSP_MatMult(ksp,Amat,X,R); /* R <- A*X */
43: VecAYPX(R,-1.0,B); /* R <- B-R == B-A*X */
44: } else {
45: VecCopy(B,R); /* R <- B (X is 0) */
46: }
47: KSP_PCApply(ksp,R,P); /* P <- B*R */
48: KSP_MatMult(ksp,Amat,P,AP); /* AP <- A*P */
49: VecCopy(P,RT); /* RT <- P */
50: VecCopy(AP,ART); /* ART <- AP */
51: VecDotBegin(RT,ART,&btop); /* (RT,ART) */
53: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
54: VecNormBegin(RT,NORM_2,&dp); /* dp <- RT'*RT */
55: VecDotEnd (RT,ART,&btop); /* (RT,ART) */
56: VecNormEnd (RT,NORM_2,&dp); /* dp <- RT'*RT */
57: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
58: VecNormBegin(R,NORM_2,&dp); /* dp <- R'*R */
59: VecDotEnd (RT,ART,&btop); /* (RT,ART) */
60: VecNormEnd (R,NORM_2,&dp); /* dp <- RT'*RT */
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: }
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: if (PetscRealPart(apq) <= 0.0) {
86: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
87: PetscInfo(ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n");
88: break;
89: }
90: ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */
92: VecAXPY(X,ai,P); /* X <- X + ai*P */
93: VecAXPY(RT,-ai,Q); /* RT <- RT - ai*Q */
94: KSP_MatMult(ksp,Amat,RT,ART); /* ART <- A*RT */
95: bbot = btop;
96: VecDotBegin(RT,ART,&btop);
98: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
99: VecNormBegin(RT,NORM_2,&dp); /* dp <- || RT || */
100: VecDotEnd (RT,ART,&btop);
101: VecNormEnd (RT,NORM_2,&dp); /* dp <- || RT || */
102: } else if (ksp->normtype == KSP_NORM_NATURAL) {
103: VecDotEnd(RT,ART,&btop);
104: dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */
105: } else if (ksp->normtype == KSP_NORM_NONE) {
106: VecDotEnd(RT,ART,&btop);
107: dp = 0.0;
108: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
109: VecAXPY(R,ai,AP); /* R <- R - ai*AP */
110: VecNormBegin(R,NORM_2,&dp); /* dp <- R'*R */
111: VecDotEnd (RT,ART,&btop);
112: VecNormEnd (R,NORM_2,&dp); /* dp <- R'*R */
113: } else SETERRQ1(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: PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");
117: break;
118: }
120: PetscObjectSAWsTakeAccess((PetscObject)ksp);
121: ksp->its++;
122: ksp->rnorm = dp;
123: PetscObjectSAWsGrantAccess((PetscObject)ksp);
125: KSPLogResidualHistory(ksp,dp);
126: KSPMonitor(ksp,i+1,dp);
127: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
128: if (ksp->reason) break;
130: bi = btop/bbot;
131: VecAYPX(P,bi,RT); /* P <- RT + Bi P */
132: 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: return(0);
137: }
140: /*MC
141: KSPCR - This code implements the (preconditioned) conjugate residuals method
143: Options Database Keys:
144: . see KSPSolve()
146: Level: beginner
148: Notes: The operator and the preconditioner must be symmetric for this method. The
149: preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE
150: Support only for left preconditioning.
152: References:
153: Methods of Conjugate Gradients for Solving Linear Systems, Magnus R. Hestenes and Eduard Stiefel,
154: Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379
155: pp. 409--436.
158: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG
159: M*/
162: PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp)
163: {
167: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
168: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
169: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
171: ksp->ops->setup = KSPSetUp_CR;
172: ksp->ops->solve = KSPSolve_CR;
173: ksp->ops->destroy = KSPDestroyDefault;
174: ksp->ops->buildsolution = KSPBuildSolutionDefault;
175: ksp->ops->buildresidual = KSPBuildResidualDefault;
176: ksp->ops->setfromoptions = 0;
177: ksp->ops->view = 0;
178: return(0);
179: }