1 2 #include <petsc-private/kspimpl.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "KSPSetUp_CR" 6 static PetscErrorCode KSPSetUp_CR(KSP ksp) 7 { 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 if (ksp->pc_side == PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"no right preconditioning for KSPCR"); 12 else if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"no symmetric preconditioning for KSPCR"); 13 ierr = KSPSetWorkVecs(ksp,6);CHKERRQ(ierr); 14 PetscFunctionReturn(0); 15 } 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "KSPSolve_CR" 19 static PetscErrorCode KSPSolve_CR(KSP ksp) 20 { 21 PetscErrorCode ierr; 22 PetscInt i = 0; 23 MatStructure pflag; 24 PetscReal dp; 25 PetscScalar ai, bi; 26 PetscScalar apq,btop, bbot; 27 Vec X,B,R,RT,P,AP,ART,Q; 28 Mat Amat, Pmat; 29 30 PetscFunctionBegin; 31 X = ksp->vec_sol; 32 B = ksp->vec_rhs; 33 R = ksp->work[0]; 34 RT = ksp->work[1]; 35 P = ksp->work[2]; 36 AP = ksp->work[3]; 37 ART = ksp->work[4]; 38 Q = ksp->work[5]; 39 40 /* R is the true residual norm, RT is the preconditioned residual norm */ 41 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr); 42 if (!ksp->guess_zero) { 43 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* R <- A*X */ 44 ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr); /* R <- B-R == B-A*X */ 45 } else { 46 ierr = VecCopy(B,R);CHKERRQ(ierr); /* R <- B (X is 0) */ 47 } 48 ierr = KSP_PCApply(ksp,R,P);CHKERRQ(ierr); /* P <- B*R */ 49 ierr = KSP_MatMult(ksp,Amat,P,AP);CHKERRQ(ierr); /* AP <- A*P */ 50 ierr = VecCopy(P,RT);CHKERRQ(ierr); /* RT <- P */ 51 ierr = VecCopy(AP,ART);CHKERRQ(ierr); /* ART <- AP */ 52 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 53 54 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 55 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 56 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 57 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 58 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 59 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 60 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 61 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 62 } else if (ksp->normtype == KSP_NORM_NATURAL) { 63 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 64 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 65 } 66 if (PetscAbsScalar(btop) < 0.0) { 67 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 68 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite matrix\n");CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 ksp->its = 0; 73 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr); 74 ierr = PetscObjectAMSTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 75 ksp->rnorm = dp; 76 ierr = PetscObjectAMSGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 77 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 78 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 79 if (ksp->reason) PetscFunctionReturn(0); 80 81 i = 0; 82 do { 83 ierr = KSP_PCApply(ksp,AP,Q);CHKERRQ(ierr); /* Q <- B* AP */ 84 85 ierr = VecDot(AP,Q,&apq);CHKERRQ(ierr); 86 if (PetscRealPart(apq) <= 0.0) { 87 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 88 ierr = PetscInfo(ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 89 break; 90 } 91 ai = btop/apq; /* ai = (RT,ART)/(AP,Q) */ 92 93 ierr = VecAXPY(X,ai,P);CHKERRQ(ierr); /* X <- X + ai*P */ 94 ierr = VecAXPY(RT,-ai,Q);CHKERRQ(ierr); /* RT <- RT - ai*Q */ 95 ierr = KSP_MatMult(ksp,Amat,RT,ART);CHKERRQ(ierr); /* ART <- A*RT */ 96 bbot = btop; 97 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); 98 99 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 100 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 101 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); 102 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- || RT || */ 103 } else if (ksp->normtype == KSP_NORM_NATURAL) { 104 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 105 dp = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 106 } else if (ksp->normtype == KSP_NORM_NONE) { 107 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 108 dp = 0.0; 109 } else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) { 110 ierr = VecAXPY(R,ai,AP);CHKERRQ(ierr); /* R <- R - ai*AP */ 111 ierr = VecNormBegin(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 112 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); 113 ierr = VecNormEnd (R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- R'*R */ 114 } else SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 115 if (PetscAbsScalar(btop) < 0.0) { 116 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 117 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 118 break; 119 } 120 121 ierr = PetscObjectAMSTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 122 ksp->its++; 123 ksp->rnorm = dp; 124 ierr = PetscObjectAMSGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 125 126 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 127 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 128 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 129 if (ksp->reason) break; 130 131 bi = btop/bbot; 132 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 133 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 134 i++; 135 } while (i<ksp->max_it); 136 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 137 PetscFunctionReturn(0); 138 } 139 140 141 /*MC 142 KSPCR - This code implements the (preconditioned) conjugate residuals method 143 144 Options Database Keys: 145 . see KSPSolve() 146 147 Level: beginner 148 149 Notes: The operator and the preconditioner must be symmetric for this method. The 150 preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE 151 Support only for left preconditioning. 152 153 References: 154 Methods of Conjugate Gradients for Solving Linear Systems, Magnus R. Hestenes and Eduard Stiefel, 155 Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379 156 pp. 409--436. 157 158 159 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 160 M*/ 161 #undef __FUNCT__ 162 #define __FUNCT__ "KSPCreate_CR" 163 PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp) 164 { 165 PetscErrorCode ierr; 166 167 PetscFunctionBegin; 168 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 169 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);CHKERRQ(ierr); 170 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);CHKERRQ(ierr); 171 172 ksp->ops->setup = KSPSetUp_CR; 173 ksp->ops->solve = KSPSolve_CR; 174 ksp->ops->destroy = KSPDestroyDefault; 175 ksp->ops->buildsolution = KSPBuildSolutionDefault; 176 ksp->ops->buildresidual = KSPBuildResidualDefault; 177 ksp->ops->setfromoptions = 0; 178 ksp->ops->view = 0; 179 PetscFunctionReturn(0); 180 } 181