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