Actual source code: cr.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/kspimpl.h>
6: static PetscErrorCode KSPSetUp_CR(KSP ksp)
7: {
11: if (ksp->pc_side == PC_RIGHT) SETERRQ(((PetscObject)ksp)->comm,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: KSPDefaultGetWork(ksp,6);
14: return(0);
15: }
19: static PetscErrorCode KSPSolve_CR(KSP ksp)
20: {
22: PetscInt i = 0;
23: MatStructure pflag;
24: PetscReal dp;
25: PetscScalar ai, bi;
26: PetscScalar apq,btop, bbot;
27: Vec X,B,R,RT,P,AP,ART,Q;
28: Mat Amat, Pmat;
31: X = ksp->vec_sol;
32: B = ksp->vec_rhs;
33: R = ksp->work[0];
34: RT = ksp->work[1];
35: P = ksp->work[2];
36: AP = ksp->work[3];
37: ART = ksp->work[4];
38: Q = ksp->work[5];
40: /* R is the true residual norm, RT is the preconditioned residual norm */
41: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
42: if (!ksp->guess_zero) {
43: KSP_MatMult(ksp,Amat,X,R); /* R <- A*X */
44: VecAYPX(R,-1.0,B); /* R <- B-R == B-A*X */
45: } else {
46: VecCopy(B,R); /* R <- B (X is 0) */
47: }
48: KSP_PCApply(ksp,R,P); /* P <- B*R */
49: KSP_MatMult(ksp,Amat,P,AP); /* AP <- A*P */
50: VecCopy(P,RT); /* RT <- P */
51: VecCopy(AP,ART); /* ART <- AP */
52: VecDotBegin(RT,ART,&btop); /* (RT,ART) */
53:
54: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
55: VecNormBegin(RT,NORM_2,&dp); /* dp <- RT'*RT */
56: VecDotEnd (RT,ART,&btop) ; /* (RT,ART) */
57: VecNormEnd (RT,NORM_2,&dp); /* dp <- RT'*RT */
58: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
59: VecNormBegin(R,NORM_2,&dp); /* dp <- R'*R */
60: VecDotEnd (RT,ART,&btop); /* (RT,ART) */
61: VecNormEnd (R,NORM_2,&dp); /* dp <- RT'*RT */
62: } else if (ksp->normtype == KSP_NORM_NATURAL) {
63: VecDotEnd (RT,ART,&btop) ; /* (RT,ART) */
64: dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */
65: }
66: if (PetscAbsScalar(btop) < 0.0) {
67: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
68: PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");
69: return(0);
70: }
72: ksp->its = 0;
73: KSPMonitor(ksp,0,dp);
74: PetscObjectTakeAccess(ksp);
75: ksp->rnorm = dp;
76: PetscObjectGrantAccess(ksp);
77: KSPLogResidualHistory(ksp,dp);
78: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
79: if (ksp->reason) return(0);
81: i = 0;
82: do {
83: KSP_PCApply(ksp,AP,Q);/* Q <- B* AP */
85: VecDot(AP,Q,&apq);
86: if (PetscAbsScalar(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: } else if (ksp->normtype == KSP_NORM_NATURAL) {
104: VecDotEnd(RT,ART,&btop);
105: dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */
106: } else if (ksp->normtype == KSP_NORM_NONE) {
107: VecDotEnd(RT,ART,&btop);
108: dp = 0.0;
109: } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
110: VecAXPY(R,ai,AP); /* R <- R - ai*AP */
111: VecNormBegin(R,NORM_2,&dp); /* dp <- R'*R */
112: VecDotEnd (RT,ART,&btop);
113: VecNormEnd (R,NORM_2,&dp); /* dp <- R'*R */
114: } else {
115: SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype);
116: }
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: PetscObjectTakeAccess(ksp);
124: ksp->its++;
125: ksp->rnorm = dp;
126: PetscObjectGrantAccess(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) {
139: ksp->reason = KSP_DIVERGED_ITS;
140: }
141: return(0);
142: }
145: /*MC
146: KSPCR - This code implements the (preconditioned) conjugate residuals method
148: Options Database Keys:
149: . see KSPSolve()
151: Level: beginner
153: Notes: The operator and the preconditioner must be symmetric for this method. The
154: preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE
155: Support only for left preconditioning.
157: References:
158: Methods of Conjugate Gradients for Solving Linear Systems, Magnus R. Hestenes and Eduard Stiefel,
159: Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379
160: pp. 409--436.
163: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG
164: M*/
165: EXTERN_C_BEGIN
168: PetscErrorCode KSPCreate_CR(KSP ksp)
169: {
173: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
174: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
175: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
177: ksp->ops->setup = KSPSetUp_CR;
178: ksp->ops->solve = KSPSolve_CR;
179: ksp->ops->destroy = KSPDefaultDestroy;
180: ksp->ops->buildsolution = KSPDefaultBuildSolution;
181: ksp->ops->buildresidual = KSPDefaultBuildResidual;
182: ksp->ops->setfromoptions = 0;
183: ksp->ops->view = 0;
184: return(0);
185: }
186: EXTERN_C_END