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; 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,-1.0,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 = PetscInfo(ksp,"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 = PetscInfo(ksp,"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 ierr = VecAXPY(RT,-ai,Q);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 90 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr);/* ART <- A*RT */ 91 bbot = btop; 92 ierr = VecDot(RT,ART,&btop);CHKERRQ(ierr); 93 if (PetscAbsScalar(btop) < 0.0) { 94 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 95 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");CHKERRQ(ierr); 96 break; 97 } 98 99 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 100 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 101 } else if (ksp->normtype == KSP_NATURAL_NORM) { 102 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 103 } else if (ksp->normtype == KSP_NO_NORM) { 104 dp = 0.0; 105 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 106 ierr = VecAXPY(R,-ai,AP);CHKERRQ(ierr); /* R <- R - ai*AP */ 107 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 108 } else { 109 SETERRQ1(PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 110 } 111 112 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 113 ksp->its++; 114 ksp->rnorm = dp; 115 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 116 117 KSPLogResidualHistory(ksp,dp); 118 KSPMonitor(ksp,i+1,dp); 119 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 120 if (ksp->reason) break; 121 122 bi = btop/bbot; 123 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 124 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 125 i++; 126 } while (i<ksp->max_it); 127 if (i >= ksp->max_it) { 128 ksp->reason = KSP_DIVERGED_ITS; 129 } 130 PetscFunctionReturn(0); 131 } 132 133 134 /*MC 135 KSPCR - This code implements the (preconditioned) conjugate residuals method 136 137 Options Database Keys: 138 . see KSPSolve() 139 140 Level: beginner 141 142 Notes: The operator and the preconditioner must be symmetric for this method. The 143 preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE 144 145 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 146 M*/ 147 EXTERN_C_BEGIN 148 #undef __FUNCT__ 149 #define __FUNCT__ "KSPCreate_CR" 150 PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_CR(KSP ksp) 151 { 152 PetscFunctionBegin; 153 ksp->pc_side = PC_LEFT; 154 ksp->ops->setup = KSPSetUp_CR; 155 ksp->ops->solve = KSPSolve_CR; 156 ksp->ops->destroy = KSPDefaultDestroy; 157 ksp->ops->buildsolution = KSPDefaultBuildSolution; 158 ksp->ops->buildresidual = KSPDefaultBuildResidual; 159 ksp->ops->setfromoptions = 0; 160 ksp->ops->view = 0; 161 PetscFunctionReturn(0); 162 } 163 EXTERN_C_END 164