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 = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 54 55 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 56 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 57 ierr = VecDotEnd (RT,ART,&btop) ;CHKERRQ(ierr); /* (RT,ART) */ 58 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 59 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 60 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 61 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 62 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 63 } else if (ksp->normtype == KSP_NORM_NATURAL) { 64 ierr = VecDotEnd (RT,ART,&btop) ;CHKERRQ(ierr); /* (RT,ART) */ 65 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 66 } 67 if (PetscAbsScalar(btop) < 0.0) { 68 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 69 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 ksp->its = 0; 74 KSPMonitor(ksp,0,dp); 75 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 76 ksp->rnorm = dp; 77 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 78 KSPLogResidualHistory(ksp,dp); 79 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 80 if (ksp->reason) PetscFunctionReturn(0); 81 82 i = 0; 83 do { 84 ierr = KSP_PCApply(ksp,AP,Q);CHKERRQ(ierr);/* Q <- B* AP */ 85 86 ierr = VecDot(AP,Q,&apq);CHKERRQ(ierr); 87 if (PetscAbsScalar(apq) <= 0.0) { 88 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 89 ierr = PetscInfo(ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 90 break; 91 } 92 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 93 94 ierr = VecAXPY(X,ai,P);CHKERRQ(ierr); /* X <- X + ai*P */ 95 ierr = VecAXPY(RT,-ai,Q);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 96 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr);/* ART <- A*RT */ 97 bbot = btop; 98 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); 99 100 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 101 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 102 ierr = VecDotEnd (RT,ART,&btop) ;CHKERRQ(ierr); 103 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 104 } else if (ksp->normtype == KSP_NORM_NATURAL) { 105 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 106 dp = sqrt(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 107 } else if (ksp->normtype == KSP_NORM_NO) { 108 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 109 dp = 0.0; 110 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 111 ierr = VecAXPY(R,ai,AP);CHKERRQ(ierr); /* R <- R - ai*AP */ 112 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 113 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); 114 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 115 } else { 116 SETERRQ1(PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 117 } 118 if (PetscAbsScalar(btop) < 0.0) { 119 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 120 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 121 break; 122 } 123 124 ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr); 125 ksp->its++; 126 ksp->rnorm = dp; 127 ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr); 128 129 KSPLogResidualHistory(ksp,dp); 130 KSPMonitor(ksp,i+1,dp); 131 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 132 if (ksp->reason) break; 133 134 bi = btop/bbot; 135 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 136 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 137 i++; 138 } while (i<ksp->max_it); 139 if (i >= ksp->max_it) { 140 ksp->reason = KSP_DIVERGED_ITS; 141 } 142 PetscFunctionReturn(0); 143 } 144 145 146 /*MC 147 KSPCR - This code implements the (preconditioned) conjugate residuals method 148 149 Options Database Keys: 150 . see KSPSolve() 151 152 Level: beginner 153 154 Notes: The operator and the preconditioner must be symmetric for this method. The 155 preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE 156 157 158 References: 159 Methods of Conjugate Gradients for Solving Linear Systems, Magnus R. Hestenes and Eduard Stiefel, 160 Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379 161 pp. 409--436. 162 163 164 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 165 M*/ 166 EXTERN_C_BEGIN 167 #undef __FUNCT__ 168 #define __FUNCT__ "KSPCreate_CR" 169 PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_CR(KSP ksp) 170 { 171 PetscFunctionBegin; 172 ksp->pc_side = PC_LEFT; 173 ksp->ops->setup = KSPSetUp_CR; 174 ksp->ops->solve = KSPSolve_CR; 175 ksp->ops->destroy = KSPDefaultDestroy; 176 ksp->ops->buildsolution = KSPDefaultBuildSolution; 177 ksp->ops->buildresidual = KSPDefaultBuildResidual; 178 ksp->ops->setfromoptions = 0; 179 ksp->ops->view = 0; 180 PetscFunctionReturn(0); 181 } 182 EXTERN_C_END 183