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