1 /* 2 This file implements the FCG (Flexible Conjugate Gradient) method 3 4 */ 5 6 #include <../src/ksp/ksp/impls/fcg/fcgimpl.h> /*I "petscksp.h" I*/ 7 extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP,PetscReal*,PetscReal*); 8 extern PetscErrorCode KSPComputeEigenvalues_CG(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt*); 9 10 #define KSPFCG_DEFAULT_MMAX 30 /* maximum number of search directions to keep */ 11 #define KSPFCG_DEFAULT_NPREALLOC 10 /* number of search directions to preallocate */ 12 #define KSPFCG_DEFAULT_VECB 5 /* number of search directions to allocate each time new direction vectors are needed */ 13 #define KSPFCG_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY 14 15 static PetscErrorCode KSPAllocateVectors_FCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize) 16 { 17 PetscErrorCode ierr; 18 PetscInt i; 19 KSP_FCG *fcg = (KSP_FCG*)ksp->data; 20 PetscInt nnewvecs, nvecsprev; 21 22 PetscFunctionBegin; 23 /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */ 24 if (fcg->nvecs < PetscMin(fcg->mmax+1,nvecsneeded)){ 25 nvecsprev = fcg->nvecs; 26 nnewvecs = PetscMin(PetscMax(nvecsneeded-fcg->nvecs,chunksize),fcg->mmax+1-fcg->nvecs); 27 ierr = KSPCreateVecs(ksp,nnewvecs,&fcg->pCvecs[fcg->nchunks],0,NULL);CHKERRQ(ierr); 28 ierr = PetscLogObjectParents((PetscObject)ksp,nnewvecs,fcg->pCvecs[fcg->nchunks]);CHKERRQ(ierr); 29 ierr = KSPCreateVecs(ksp,nnewvecs,&fcg->pPvecs[fcg->nchunks],0,NULL);CHKERRQ(ierr); 30 ierr = PetscLogObjectParents((PetscObject)ksp,nnewvecs,fcg->pPvecs[fcg->nchunks]);CHKERRQ(ierr); 31 fcg->nvecs += nnewvecs; 32 for (i=0;i<nnewvecs;++i){ 33 fcg->Cvecs[nvecsprev + i] = fcg->pCvecs[fcg->nchunks][i]; 34 fcg->Pvecs[nvecsprev + i] = fcg->pPvecs[fcg->nchunks][i]; 35 } 36 fcg->chunksizes[fcg->nchunks] = nnewvecs; 37 ++fcg->nchunks; 38 } 39 PetscFunctionReturn(0); 40 } 41 42 static PetscErrorCode KSPSetUp_FCG(KSP ksp) 43 { 44 PetscErrorCode ierr; 45 KSP_FCG *fcg = (KSP_FCG*)ksp->data; 46 PetscInt maxit = ksp->max_it; 47 const PetscInt nworkstd = 2; 48 49 PetscFunctionBegin; 50 51 /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */ 52 ierr = KSPSetWorkVecs(ksp,nworkstd);CHKERRQ(ierr); 53 54 /* Allocated space for pointers to additional work vectors 55 note that mmax is the number of previous directions, so we add 1 for the current direction, 56 and an extra 1 for the prealloc (which might be empty) */ 57 ierr = PetscMalloc5(fcg->mmax+1,&fcg->Pvecs,fcg->mmax+1,&fcg->Cvecs,fcg->mmax+1,&fcg->pPvecs,fcg->mmax+1,&fcg->pCvecs,fcg->mmax+2,&fcg->chunksizes);CHKERRQ(ierr); 58 ierr = PetscLogObjectMemory((PetscObject)ksp,2*(fcg->mmax+1)*sizeof(Vec*) + 2*(fcg->mmax + 1)*sizeof(Vec**) + (fcg->mmax + 2)*sizeof(PetscInt));CHKERRQ(ierr); 59 60 /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */ 61 if(fcg->nprealloc > fcg->mmax+1){ 62 ierr = PetscInfo2(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",fcg->nprealloc, fcg->mmax+1);CHKERRQ(ierr); 63 } 64 65 /* Preallocate additional work vectors */ 66 ierr = KSPAllocateVectors_FCG(ksp,fcg->nprealloc,fcg->nprealloc);CHKERRQ(ierr); 67 /* 68 If user requested computations of eigenvalues then allocate work 69 work space needed 70 */ 71 if (ksp->calc_sings) { 72 /* get space to store tridiagonal matrix for Lanczos */ 73 ierr = PetscMalloc4(maxit,&fcg->e,maxit,&fcg->d,maxit,&fcg->ee,maxit,&fcg->dd);CHKERRQ(ierr); 74 ierr = PetscLogObjectMemory((PetscObject)ksp,2*(maxit+1)*(sizeof(PetscScalar)+sizeof(PetscReal)));CHKERRQ(ierr); 75 76 ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG; 77 ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG; 78 } 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode KSPSolve_FCG(KSP ksp) 83 { 84 PetscErrorCode ierr; 85 PetscInt i,k,idx,mi; 86 KSP_FCG *fcg = (KSP_FCG*)ksp->data; 87 PetscScalar alpha=0.0,beta = 0.0,dpi,s; 88 PetscReal dp=0.0; 89 Vec B,R,Z,X,Pcurr,Ccurr; 90 Mat Amat,Pmat; 91 PetscInt eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/ 92 PetscInt stored_max_it = ksp->max_it; 93 PetscScalar alphaold = 0,betaold = 1.0,*e = 0,*d = 0;/* Variables for eigen estimation - FINISH */ 94 95 PetscFunctionBegin; 96 97 #define VecXDot(x,y,a) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a)) 98 #define VecXMDot(a,b,c,d) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(a,b,c,d) : VecMTDot(a,b,c,d)) 99 100 X = ksp->vec_sol; 101 B = ksp->vec_rhs; 102 R = ksp->work[0]; 103 Z = ksp->work[1]; 104 105 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr); 106 if (eigs) {e = fcg->e; d = fcg->d; e[0] = 0.0; } 107 /* Compute initial residual needed for convergence check*/ 108 ksp->its = 0; 109 if (!ksp->guess_zero) { 110 ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); 111 ierr = VecAYPX(R,-1.0,B);CHKERRQ(ierr); /* r <- b - Ax */ 112 } else { 113 ierr = VecCopy(B,R);CHKERRQ(ierr); /* r <- b (x is 0) */ 114 } 115 switch (ksp->normtype) { 116 case KSP_NORM_PRECONDITIONED: 117 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 118 ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr); /* dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */ 119 break; 120 case KSP_NORM_UNPRECONDITIONED: 121 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */ 122 break; 123 case KSP_NORM_NATURAL: 124 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 125 ierr = VecXDot(R,Z,&s);CHKERRQ(ierr); 126 dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */ 127 break; 128 case KSP_NORM_NONE: 129 dp = 0.0; 130 break; 131 default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]); 132 } 133 134 /* Initial Convergence Check */ 135 ierr = KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 136 ierr = KSPMonitor(ksp,0,dp);CHKERRQ(ierr); 137 ksp->rnorm = dp; 138 if (ksp->normtype == KSP_NORM_NONE) { 139 ierr = KSPConvergedSkip(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 140 } else { 141 ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 142 } 143 if (ksp->reason) PetscFunctionReturn(0); 144 145 /* Apply PC if not already done for convergence check */ 146 if(ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE){ 147 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 148 } 149 150 i = 0; 151 do { 152 ksp->its = i+1; 153 154 /* If needbe, allocate a new chunk of vectors in P and C */ 155 ierr = KSPAllocateVectors_FCG(ksp,i+1,fcg->vecb);CHKERRQ(ierr); 156 157 /* Note that we wrap around and start clobbering old vectors */ 158 idx = i % (fcg->mmax+1); 159 Pcurr = fcg->Pvecs[idx]; 160 Ccurr = fcg->Cvecs[idx]; 161 162 /* number of old directions to orthogonalize against */ 163 switch(fcg->truncstrat){ 164 case KSP_FCD_TRUNC_TYPE_STANDARD: 165 mi = fcg->mmax; 166 break; 167 case KSP_FCD_TRUNC_TYPE_NOTAY: 168 mi = ((i-1) % fcg->mmax)+1; 169 break; 170 default: 171 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized Truncation Strategy"); 172 } 173 174 /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/ 175 ierr = VecCopy(Z,Pcurr);CHKERRQ(ierr); 176 177 { 178 PetscInt l,ndots; 179 180 l = PetscMax(0,i-mi); 181 ndots = i-l; 182 if (ndots){ 183 PetscInt j; 184 Vec *Pold, *Cold; 185 PetscScalar *dots; 186 187 ierr = PetscMalloc3(ndots,&dots,ndots,&Cold,ndots,&Pold);CHKERRQ(ierr); 188 for(k=l,j=0;j<ndots;++k,++j){ 189 idx = k % (fcg->mmax+1); 190 Cold[j] = fcg->Cvecs[idx]; 191 Pold[j] = fcg->Pvecs[idx]; 192 } 193 ierr = VecXMDot(Z,ndots,Cold,dots);CHKERRQ(ierr); 194 for(k=0;k<ndots;++k){ 195 dots[k] = -dots[k]; 196 } 197 ierr = VecMAXPY(Pcurr,ndots,dots,Pold);CHKERRQ(ierr); 198 ierr = PetscFree3(dots,Cold,Pold);CHKERRQ(ierr); 199 } 200 } 201 202 /* Update X and R */ 203 betaold = beta; 204 ierr = VecXDot(Pcurr,R,&beta);CHKERRQ(ierr); /* beta <- pi'*r */ 205 ierr = KSP_MatMult(ksp,Amat,Pcurr,Ccurr);CHKERRQ(ierr); /* w <- A*pi (stored in ci) */ 206 ierr = VecXDot(Pcurr,Ccurr,&dpi);CHKERRQ(ierr); /* dpi <- pi'*w */ 207 alphaold = alpha; 208 alpha = beta / dpi; /* alpha <- beta/dpi */ 209 ierr = VecAXPY(X,alpha,Pcurr);CHKERRQ(ierr); /* x <- x + alpha * pi */ 210 ierr = VecAXPY(R,-alpha,Ccurr);CHKERRQ(ierr); /* r <- r - alpha * wi */ 211 212 /* Compute norm for convergence check */ 213 switch (ksp->normtype) { 214 case KSP_NORM_PRECONDITIONED: 215 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 216 ierr = VecNorm(Z,NORM_2,&dp);CHKERRQ(ierr); /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */ 217 break; 218 case KSP_NORM_UNPRECONDITIONED: 219 ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */ 220 break; 221 case KSP_NORM_NATURAL: 222 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 223 ierr = VecXDot(R,Z,&s);CHKERRQ(ierr); 224 dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */ 225 break; 226 case KSP_NORM_NONE: 227 dp = 0.0; 228 break; 229 default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]); 230 } 231 232 /* Check for convergence */ 233 ksp->rnorm = dp; 234 KSPLogResidualHistory(ksp,dp);CHKERRQ(ierr); 235 ierr = KSPMonitor(ksp,i+1,dp);CHKERRQ(ierr); 236 ierr = (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); 237 if (ksp->reason) break; 238 239 /* Apply PC if not already done for convergence check */ 240 if(ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE){ 241 ierr = KSP_PCApply(ksp,R,Z);CHKERRQ(ierr); /* z <- Br */ 242 } 243 244 /* Compute current C (which is W/dpi) */ 245 ierr = VecScale(Ccurr,1.0/dpi);CHKERRQ(ierr); /* w <- ci/dpi */ 246 247 if (eigs) { 248 if (i > 0) { 249 if (ksp->max_it != stored_max_it) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Can not change maxit AND calculate eigenvalues"); 250 e[i] = PetscSqrtReal(PetscAbsScalar(beta/betaold))/alphaold; 251 d[i] = PetscSqrtReal(PetscAbsScalar(beta/betaold))*e[i] + 1.0/alpha; 252 } else { 253 d[i] = PetscSqrtReal(PetscAbsScalar(beta))*e[i] + 1.0/alpha; 254 } 255 fcg->ned = ksp->its-1; 256 } 257 ++i; 258 } while (i<ksp->max_it); 259 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS; 260 PetscFunctionReturn(0); 261 } 262 263 static PetscErrorCode KSPDestroy_FCG(KSP ksp) 264 { 265 PetscErrorCode ierr; 266 PetscInt i; 267 KSP_FCG *fcg = (KSP_FCG*)ksp->data; 268 269 PetscFunctionBegin; 270 271 /* Destroy "standard" work vecs */ 272 VecDestroyVecs(ksp->nwork,&ksp->work); 273 274 /* Destroy P and C vectors and the arrays that manage pointers to them */ 275 if (fcg->nvecs){ 276 for (i=0;i<fcg->nchunks;++i){ 277 ierr = VecDestroyVecs(fcg->chunksizes[i],&fcg->pPvecs[i]);CHKERRQ(ierr); 278 ierr = VecDestroyVecs(fcg->chunksizes[i],&fcg->pCvecs[i]);CHKERRQ(ierr); 279 } 280 } 281 ierr = PetscFree5(fcg->Pvecs,fcg->Cvecs,fcg->pPvecs,fcg->pCvecs,fcg->chunksizes);CHKERRQ(ierr); 282 /* free space used for singular value calculations */ 283 if (ksp->calc_sings) { 284 ierr = PetscFree4(fcg->e,fcg->d,fcg->ee,fcg->dd);CHKERRQ(ierr); 285 } 286 ierr = KSPDestroyDefault(ksp);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 static PetscErrorCode KSPView_FCG(KSP ksp,PetscViewer viewer) 291 { 292 KSP_FCG *fcg = (KSP_FCG*)ksp->data; 293 PetscErrorCode ierr; 294 PetscBool iascii,isstring; 295 const char *truncstr; 296 297 PetscFunctionBegin; 298 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 299 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 300 301 if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) truncstr = "Using standard truncation strategy"; 302 else if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) truncstr = "Using Notay's truncation strategy"; 303 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Undefined FCG truncation strategy"); 304 305 if (iascii) { 306 ierr = PetscViewerASCIIPrintf(viewer," m_max=%D\n",fcg->mmax);CHKERRQ(ierr); 307 ierr = PetscViewerASCIIPrintf(viewer," preallocated %D directions\n",PetscMin(fcg->nprealloc,fcg->mmax+1));CHKERRQ(ierr); 308 ierr = PetscViewerASCIIPrintf(viewer," %s\n",truncstr);CHKERRQ(ierr); 309 } else if (isstring) { 310 ierr = PetscViewerStringSPrintf(viewer,"m_max %D nprealloc %D %s",fcg->mmax,fcg->nprealloc,truncstr);CHKERRQ(ierr); 311 } 312 PetscFunctionReturn(0); 313 } 314 315 /*@ 316 KSPFCGSetMmax - set the maximum number of previous directions FCG will store for orthogonalization 317 318 Note: mmax + 1 directions are stored (mmax previous ones along with a current one) 319 and whether all are used in each iteration also depends on the truncation strategy 320 (see KSPFCGSetTruncationType()) 321 322 Logically Collective on KSP 323 324 Input Parameters: 325 + ksp - the Krylov space context 326 - mmax - the maximum number of previous directions to orthogonalize againt 327 328 Level: intermediate 329 330 Options Database: 331 . -ksp_fcg_mmax <N> 332 333 .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc() 334 @*/ 335 PetscErrorCode KSPFCGSetMmax(KSP ksp,PetscInt mmax) 336 { 337 KSP_FCG *fcg = (KSP_FCG*)ksp->data; 338 339 PetscFunctionBegin; 340 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 341 PetscValidLogicalCollectiveEnum(ksp,mmax,2); 342 fcg->mmax = mmax; 343 PetscFunctionReturn(0); 344 } 345 346 /*@ 347 KSPFCGGetMmax - get the maximum number of previous directions FCG will store 348 349 Note: FCG stores mmax+1 directions at most (mmax previous ones, and one current one) 350 351 Not Collective 352 353 Input Parameter: 354 . ksp - the Krylov space context 355 356 Output Parameter: 357 . mmax - the maximum number of previous directons allowed for orthogonalization 358 359 Options Database: 360 . -ksp_fcg_mmax <N> 361 362 Level: intermediate 363 364 .keywords: KSP, FCG, truncation 365 366 .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc(), KSPFCGSetMmax() 367 @*/ 368 369 PetscErrorCode KSPFCGGetMmax(KSP ksp,PetscInt *mmax) 370 { 371 KSP_FCG *fcg=(KSP_FCG*)ksp->data; 372 373 PetscFunctionBegin; 374 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 375 *mmax = fcg->mmax; 376 PetscFunctionReturn(0); 377 } 378 379 /*@ 380 KSPFCGSetNprealloc - set the number of directions to preallocate with FCG 381 382 Logically Collective on KSP 383 384 Input Parameters: 385 + ksp - the Krylov space context 386 - nprealloc - the number of vectors to preallocate 387 388 Level: advanced 389 390 Options Database: 391 . -ksp_fcg_nprealloc <N> - number of directions to preallocate 392 393 .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc() 394 @*/ 395 PetscErrorCode KSPFCGSetNprealloc(KSP ksp,PetscInt nprealloc) 396 { 397 KSP_FCG *fcg=(KSP_FCG*)ksp->data; 398 399 PetscFunctionBegin; 400 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 401 PetscValidLogicalCollectiveEnum(ksp,nprealloc,2); 402 if(nprealloc > fcg->mmax+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Cannot preallocate more than m_max+1 vectors"); 403 fcg->nprealloc = nprealloc; 404 PetscFunctionReturn(0); 405 } 406 407 /*@ 408 KSPFCGGetNprealloc - get the number of directions preallocate by FCG 409 410 Not Collective 411 412 Input Parameter: 413 . ksp - the Krylov space context 414 415 Output Parameter: 416 . nprealloc - the number of directions preallocated 417 418 Level: advanced 419 420 .keywords: KSP, FCG, truncation 421 422 .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGSetNprealloc() 423 @*/ 424 PetscErrorCode KSPFCGGetNprealloc(KSP ksp,PetscInt *nprealloc) 425 { 426 KSP_FCG *fcg=(KSP_FCG*)ksp->data; 427 428 PetscFunctionBegin; 429 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 430 *nprealloc = fcg->nprealloc; 431 PetscFunctionReturn(0); 432 } 433 434 /*@ 435 KSPFCGSetTruncationType - specify how many of its stored previous directions FCG uses during orthoganalization 436 437 Logically Collective on KSP 438 439 KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions 440 KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,.. 441 442 Input Parameters: 443 + ksp - the Krylov space context 444 - truncstrat - the choice of strategy 445 446 Level: intermediate 447 448 Options Database: 449 . -ksp_fcg_truncation_type <standard, notay> - specify how many of its stored previous directions FCG uses during orthoganalization 450 451 .seealso: KSPFCDTruncationType, KSPFCGGetTruncationType 452 @*/ 453 PetscErrorCode KSPFCGSetTruncationType(KSP ksp,KSPFCDTruncationType truncstrat) 454 { 455 KSP_FCG *fcg=(KSP_FCG*)ksp->data; 456 457 PetscFunctionBegin; 458 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 459 PetscValidLogicalCollectiveEnum(ksp,truncstrat,2); 460 fcg->truncstrat=truncstrat; 461 PetscFunctionReturn(0); 462 } 463 464 /*@ 465 KSPFCGGetTruncationType - get the truncation strategy employed by FCG 466 467 Not Collective 468 469 Input Parameter: 470 . ksp - the Krylov space context 471 472 Output Parameter: 473 . truncstrat - the strategy type 474 475 Level: intermediate 476 477 .keywords: KSP, FCG, truncation 478 479 .seealso: KSPFCG, KSPFCGSetTruncationType, KSPFCDTruncationType 480 @*/ 481 PetscErrorCode KSPFCGGetTruncationType(KSP ksp,KSPFCDTruncationType *truncstrat) 482 { 483 KSP_FCG *fcg=(KSP_FCG*)ksp->data; 484 485 PetscFunctionBegin; 486 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 487 *truncstrat=fcg->truncstrat; 488 PetscFunctionReturn(0); 489 } 490 491 static PetscErrorCode KSPSetFromOptions_FCG(PetscOptionItems *PetscOptionsObject,KSP ksp) 492 { 493 PetscErrorCode ierr; 494 KSP_FCG *fcg=(KSP_FCG*)ksp->data; 495 PetscInt mmax,nprealloc; 496 PetscBool flg; 497 498 PetscFunctionBegin; 499 ierr = PetscOptionsHead(PetscOptionsObject,"KSP FCG Options");CHKERRQ(ierr); 500 ierr = PetscOptionsInt("-ksp_fcg_mmax","Maximum number of search directions to store","KSPFCGSetMmax",fcg->mmax,&mmax,&flg);CHKERRQ(ierr); 501 if (flg) { 502 ierr = KSPFCGSetMmax(ksp,mmax);CHKERRQ(ierr); 503 } 504 ierr = PetscOptionsInt("-ksp_fcg_nprealloc","Number of directions to preallocate","KSPFCGSetNprealloc",fcg->nprealloc,&nprealloc,&flg);CHKERRQ(ierr); 505 if (flg) { 506 ierr = KSPFCGSetNprealloc(ksp,nprealloc);CHKERRQ(ierr); 507 } 508 ierr = PetscOptionsEnum("-ksp_fcg_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCDTruncationTypes,(PetscEnum)fcg->truncstrat,(PetscEnum*)&fcg->truncstrat,NULL);CHKERRQ(ierr); 509 ierr = PetscOptionsTail();CHKERRQ(ierr); 510 PetscFunctionReturn(0); 511 } 512 513 /*MC 514 KSPFCG - Implements the Flexible Conjugate Gradient method (FCG) 515 516 Options Database Keys: 517 + -ksp_fcg_mmax <N> - maximum number of search directions 518 . -ksp_fcg_nprealloc <N> - number of directions to preallocate 519 - -ksp_fcg_truncation_type <standard,notay> - truncation approach for directions 520 521 Contributed by Patrick Sanan 522 523 Notes: 524 Supports left preconditioning only. 525 526 Level: beginner 527 528 References: 529 + 1. - Notay, Y."Flexible Conjugate Gradients", SIAM J. Sci. Comput. 22:4, 2000 530 - 2. - Axelsson, O. and Vassilevski, P. S. "A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable step Preconditioning", 531 SIAM J. Matrix Anal. Appl. 12:4, 1991 532 533 .seealso : KSPGCR, KSPFGMRES, KSPCG, KSPFCGSetMmax(), KSPFCGGetMmax(), KSPFCGSetNprealloc(), KSPFCGGetNprealloc(), KSPFCGSetTruncationType(), KSPFCGGetTruncationType() 534 535 M*/ 536 PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp) 537 { 538 PetscErrorCode ierr; 539 KSP_FCG *fcg; 540 541 PetscFunctionBegin; 542 ierr = PetscNewLog(ksp,&fcg);CHKERRQ(ierr); 543 #if !defined(PETSC_USE_COMPLEX) 544 fcg->type = KSP_CG_SYMMETRIC; 545 #else 546 fcg->type = KSP_CG_HERMITIAN; 547 #endif 548 fcg->mmax = KSPFCG_DEFAULT_MMAX; 549 fcg->nprealloc = KSPFCG_DEFAULT_NPREALLOC; 550 fcg->nvecs = 0; 551 fcg->vecb = KSPFCG_DEFAULT_VECB; 552 fcg->nchunks = 0; 553 fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT; 554 555 ksp->data = (void*)fcg; 556 557 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr); 558 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);CHKERRQ(ierr); 559 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);CHKERRQ(ierr); 560 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr); 561 562 ksp->ops->setup = KSPSetUp_FCG; 563 ksp->ops->solve = KSPSolve_FCG; 564 ksp->ops->destroy = KSPDestroy_FCG; 565 ksp->ops->view = KSPView_FCG; 566 ksp->ops->setfromoptions = KSPSetFromOptions_FCG; 567 ksp->ops->buildsolution = KSPBuildSolutionDefault; 568 ksp->ops->buildresidual = KSPBuildResidualDefault; 569 PetscFunctionReturn(0); 570 } 571 572