1 #define PETSCKSP_DLL 2 3 #include "include/private/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_NORM_PRECONDITIONED) { 61 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 62 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 63 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 64 } else if (ksp->normtype == KSP_NORM_NATURAL) { 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_NORM_PRECONDITIONED) { 100 ierr = VecNorm(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 101 } else if (ksp->normtype == KSP_NORM_NATURAL) { 102 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 103 } else if (ksp->normtype == KSP_NORM_NO) { 104 dp = 0.0; 105 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 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 146 References: 147 Methods of Conjugate Gradients for Solving Linear Systems, Magnus R. Hestenes and Eduard Stiefel, 148 Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379 149 pp. 409--436. 150 151 152 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 153 M*/ 154 EXTERN_C_BEGIN 155 #undef __FUNCT__ 156 #define __FUNCT__ "KSPCreate_CR" 157 PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_CR(KSP ksp) 158 { 159 PetscFunctionBegin; 160 ksp->pc_side = PC_LEFT; 161 ksp->ops->setup = KSPSetUp_CR; 162 ksp->ops->solve = KSPSolve_CR; 163 ksp->ops->destroy = KSPDefaultDestroy; 164 ksp->ops->buildsolution = KSPDefaultBuildSolution; 165 ksp->ops->buildresidual = KSPDefaultBuildResidual; 166 ksp->ops->setfromoptions = 0; 167 ksp->ops->view = 0; 168 PetscFunctionReturn(0); 169 } 170 EXTERN_C_END 171