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