1 2 #include "src/ksp/ksp/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(2,"no right preconditioning for KSPCR");} 12 else if (ksp->pc_side == PC_SYMMETRIC) {SETERRQ(2,"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, mone = -1.0; 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(&mone,B,R);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 = VecDot(RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 53 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 54 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 55 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 56 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 57 } else if (ksp->normtype == KSP_NATURAL_NORM) { 58 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 59 } 60 ksp->its = 0; 61 KSPMonitor(ksp,0,dp); 62 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 63 ksp->rnorm = dp; 64 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 65 KSPLogResidualHistory(ksp,dp); 66 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 67 if (ksp->reason) PetscFunctionReturn(0); 68 69 i = 0; 70 do { 71 ierr = KSP_PCApply(ksp,AP,Q);CHKERRQ(ierr);/* Q <- B* AP */ 72 /* Step 3 */ 73 74 ierr = VecDot(AP,Q,&apq);CHKERRQ(ierr); 75 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 76 77 ierr = VecAXPY(&ai,P,X);CHKERRQ(ierr); /* X <- X + ai*P */ 78 ai = -ai; 79 ierr = VecAXPY(&ai,Q,RT);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 80 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr);/* ART <- A*RT */ 81 bbot = btop; 82 ierr = VecDot(RT,ART,&btop);CHKERRQ(ierr); 83 84 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 85 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 86 } else if (ksp->normtype == KSP_NATURAL_NORM) { 87 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 88 } else if (ksp->normtype == KSP_NO_NORM) { 89 dp = 0.0; 90 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 91 ierr = VecAXPY(&ai,AP,R);CHKERRQ(ierr); /* R <- R - ai*AP */ 92 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 93 } else { 94 SETERRQ1(PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 95 } 96 97 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 98 ksp->its++; 99 ksp->rnorm = dp; 100 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 101 102 KSPLogResidualHistory(ksp,dp); 103 KSPMonitor(ksp,i+1,dp); 104 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 105 if (ksp->reason) break; 106 107 bi = btop/bbot; 108 ierr = VecAYPX(&bi,RT,P);CHKERRQ(ierr); /* P <- RT + Bi P */ 109 ierr = VecAYPX(&bi,ART,AP);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 110 i++; 111 } while (i<ksp->max_it); 112 if (i >= ksp->max_it) { 113 ksp->reason = KSP_DIVERGED_ITS; 114 } 115 PetscFunctionReturn(0); 116 } 117 118 119 /*MC 120 KSPCR - This code implements the (preconditioned) conjugate residuals method 121 122 Options Database Keys: 123 . see KSPSolve() 124 125 Level: beginner 126 127 Notes: The operator and the preconditioner must be symmetric for this method 128 129 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 130 M*/ 131 EXTERN_C_BEGIN 132 #undef __FUNCT__ 133 #define __FUNCT__ "KSPCreate_CR" 134 PetscErrorCode KSPCreate_CR(KSP ksp) 135 { 136 PetscFunctionBegin; 137 ksp->pc_side = PC_LEFT; 138 ksp->ops->setup = KSPSetUp_CR; 139 ksp->ops->solve = KSPSolve_CR; 140 ksp->ops->destroy = KSPDefaultDestroy; 141 ksp->ops->buildsolution = KSPDefaultBuildSolution; 142 ksp->ops->buildresidual = KSPDefaultBuildResidual; 143 ksp->ops->setfromoptions = 0; 144 ksp->ops->view = 0; 145 PetscFunctionReturn(0); 146 } 147 EXTERN_C_END 148