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