1 /*$Id: cr.c,v 1.64 2001/08/07 03:03:49 balay Exp $*/ 2 3 /* 4 This implements Preconditioned Conjugate Residuals. 5 */ 6 #include "src/ksp/ksp/kspimpl.h" 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "KSPSetUp_CR" 10 static int KSPSetUp_CR(KSP ksp) 11 { 12 int ierr; 13 14 PetscFunctionBegin; 15 if (ksp->pc_side == PC_RIGHT) {SETERRQ(2,"no right preconditioning for KSPCR");} 16 else if (ksp->pc_side == PC_SYMMETRIC) {SETERRQ(2,"no symmetric preconditioning for KSPCR");} 17 ierr = KSPDefaultGetWork(ksp,6);CHKERRQ(ierr); 18 PetscFunctionReturn(0); 19 } 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "KSPSolve_CR" 23 static int KSPSolve_CR(KSP ksp) 24 { 25 int i = 0, ierr; 26 MatStructure pflag; 27 PetscReal dp; 28 PetscScalar ai, bi; 29 PetscScalar apq,btop, bbot, mone = -1.0; 30 Vec X,B,R,RT,P,AP,ART,Q; 31 Mat Amat, Pmat; 32 33 PetscFunctionBegin; 34 35 X = ksp->vec_sol; 36 B = ksp->vec_rhs; 37 R = ksp->work[0]; 38 RT = ksp->work[1]; 39 P = ksp->work[2]; 40 AP = ksp->work[3]; 41 ART = ksp->work[4]; 42 Q = ksp->work[5]; 43 44 /* R is the true residual norm, RT is the preconditioned residual norm */ 45 ierr = PCGetOperators(ksp->B,&Amat,&Pmat,&pflag);CHKERRQ(ierr); 46 if (!ksp->guess_zero) { 47 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* R <- A*X */ 48 ierr = VecAYPX(&mone,B,R);CHKERRQ(ierr); /* R <- B-R == B-A*X */ 49 } else { 50 ierr = VecCopy(B,R);CHKERRQ(ierr); /* R <- B (X is 0) */ 51 } 52 ierr = KSP_PCApply(ksp,ksp->B,R,P);CHKERRQ(ierr); /* P <- B*R */ 53 ierr = KSP_MatMult(ksp,Amat,P,AP);CHKERRQ(ierr); /* AP <- A*P */ 54 ierr = VecCopy(P,RT);CHKERRQ(ierr); /* RT <- P */ 55 ierr = VecCopy(AP,ART);CHKERRQ(ierr); /* ART <- AP */ 56 ierr = VecDot(RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 57 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 58 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 59 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 60 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 61 } else if (ksp->normtype == KSP_NATURAL_NORM) { 62 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 63 } 64 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 65 ksp->its = 0; 66 if (ksp->reason) PetscFunctionReturn(0); 67 KSPMonitor(ksp,0,dp); 68 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 69 ksp->rnorm = dp; 70 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 71 KSPLogResidualHistory(ksp,dp); 72 73 i = 0; 74 do { 75 ierr = KSP_PCApply(ksp,ksp->B,AP,Q);CHKERRQ(ierr);/* Q <- B* AP */ 76 /* Step 3 */ 77 78 ierr = VecDot(AP,Q,&apq);CHKERRQ(ierr); 79 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 80 81 ierr = VecAXPY(&ai,P,X);CHKERRQ(ierr); /* X <- X + ai*P */ 82 ai = -ai; 83 ierr = VecAXPY(&ai,Q,RT);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 84 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr);/* ART <- A*RT */ 85 bbot = btop; 86 ierr = VecDot(RT,ART,&btop);CHKERRQ(ierr); 87 88 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 89 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 90 } else if (ksp->normtype == KSP_NATURAL_NORM) { 91 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 92 } else if (ksp->normtype == KSP_NO_NORM) { 93 dp = 0.0; 94 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 95 ierr = VecAXPY(&ai,AP,R);CHKERRQ(ierr); /* R <- R - ai*AP */ 96 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 97 } else { 98 SETERRQ1(1,"KSPNormType of %d not supported",(int)ksp->normtype); 99 } 100 101 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 102 ksp->its++; 103 ksp->rnorm = dp; 104 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 105 106 KSPLogResidualHistory(ksp,dp); 107 KSPMonitor(ksp,i+1,dp); 108 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 109 if (ksp->reason) break; 110 111 bi = btop/bbot; 112 ierr = VecAYPX(&bi,RT,P);CHKERRQ(ierr); /* P <- RT + Bi P */ 113 ierr = VecAYPX(&bi,ART,AP);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 114 i++; 115 } while (i<ksp->max_it); 116 if (i == ksp->max_it) { 117 ksp->reason = KSP_DIVERGED_ITS; 118 } 119 PetscFunctionReturn(0); 120 } 121 122 123 EXTERN_C_BEGIN 124 #undef __FUNCT__ 125 #define __FUNCT__ "KSPCreate_CR" 126 int KSPCreate_CR(KSP ksp) 127 { 128 PetscFunctionBegin; 129 ksp->pc_side = PC_LEFT; 130 ksp->ops->setup = KSPSetUp_CR; 131 ksp->ops->solve = KSPSolve_CR; 132 ksp->ops->destroy = KSPDefaultDestroy; 133 ksp->ops->buildsolution = KSPDefaultBuildSolution; 134 ksp->ops->buildresidual = KSPDefaultBuildResidual; 135 ksp->ops->setfromoptions = 0; 136 ksp->ops->view = 0; 137 PetscFunctionReturn(0); 138 } 139 EXTERN_C_END 140