1 #define PETSCKSP_DLL 2 3 #include "src/ksp/ksp/kspimpl.h" 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "KSPSetUp_CR" 7 static PetscErrorCode KSPSetUp_CR(KSP ksp) 8 { 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 if (ksp->pc_side == PC_RIGHT) {SETERRQ(PETSC_ERR_SUP,"no right preconditioning for KSPCR");} 13 else if (ksp->pc_side == PC_SYMMETRIC) {SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPCR");} 14 ierr = KSPDefaultGetWork(ksp,6);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "KSPSolve_CR" 20 static PetscErrorCode KSPSolve_CR(KSP ksp) 21 { 22 PetscErrorCode ierr; 23 PetscInt i = 0; 24 MatStructure pflag; 25 PetscReal dp; 26 PetscScalar ai, bi; 27 PetscScalar apq,btop, bbot, mone = -1.0; 28 Vec X,B,R,RT,P,AP,ART,Q; 29 Mat Amat, Pmat; 30 31 PetscFunctionBegin; 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(R,mone,B);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 (PetscAbsScalar(btop) < 0.0) { 55 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 56 ierr = PetscLogInfo((ksp,"KSPSolve_CR:diverging due to indefinite or negative definite matrix\n"));CHKERRQ(ierr); 57 PetscFunctionReturn(0); 58 } 59 60 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 61 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 62 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 63 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 64 } else if (ksp->normtype == KSP_NATURAL_NORM) { 65 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 66 } 67 ksp->its = 0; 68 KSPMonitor(ksp,0,dp); 69 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 70 ksp->rnorm = dp; 71 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 72 KSPLogResidualHistory(ksp,dp); 73 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 74 if (ksp->reason) PetscFunctionReturn(0); 75 76 i = 0; 77 do { 78 ierr = KSP_PCApply(ksp,AP,Q);CHKERRQ(ierr);/* Q <- B* AP */ 79 80 ierr = VecDot(AP,Q,&apq);CHKERRQ(ierr); 81 if (PetscAbsScalar(apq) <= 0.0) { 82 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 83 ierr = PetscLogInfo((ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n"));CHKERRQ(ierr); 84 break; 85 } 86 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 87 88 ierr = VecAXPY(X,ai,P);CHKERRQ(ierr); /* X <- X + ai*P */ 89 ai = -ai; 90 ierr = VecAXPY(RT,ai,Q);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 91 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr);/* ART <- A*RT */ 92 bbot = btop; 93 ierr = VecDot(RT,ART,&btop);CHKERRQ(ierr); 94 if (PetscAbsScalar(btop) < 0.0) { 95 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 96 ierr = PetscLogInfo((ksp,"KSPSolve_CR:diverging due to indefinite or negative definite matrix\n"));CHKERRQ(ierr); 97 break; 98 } 99 100 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 101 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 102 } else if (ksp->normtype == KSP_NATURAL_NORM) { 103 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 104 } else if (ksp->normtype == KSP_NO_NORM) { 105 dp = 0.0; 106 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 107 ierr = VecAXPY(R,ai,AP);CHKERRQ(ierr); /* R <- R - ai*AP */ 108 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 109 } else { 110 SETERRQ1(PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 111 } 112 113 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 114 ksp->its++; 115 ksp->rnorm = dp; 116 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 117 118 KSPLogResidualHistory(ksp,dp); 119 KSPMonitor(ksp,i+1,dp); 120 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 121 if (ksp->reason) break; 122 123 bi = btop/bbot; 124 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 125 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 126 i++; 127 } while (i<ksp->max_it); 128 if (i >= ksp->max_it) { 129 ksp->reason = KSP_DIVERGED_ITS; 130 } 131 PetscFunctionReturn(0); 132 } 133 134 135 /*MC 136 KSPCR - This code implements the (preconditioned) conjugate residuals method 137 138 Options Database Keys: 139 . see KSPSolve() 140 141 Level: beginner 142 143 Notes: The operator and the preconditioner must be symmetric for this method. The 144 preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE 145 146 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 147 M*/ 148 EXTERN_C_BEGIN 149 #undef __FUNCT__ 150 #define __FUNCT__ "KSPCreate_CR" 151 PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_CR(KSP ksp) 152 { 153 PetscFunctionBegin; 154 ksp->pc_side = PC_LEFT; 155 ksp->ops->setup = KSPSetUp_CR; 156 ksp->ops->solve = KSPSolve_CR; 157 ksp->ops->destroy = KSPDefaultDestroy; 158 ksp->ops->buildsolution = KSPDefaultBuildSolution; 159 ksp->ops->buildresidual = KSPDefaultBuildResidual; 160 ksp->ops->setfromoptions = 0; 161 ksp->ops->view = 0; 162 PetscFunctionReturn(0); 163 } 164 EXTERN_C_END 165