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