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