1 2 #include <petsc-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(PetscObjectComm((PetscObject)ksp),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 = KSPSetWorkVecs(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 PetscReal dp; 24 PetscScalar ai, bi; 25 PetscScalar apq,btop, bbot; 26 Vec X,B,R,RT,P,AP,ART,Q; 27 Mat Amat, Pmat; 28 29 PetscFunctionBegin; 30 X = ksp->vec_sol; 31 B = ksp->vec_rhs; 32 R = ksp->work[0]; 33 RT = ksp->work[1]; 34 P = ksp->work[2]; 35 AP = ksp->work[3]; 36 ART = ksp->work[4]; 37 Q = ksp->work[5]; 38 39 /* R is the true residual norm, RT is the preconditioned residual norm */ 40 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr); 41 if (!ksp->guess_zero) { 42 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* R <- A*X */ 43 ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr); /* R <- B-R == B-A*X */ 44 } else { 45 ierr = VecCopy(B,R);CHKERRQ(ierr); /* R <- B (X is 0) */ 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 } 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 } else if (ksp->normtype == KSP_NORM_NATURAL) { 62 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 63 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 64 } 65 if (PetscAbsScalar(btop) < 0.0) { 66 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 67 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");CHKERRQ(ierr); 68 PetscFunctionReturn(0); 69 } 70 71 ksp->its = 0; 72 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr); 73 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 74 ksp->rnorm = dp; 75 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 76 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 77 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 78 if (ksp->reason) PetscFunctionReturn(0); 79 80 i = 0; 81 do { 82 ierr = KSP_PCApply(ksp,AP,Q);CHKERRQ(ierr); /* Q <- B* AP */ 83 84 ierr = VecDot(AP,Q,&apq);CHKERRQ(ierr); 85 if (PetscRealPart(apq) <= 0.0) { 86 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 87 ierr = PetscInfo(ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 88 break; 89 } 90 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 91 92 ierr = VecAXPY(X,ai,P);CHKERRQ(ierr); /* X <- X + ai*P */ 93 ierr = VecAXPY(RT,-ai,Q);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 94 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr); /* ART <- A*RT */ 95 bbot = btop; 96 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); 97 98 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 99 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 100 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); 101 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 102 } else if (ksp->normtype == KSP_NORM_NATURAL) { 103 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 104 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 105 } else if (ksp->normtype == KSP_NORM_NONE) { 106 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 107 dp = 0.0; 108 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 109 ierr = VecAXPY(R,ai,AP);CHKERRQ(ierr); /* R <- R - ai*AP */ 110 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 111 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); 112 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 113 } else SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 114 if (PetscAbsScalar(btop) < 0.0) { 115 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 116 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 117 break; 118 } 119 120 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 121 ksp->its++; 122 ksp->rnorm = dp; 123 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 124 125 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 126 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 127 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 128 if (ksp->reason) break; 129 130 bi = btop/bbot; 131 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 132 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 133 i++; 134 } while (i<ksp->max_it); 135 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 136 PetscFunctionReturn(0); 137 } 138 139 140 /*MC 141 KSPCR - This code implements the (preconditioned) conjugate residuals method 142 143 Options Database Keys: 144 . see KSPSolve() 145 146 Level: beginner 147 148 Notes: The operator and the preconditioner must be symmetric for this method. The 149 preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE 150 Support only for left preconditioning. 151 152 References: 153 Methods of Conjugate Gradients for Solving Linear Systems, Magnus R. Hestenes and Eduard Stiefel, 154 Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379 155 pp. 409--436. 156 157 158 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 159 M*/ 160 #undef __FUNCT__ 161 #define __FUNCT__ "KSPCreate_CR" 162 PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr); 168 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 169 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr); 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 = 0; 177 ksp->ops->view = 0; 178 PetscFunctionReturn(0); 179 } 180