1 2 /* 3 This file implements the conjugate gradient method in PETSc as part of 4 KSP. You can use this as a starting point for implementing your own 5 Krylov method that is not provided with PETSc. 6 7 The following basic routines are required for each Krylov method. 8 KSPCreate_XXX() - Creates the Krylov context 9 KSPSetFromOptions_XXX() - Sets runtime options 10 KSPSolve_XXX() - Runs the Krylov method 11 KSPDestroy_XXX() - Destroys the Krylov context, freeing all 12 memory it needed 13 Here the "_XXX" denotes a particular implementation, in this case 14 we use _CG (e.g. KSPCreate_CG, KSPDestroy_CG). These routines are 15 are actually called vai the common user interface routines 16 KSPSetType(), KSPSetFromOptions(), KSPSolve(), and KSPDestroy() so the 17 application code interface remains identical for all preconditioners. 18 19 Other basic routines for the KSP objects include 20 KSPSetUp_XXX() 21 KSPView_XXX() - Prints details of solver being used. 22 23 Detailed notes: 24 By default, this code implements the CG (Conjugate Gradient) method, 25 which is valid for real symmetric (and complex Hermitian) positive 26 definite matrices. Note that for the complex Hermitian case, the 27 VecDot() arguments within the code MUST remain in the order given 28 for correct computation of inner products. 29 30 Reference: Hestenes and Steifel, 1952. 31 32 By switching to the indefinite vector inner product, VecTDot(), the 33 same code is used for the complex symmetric case as well. The user 34 must call KSPCGSetType(ksp,KSP_CG_SYMMETRIC) or use the option 35 -ksp_cg_symmetric to invoke this variant for the complex case. 36 Note, however, that the complex symmetric code is NOT valid for 37 all such matrices ... and thus we don't recommend using this method. 38 */ 39 /* 40 cgctx.h defines the simple data structured used to store information 41 related to the type of matrix (e.g. complex symmetric) being solved and 42 data used during the optional Lanczo process used to compute eigenvalues 43 */ 44 #include "src/ksp/ksp/impls/cg/cgctx.h" /*I "petscksp.h" I*/ 45 EXTERN PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP,PetscReal *,PetscReal *); 46 EXTERN PetscErrorCode KSPComputeEigenvalues_CG(KSP,int,PetscReal *,PetscReal *,int *); 47 48 /* 49 KSPSetUp_CG - Sets up the workspace needed by the CG method. 50 51 This is called once, usually automatically by KSPSolve() or KSPSetUp() 52 but can be called directly by KSPSetUp() 53 */ 54 #undef __FUNCT__ 55 #define __FUNCT__ "KSPSetUp_CG" 56 PetscErrorCode KSPSetUp_CG(KSP ksp) 57 { 58 KSP_CG *cgP = (KSP_CG*)ksp->data; 59 PetscErrorCode ierr; 60 int maxit = ksp->max_it; 61 62 PetscFunctionBegin; 63 /* 64 This implementation of CG only handles left preconditioning 65 so generate an error otherwise. 66 */ 67 if (ksp->pc_side == PC_RIGHT) { 68 SETERRQ(2,"No right preconditioning for KSPCG"); 69 } else if (ksp->pc_side == PC_SYMMETRIC) { 70 SETERRQ(2,"No symmetric preconditioning for KSPCG"); 71 } 72 73 /* get work vectors needed by CG */ 74 ierr = KSPDefaultGetWork(ksp,3);CHKERRQ(ierr); 75 76 /* 77 If user requested computations of eigenvalues then allocate work 78 work space needed 79 */ 80 if (ksp->calc_sings) { 81 /* get space to store tridiagonal matrix for Lanczos */ 82 ierr = PetscMalloc(2*(maxit+1)*sizeof(PetscScalar),&cgP->e);CHKERRQ(ierr); 83 PetscLogObjectMemory(ksp,2*(maxit+1)*sizeof(PetscScalar)); 84 cgP->d = cgP->e + maxit + 1; 85 ierr = PetscMalloc(2*(maxit+1)*sizeof(PetscReal),&cgP->ee);CHKERRQ(ierr); 86 PetscLogObjectMemory(ksp,2*(maxit+1)*sizeof(PetscScalar)); 87 cgP->dd = cgP->ee + maxit + 1; 88 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG; 89 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG; 90 } 91 PetscFunctionReturn(0); 92 } 93 94 /* 95 KSPSolve_CG - This routine actually applies the conjugate gradient 96 method 97 98 Input Parameter: 99 . ksp - the Krylov space object that was set to use conjugate gradient, by, for 100 example, KSPCreate(MPI_Comm,KSP *ksp); KSPSetType(ksp,KSPCG); 101 102 Output Parameter: 103 . its - number of iterations used 104 105 */ 106 #undef __FUNCT__ 107 #define __FUNCT__ "KSPSolve_CG" 108 PetscErrorCode KSPSolve_CG(KSP ksp) 109 { 110 PetscErrorCode ierr; 111 int i,stored_max_it,eigs; 112 PetscScalar dpi,a = 1.0,beta,betaold = 1.0,b,*e = 0,*d = 0,mone = -1.0,ma; 113 PetscReal dp = 0.0; 114 Vec X,B,Z,R,P; 115 KSP_CG *cg; 116 Mat Amat,Pmat; 117 MatStructure pflag; 118 PetscTruth diagonalscale; 119 120 PetscFunctionBegin; 121 ierr = PCDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr); 122 if (diagonalscale) SETERRQ1(1,"Krylov method %s does not support diagonal scaling",ksp->type_name); 123 124 cg = (KSP_CG*)ksp->data; 125 eigs = ksp->calc_sings; 126 stored_max_it = ksp->max_it; 127 X = ksp->vec_sol; 128 B = ksp->vec_rhs; 129 R = ksp->work[0]; 130 Z = ksp->work[1]; 131 P = ksp->work[2]; 132 133 #if !defined(PETSC_USE_COMPLEX) 134 #define VecXDot(x,y,a) VecDot(x,y,a) 135 #else 136 #define VecXDot(x,y,a) (((cg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a)) 137 #endif 138 139 if (eigs) {e = cg->e; d = cg->d; e[0] = 0.0; b = 0.0; } 140 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr); 141 142 ksp->its = 0; 143 if (!ksp->guess_zero) { 144 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /* r <- b - Ax */ 145 ierr = VecAYPX(&mone,B,R);CHKERRQ(ierr); 146 } else { 147 ierr = VecCopy(B,R);CHKERRQ(ierr); /* r <- b (x is 0) */ 148 } 149 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 150 ierr = VecXDot(Z,R,&beta);CHKERRQ(ierr); 151 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 152 ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr); /* dp <- z'*z */ 153 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 154 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- r'*r */ 155 } else if (ksp->normtype == KSP_NATURAL_NORM) { 156 dp = sqrt(PetscAbsScalar(beta)); 157 } else dp = 0.0; 158 KSPLogResidualHistory(ksp,dp); 159 KSPMonitor(ksp,0,dp); /* call any registered monitor routines */ 160 ksp->rnorm = dp; 161 162 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */ 163 if (ksp->reason) PetscFunctionReturn(0); 164 165 i = 0; 166 do { 167 ksp->its = i+1; 168 ierr = VecXDot(Z,R,&beta);CHKERRQ(ierr); /* beta <- r'z */ 169 if (beta == 0.0) { 170 ksp->reason = KSP_CONVERGED_ATOL; 171 PetscLogInfo(ksp,"KSPSolve_CG:converged due to beta = 0"); 172 break; 173 #if !defined(PETSC_USE_COMPLEX) 174 } else if (beta < 0.0) { 175 ksp->reason = KSP_DIVERGED_INDEFINITE_PC; 176 PetscLogInfo(ksp,"KSPSolve_CG:diverging due to indefinite preconditioner"); 177 break; 178 #endif 179 } 180 if (!i) { 181 ierr = VecCopy(Z,P);CHKERRQ(ierr); /* p <- z */ 182 } else { 183 b = beta/betaold; 184 if (eigs) { 185 if (ksp->max_it != stored_max_it) { 186 SETERRQ(1,"Can not change maxit AND calculate eigenvalues"); 187 } 188 e[i] = sqrt(PetscAbsScalar(b))/a; 189 } 190 ierr = VecAYPX(&b,Z,P);CHKERRQ(ierr); /* p <- z + b* p */ 191 } 192 betaold = beta; 193 ierr = KSP_MatMult(ksp,Amat,P,Z);CHKERRQ(ierr); /* z <- Kp */ 194 ierr = VecXDot(P,Z,&dpi);CHKERRQ(ierr); /* dpi <- z'p */ 195 a = beta/dpi; /* a = beta/p'z */ 196 if (eigs) { 197 d[i] = sqrt(PetscAbsScalar(b))*e[i] + 1.0/a; 198 } 199 ierr = VecAXPY(&a,P,X);CHKERRQ(ierr); /* x <- x + ap */ 200 ma = -a; VecAXPY(&ma,Z,R); /* r <- r - az */ 201 if (ksp->normtype == KSP_PRECONDITIONED_NORM) { 202 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 203 ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr); /* dp <- z'*z */ 204 } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) { 205 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- r'*r */ 206 } else if (ksp->normtype == KSP_NATURAL_NORM) { 207 dp = sqrt(PetscAbsScalar(beta)); 208 } else { 209 dp = 0.0; 210 } 211 ksp->rnorm = dp; 212 KSPLogResidualHistory(ksp,dp); 213 KSPMonitor(ksp,i+1,dp); 214 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 215 if (ksp->reason) break; 216 if (ksp->normtype != KSP_PRECONDITIONED_NORM) { 217 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 218 } 219 i++; 220 } while (i<ksp->max_it); 221 if (i == ksp->max_it) { 222 ksp->reason = KSP_DIVERGED_ITS; 223 } 224 PetscFunctionReturn(0); 225 } 226 /* 227 KSPDestroy_CG - Frees all memory space used by the Krylov method 228 229 */ 230 #undef __FUNCT__ 231 #define __FUNCT__ "KSPDestroy_CG" 232 PetscErrorCode KSPDestroy_CG(KSP ksp) 233 { 234 KSP_CG *cg = (KSP_CG*)ksp->data; 235 PetscErrorCode ierr; 236 237 PetscFunctionBegin; 238 /* free space used for singular value calculations */ 239 if (ksp->calc_sings) { 240 ierr = PetscFree(cg->e);CHKERRQ(ierr); 241 ierr = PetscFree(cg->ee);CHKERRQ(ierr); 242 } 243 244 ierr = KSPDefaultFreeWork(ksp);CHKERRQ(ierr); 245 246 /* free the context variable */ 247 ierr = PetscFree(cg);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 /* 252 KSPView_CG - Prints information about the current Krylov method being used 253 254 Currently this only prints information to a file (or stdout) about the 255 symmetry of the problem. If your Krylov method has special options or 256 flags that information should be printed here. 257 258 */ 259 #undef __FUNCT__ 260 #define __FUNCT__ "KSPView_CG" 261 PetscErrorCode KSPView_CG(KSP ksp,PetscViewer viewer) 262 { 263 #if defined(PETSC_USE_COMPLEX) 264 KSP_CG *cg = (KSP_CG *)ksp->data; 265 PetscErrorCode ierr; 266 PetscTruth iascii; 267 268 PetscFunctionBegin; 269 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 270 if (iascii) { 271 if (cg->type == KSP_CG_HERMITIAN) { 272 ierr = PetscViewerASCIIPrintf(viewer," CG: variant for complex, Hermitian system\n");CHKERRQ(ierr); 273 } else if (cg->type == KSP_CG_SYMMETRIC) { 274 ierr = PetscViewerASCIIPrintf(viewer," CG: variant for complex, symmetric system\n");CHKERRQ(ierr); 275 } else { 276 ierr = PetscViewerASCIIPrintf(viewer," CG: unknown variant\n");CHKERRQ(ierr); 277 } 278 } else { 279 SETERRQ1(1,"Viewer type %s not supported for KSP cg",((PetscObject)viewer)->type_name); 280 } 281 #endif 282 PetscFunctionReturn(0); 283 } 284 285 /* 286 KSPSetFromOptions_CG - Checks the options database for options related to the 287 conjugate gradient method. 288 */ 289 #undef __FUNCT__ 290 #define __FUNCT__ "KSPSetFromOptions_CG" 291 PetscErrorCode KSPSetFromOptions_CG(KSP ksp) 292 { 293 #if defined(PETSC_USE_COMPLEX) 294 PetscErrorCode ierr; 295 PetscTruth flg; 296 #endif 297 298 PetscFunctionBegin; 299 #if defined(PETSC_USE_COMPLEX) 300 ierr = PetscOptionsHead("KSP CG options");CHKERRQ(ierr); 301 ierr = PetscOptionsLogicalGroupBegin("-ksp_cg_Hermitian","Matrix is Hermitian","KSPCGSetType",&flg);CHKERRQ(ierr); 302 if (flg) { ierr = KSPCGSetType(ksp,KSP_CG_HERMITIAN);CHKERRQ(ierr); } 303 ierr = PetscOptionsLogicalGroupEnd("-ksp_cg_symmetric","Matrix is complex symmetric, not Hermitian","KSPCGSetType",&flg);CHKERRQ(ierr); 304 if (flg) { ierr = KSPCGSetType(ksp,KSP_CG_SYMMETRIC);CHKERRQ(ierr); } 305 ierr = PetscOptionsTail();CHKERRQ(ierr); 306 #endif 307 PetscFunctionReturn(0); 308 } 309 310 /* 311 KSPCGSetType_CG - This is an option that is SPECIFIC to this particular Krylov method. 312 This routine is registered below in KSPCreate_CG() and called from the 313 routine KSPCGSetType() (see the file cgtype.c). 314 315 This must be wrapped in an EXTERN_C_BEGIN to be dynamically linkable in C++ 316 */ 317 EXTERN_C_BEGIN 318 #undef __FUNCT__ 319 #define __FUNCT__ "KSPCGSetType_CG" 320 PetscErrorCode KSPCGSetType_CG(KSP ksp,KSPCGType type) 321 { 322 KSP_CG *cg; 323 324 PetscFunctionBegin; 325 cg = (KSP_CG *)ksp->data; 326 cg->type = type; 327 PetscFunctionReturn(0); 328 } 329 EXTERN_C_END 330 331 /* 332 KSPCreate_CG - Creates the data structure for the Krylov method CG and sets the 333 function pointers for all the routines it needs to call (KSPSolve_CG() etc) 334 335 It must be wrapped in EXTERN_C_BEGIN to be dynamically linkable in C++ 336 */ 337 /*MC 338 KSPCG - The preconditioned conjugate gradient (PCG) iterative method 339 340 Options Database Keys: 341 + -ksp_cg_Hermitian - (for complex matrices only) indicates the matrix is Hermitian 342 - -ksp_cg_symmetric - (for complex matrices only) indicates the matrix is symmetric 343 344 Level: beginner 345 346 Notes: The PCG method requires both the matrix and preconditioner to 347 be symmetric positive (semi) definite 348 349 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, 350 KSPCGSetType() 351 352 M*/ 353 EXTERN_C_BEGIN 354 #undef __FUNCT__ 355 #define __FUNCT__ "KSPCreate_CG" 356 PetscErrorCode KSPCreate_CG(KSP ksp) 357 { 358 PetscErrorCode ierr; 359 KSP_CG *cg; 360 361 PetscFunctionBegin; 362 ierr = PetscNew(KSP_CG,&cg);CHKERRQ(ierr); 363 ierr = PetscMemzero(cg,sizeof(KSP_CG));CHKERRQ(ierr); 364 PetscLogObjectMemory(ksp,sizeof(KSP_CG)); 365 #if !defined(PETSC_USE_COMPLEX) 366 cg->type = KSP_CG_SYMMETRIC; 367 #else 368 cg->type = KSP_CG_HERMITIAN; 369 #endif 370 ksp->data = (void*)cg; 371 ksp->pc_side = PC_LEFT; 372 373 /* 374 Sets the functions that are associated with this data structure 375 (in C++ this is the same as defining virtual functions) 376 */ 377 ksp->ops->setup = KSPSetUp_CG; 378 ksp->ops->solve = KSPSolve_CG; 379 ksp->ops->destroy = KSPDestroy_CG; 380 ksp->ops->view = KSPView_CG; 381 ksp->ops->setfromoptions = KSPSetFromOptions_CG; 382 ksp->ops->buildsolution = KSPDefaultBuildSolution; 383 ksp->ops->buildresidual = KSPDefaultBuildResidual; 384 385 /* 386 Attach the function KSPCGSetType_CG() to this object. The routine 387 KSPCGSetType() checks for this attached function and calls it if it finds 388 it. (Sort of like a dynamic member function that can be added at run time 389 */ 390 ierr = PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPCGSetType_C","KSPCGSetType_CG", 391 KSPCGSetType_CG);CHKERRQ(ierr); 392 PetscFunctionReturn(0); 393 } 394 EXTERN_C_END 395 396 397 398 399