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