1 2 #include <private/kspimpl.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "KSPSetUp_CR" 6 static PetscErrorCode KSPSetUp_CR(KSP ksp) 7 { 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 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 ierr = KSPDefaultGetWork(ksp,6);CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "KSPSolve_CR" 19 static PetscErrorCode KSPSolve_CR(KSP ksp) 20 { 21 PetscErrorCode ierr; 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; 29 30 PetscFunctionBegin; 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]; 39 40 /* R is the true residual norm, RT is the preconditioned residual norm */ 41 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr); 42 if (!ksp->guess_zero) { 43 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* R <- A*X */ 44 ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr); /* R <- B-R == B-A*X */ 45 } else { 46 ierr = VecCopy(B,R);CHKERRQ(ierr); /* R <- B (X is 0) */ 47 } 48 ierr = KSP_PCApply(ksp,R,P);CHKERRQ(ierr); /* P <- B*R */ 49 ierr = KSP_MatMult(ksp,Amat,P,AP);CHKERRQ(ierr); /* AP <- A*P */ 50 ierr = VecCopy(P,RT);CHKERRQ(ierr); /* RT <- P */ 51 ierr = VecCopy(AP,ART);CHKERRQ(ierr); /* ART <- AP */ 52 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 53 54 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 55 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 56 ierr = VecDotEnd (RT,ART,&btop) ;CHKERRQ(ierr); /* (RT,ART) */ 57 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 58 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 59 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 60 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 61 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 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 } 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 = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 75 ksp->rnorm = dp; 76 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 77 KSPLogResidualHistory(ksp,dp); 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 if (PetscAbsScalar(apq) <= 0.0) { 87 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 88 ierr = PetscInfo(ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 89 break; 90 } 91 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 92 93 ierr = VecAXPY(X,ai,P);CHKERRQ(ierr); /* X <- X + ai*P */ 94 ierr = VecAXPY(RT,-ai,Q);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 95 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr);/* ART <- A*RT */ 96 bbot = btop; 97 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); 98 99 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 100 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 101 ierr = VecDotEnd (RT,ART,&btop) ;CHKERRQ(ierr); 102 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 103 } else if (ksp->normtype == KSP_NORM_NATURAL) { 104 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 105 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 106 } else if (ksp->normtype == KSP_NORM_NONE) { 107 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 108 dp = 0.0; 109 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 110 ierr = VecAXPY(R,ai,AP);CHKERRQ(ierr); /* R <- R - ai*AP */ 111 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 112 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); 113 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* 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 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 120 break; 121 } 122 123 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 124 ksp->its++; 125 ksp->rnorm = dp; 126 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 127 128 KSPLogResidualHistory(ksp,dp); 129 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 130 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 131 if (ksp->reason) break; 132 133 bi = btop/bbot; 134 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 135 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* 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 PetscFunctionReturn(0); 142 } 143 144 145 /*MC 146 KSPCR - This code implements the (preconditioned) conjugate residuals method 147 148 Options Database Keys: 149 . see KSPSolve() 150 151 Level: beginner 152 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. 156 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. 161 162 163 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 164 M*/ 165 EXTERN_C_BEGIN 166 #undef __FUNCT__ 167 #define __FUNCT__ "KSPCreate_CR" 168 PetscErrorCode KSPCreate_CR(KSP ksp) 169 { 170 PetscErrorCode ierr; 171 172 PetscFunctionBegin; 173 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 174 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);CHKERRQ(ierr); 175 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);CHKERRQ(ierr); 176 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 PetscFunctionReturn(0); 185 } 186 EXTERN_C_END 187