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 PetscReal dp; 24 PetscScalar ai, bi; 25 PetscScalar apq,btop, bbot; 26 Vec X,B,R,RT,P,AP,ART,Q; 27 Mat Amat, Pmat; 28 29 PetscFunctionBegin; 30 X = ksp->vec_sol; 31 B = ksp->vec_rhs; 32 R = ksp->work[0]; 33 RT = ksp->work[1]; 34 P = ksp->work[2]; 35 AP = ksp->work[3]; 36 ART = ksp->work[4]; 37 Q = ksp->work[5]; 38 39 /* R is the true residual norm, RT is the preconditioned residual norm */ 40 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr); 41 if (!ksp->guess_zero) { 42 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* R <- A*X */ 43 ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr); /* R <- B-R == B-A*X */ 44 } else { 45 ierr = VecCopy(B,R);CHKERRQ(ierr); /* R <- B (X is 0) */ 46 } 47 ierr = KSP_PCApply(ksp,R,P);CHKERRQ(ierr); /* P <- B*R */ 48 ierr = KSP_MatMult(ksp,Amat,P,AP);CHKERRQ(ierr); /* AP <- A*P */ 49 ierr = VecCopy(P,RT);CHKERRQ(ierr); /* RT <- P */ 50 ierr = VecCopy(AP,ART);CHKERRQ(ierr); /* ART <- AP */ 51 ierr = VecDotBegin(RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 52 53 if (ksp->normtype == KSP_NORM_PRECONDITIONED) { 54 ierr = VecNormBegin(RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 55 ierr = VecDotEnd (RT,ART,&btop);CHKERRQ(ierr); /* (RT,ART) */ 56 ierr = VecNormEnd (RT,NORM_2,&dp);CHKERRQ(ierr); /* dp <- RT'*RT */ 57 } else if (ksp->normtype == KSP_NORM_NONE) { 58 dp = 0.0; /* meaningless value that is passed to monitor and convergence test */ 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 = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 66 } else SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 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 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr); 75 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 76 ksp->rnorm = dp; 77 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 78 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 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 (PetscRealPart(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 = PetscSqrtReal(PetscAbsScalar(btop)); /* dp = sqrt(R,AR) */ 107 } else if (ksp->normtype == KSP_NORM_NONE) { 108 ierr = VecDotEnd(RT,ART,&btop);CHKERRQ(ierr); 109 dp = 0.0; /* meaningless value that is passed to monitor and convergence test */ 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 SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype); 116 if (PetscAbsScalar(btop) < 0.0) { 117 ksp->reason = KSP_DIVERGED_INDEFINITE_MAT; 118 ierr = PetscInfo(ksp,"diverging due to indefinite or negative definite PC\n");CHKERRQ(ierr); 119 break; 120 } 121 122 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 123 ksp->its++; 124 ksp->rnorm = dp; 125 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 126 127 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 128 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 129 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 130 if (ksp->reason) break; 131 132 bi = btop/bbot; 133 ierr = VecAYPX(P,bi,RT);CHKERRQ(ierr); /* P <- RT + Bi P */ 134 ierr = VecAYPX(AP,bi,ART);CHKERRQ(ierr); /* AP <- ART + Bi AP */ 135 i++; 136 } while (i<ksp->max_it); 137 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 138 PetscFunctionReturn(0); 139 } 140 141 142 /*MC 143 KSPCR - This code implements the (preconditioned) conjugate residuals method 144 145 Options Database Keys: 146 . see KSPSolve() 147 148 Level: beginner 149 150 Notes: The operator and the preconditioner must be symmetric for this method. The 151 preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE. 152 Support only for left preconditioning. 153 154 References: 155 . 1. - Magnus R. Hestenes and Eduard Stiefel, Methods of Conjugate Gradients for Solving Linear Systems, 156 Journal of Research of the National Bureau of Standards Vol. 49, No. 6, December 1952 Research Paper 2379 157 158 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG 159 M*/ 160 #undef __FUNCT__ 161 #define __FUNCT__ "KSPCreate_CR" 162 PETSC_EXTERN PetscErrorCode KSPCreate_CR(KSP ksp) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr); 168 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 169 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr); 170 171 ksp->ops->setup = KSPSetUp_CR; 172 ksp->ops->solve = KSPSolve_CR; 173 ksp->ops->destroy = KSPDestroyDefault; 174 ksp->ops->buildsolution = KSPBuildSolutionDefault; 175 ksp->ops->buildresidual = KSPBuildResidualDefault; 176 ksp->ops->setfromoptions = 0; 177 ksp->ops->view = 0; 178 PetscFunctionReturn(0); 179 } 180