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
```