Actual source code: gd.c

slepc-3.20.2 2024-03-15
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc eigensolver: "gd"

 13:    Method: Generalized Davidson

 15:    Algorithm:

 17:        Generalized Davidson with various subspace extraction and
 18:        restart techniques.

 20:    References:

 22:        [1] E.R. Davidson, "Super-matrix methods", Comput. Phys. Commun.
 23:            53(2):49-60, 1989.

 25:        [2] E. Romero and J.E. Roman, "A parallel implementation of
 26:            Davidson methods for large-scale eigenvalue problems in
 27:            SLEPc", ACM Trans. Math. Software 40(2), Article 13, 2014.
 28: */

 30: #include <slepc/private/epsimpl.h>
 31: #include <../src/eps/impls/davidson/davidson.h>

 33: static PetscErrorCode EPSSetFromOptions_GD(EPS eps,PetscOptionItems *PetscOptionsObject)
 34: {
 35:   PetscBool      flg,flg2,op,orth;
 36:   PetscInt       opi,opi0;

 38:   PetscFunctionBegin;
 39:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Generalized Davidson (GD) Options");

 41:     PetscCall(EPSGDGetKrylovStart(eps,&op));
 42:     PetscCall(PetscOptionsBool("-eps_gd_krylov_start","Start the search subspace with a Krylov basis","EPSGDSetKrylovStart",op,&op,&flg));
 43:     if (flg) PetscCall(EPSGDSetKrylovStart(eps,op));

 45:     PetscCall(EPSGDGetBOrth(eps,&orth));
 46:     PetscCall(PetscOptionsBool("-eps_gd_borth","Use B-orthogonalization in the search subspace","EPSGDSetBOrth",op,&op,&flg));
 47:     if (flg) PetscCall(EPSGDSetBOrth(eps,op));

 49:     PetscCall(EPSGDGetBlockSize(eps,&opi));
 50:     PetscCall(PetscOptionsInt("-eps_gd_blocksize","Number of vectors to add to the search subspace","EPSGDSetBlockSize",opi,&opi,&flg));
 51:     if (flg) PetscCall(EPSGDSetBlockSize(eps,opi));

 53:     PetscCall(EPSGDGetRestart(eps,&opi,&opi0));
 54:     PetscCall(PetscOptionsInt("-eps_gd_minv","Size of the search subspace after restarting","EPSGDSetRestart",opi,&opi,&flg));
 55:     PetscCall(PetscOptionsInt("-eps_gd_plusk","Number of eigenvectors saved from the previous iteration when restarting","EPSGDSetRestart",opi0,&opi0,&flg2));
 56:     if (flg || flg2) PetscCall(EPSGDSetRestart(eps,opi,opi0));

 58:     PetscCall(EPSGDGetInitialSize(eps,&opi));
 59:     PetscCall(PetscOptionsInt("-eps_gd_initial_size","Initial size of the search subspace","EPSGDSetInitialSize",opi,&opi,&flg));
 60:     if (flg) PetscCall(EPSGDSetInitialSize(eps,opi));

 62:     PetscCall(PetscOptionsBool("-eps_gd_double_expansion","Use the doble-expansion variant of GD","EPSGDSetDoubleExpansion",PETSC_FALSE,&op,&flg));
 63:     if (flg) PetscCall(EPSGDSetDoubleExpansion(eps,op));

 65:   PetscOptionsHeadEnd();
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: static PetscErrorCode EPSSetUp_GD(EPS eps)
 70: {
 71:   PetscBool      t;
 72:   KSP            ksp;

 74:   PetscFunctionBegin;
 75:   /* Setup common for all davidson solvers */
 76:   PetscCall(EPSSetUp_XD(eps));

 78:   /* Check some constraints */
 79:   PetscCall(STGetKSP(eps->st,&ksp));
 80:   PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&t));
 81:   PetscCheck(t,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"EPSGD only works with KSPPREONLY");
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: static PetscErrorCode EPSView_GD(EPS eps,PetscViewer viewer)
 86: {
 87:   PetscBool      isascii,opb;
 88:   PetscInt       opi,opi0;
 89:   PetscBool      borth;
 90:   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;

 92:   PetscFunctionBegin;
 93:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
 94:   if (isascii) {
 95:     if (data->doubleexp) PetscCall(PetscViewerASCIIPrintf(viewer,"  using double expansion variant (GD2)\n"));
 96:     PetscCall(EPSXDGetBOrth_XD(eps,&borth));
 97:     if (borth) PetscCall(PetscViewerASCIIPrintf(viewer,"  search subspace is B-orthogonalized\n"));
 98:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  search subspace is orthogonalized\n"));
 99:     PetscCall(EPSXDGetBlockSize_XD(eps,&opi));
100:     PetscCall(PetscViewerASCIIPrintf(viewer,"  block size=%" PetscInt_FMT "\n",opi));
101:     PetscCall(EPSXDGetKrylovStart_XD(eps,&opb));
102:     if (!opb) PetscCall(PetscViewerASCIIPrintf(viewer,"  type of the initial subspace: non-Krylov\n"));
103:     else PetscCall(PetscViewerASCIIPrintf(viewer,"  type of the initial subspace: Krylov\n"));
104:     PetscCall(EPSXDGetRestart_XD(eps,&opi,&opi0));
105:     PetscCall(PetscViewerASCIIPrintf(viewer,"  size of the subspace after restarting: %" PetscInt_FMT "\n",opi));
106:     PetscCall(PetscViewerASCIIPrintf(viewer,"  number of vectors after restarting from the previous iteration: %" PetscInt_FMT "\n",opi0));
107:   }
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: static PetscErrorCode EPSDestroy_GD(EPS eps)
112: {
113:   PetscFunctionBegin;
114:   PetscCall(PetscFree(eps->data));
115:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",NULL));
116:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",NULL));
117:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",NULL));
118:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",NULL));
119:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",NULL));
120:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",NULL));
121:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",NULL));
122:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",NULL));
123:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",NULL));
124:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",NULL));
125:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",NULL));
126:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",NULL));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*@
131:    EPSGDSetKrylovStart - Activates or deactivates starting the searching
132:    subspace with a Krylov basis.

134:    Logically Collective

136:    Input Parameters:
137: +  eps - the eigenproblem solver context
138: -  krylovstart - boolean flag

140:    Options Database Key:
141: .  -eps_gd_krylov_start - Activates starting the searching subspace with a
142:     Krylov basis

144:    Level: advanced

146: .seealso: EPSGDGetKrylovStart()
147: @*/
148: PetscErrorCode EPSGDSetKrylovStart(EPS eps,PetscBool krylovstart)
149: {
150:   PetscFunctionBegin;
153:   PetscTryMethod(eps,"EPSGDSetKrylovStart_C",(EPS,PetscBool),(eps,krylovstart));
154:   PetscFunctionReturn(PETSC_SUCCESS);
155: }

157: /*@
158:    EPSGDGetKrylovStart - Returns a flag indicating if the search subspace is started with a
159:    Krylov basis.

161:    Not Collective

163:    Input Parameter:
164: .  eps - the eigenproblem solver context

166:    Output Parameters:
167: .  krylovstart - boolean flag indicating if the search subspace is started
168:    with a Krylov basis

170:    Level: advanced

172: .seealso: EPSGDSetKrylovStart()
173: @*/
174: PetscErrorCode EPSGDGetKrylovStart(EPS eps,PetscBool *krylovstart)
175: {
176:   PetscFunctionBegin;
178:   PetscAssertPointer(krylovstart,2);
179:   PetscUseMethod(eps,"EPSGDGetKrylovStart_C",(EPS,PetscBool*),(eps,krylovstart));
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: /*@
184:    EPSGDSetBlockSize - Sets the number of vectors to be added to the searching space
185:    in every iteration.

187:    Logically Collective

189:    Input Parameters:
190: +  eps - the eigenproblem solver context
191: -  blocksize - number of vectors added to the search space in every iteration

193:    Options Database Key:
194: .  -eps_gd_blocksize - number of vectors added to the search space in every iteration

196:    Level: advanced

198: .seealso: EPSGDSetKrylovStart()
199: @*/
200: PetscErrorCode EPSGDSetBlockSize(EPS eps,PetscInt blocksize)
201: {
202:   PetscFunctionBegin;
205:   PetscTryMethod(eps,"EPSGDSetBlockSize_C",(EPS,PetscInt),(eps,blocksize));
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

209: /*@
210:    EPSGDGetBlockSize - Returns the number of vectors to be added to the searching space
211:    in every iteration.

213:    Not Collective

215:    Input Parameter:
216: .  eps - the eigenproblem solver context

218:    Output Parameter:
219: .  blocksize - number of vectors added to the search space in every iteration

221:    Level: advanced

223: .seealso: EPSGDSetBlockSize()
224: @*/
225: PetscErrorCode EPSGDGetBlockSize(EPS eps,PetscInt *blocksize)
226: {
227:   PetscFunctionBegin;
229:   PetscAssertPointer(blocksize,2);
230:   PetscUseMethod(eps,"EPSGDGetBlockSize_C",(EPS,PetscInt*),(eps,blocksize));
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

234: /*@
235:    EPSGDSetRestart - Sets the number of vectors of the searching space after
236:    restarting and the number of vectors saved from the previous iteration.

238:    Logically Collective

240:    Input Parameters:
241: +  eps - the eigenproblem solver context
242: .  minv - number of vectors of the searching subspace after restarting
243: -  plusk - number of vectors saved from the previous iteration

245:    Options Database Keys:
246: +  -eps_gd_minv - number of vectors of the searching subspace after restarting
247: -  -eps_gd_plusk - number of vectors saved from the previous iteration

249:    Level: advanced

251: .seealso: EPSGDGetRestart()
252: @*/
253: PetscErrorCode EPSGDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)
254: {
255:   PetscFunctionBegin;
259:   PetscTryMethod(eps,"EPSGDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: /*@
264:    EPSGDGetRestart - Gets the number of vectors of the searching space after
265:    restarting and the number of vectors saved from the previous iteration.

267:    Not Collective

269:    Input Parameter:
270: .  eps - the eigenproblem solver context

272:    Output Parameters:
273: +  minv - number of vectors of the searching subspace after restarting
274: -  plusk - number of vectors saved from the previous iteration

276:    Level: advanced

278: .seealso: EPSGDSetRestart()
279: @*/
280: PetscErrorCode EPSGDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)
281: {
282:   PetscFunctionBegin;
284:   PetscUseMethod(eps,"EPSGDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

288: /*@
289:    EPSGDSetInitialSize - Sets the initial size of the searching space.

291:    Logically Collective

293:    Input Parameters:
294: +  eps - the eigenproblem solver context
295: -  initialsize - number of vectors of the initial searching subspace

297:    Options Database Key:
298: .  -eps_gd_initial_size - number of vectors of the initial searching subspace

300:    Notes:
301:    If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
302:    EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
303:    provided vectors are not enough, the solver completes the subspace with
304:    random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
305:    gets the first vector provided by the user or, if not available, a random vector,
306:    and expands the Krylov basis up to initialsize vectors.

308:    Level: advanced

310: .seealso: EPSGDGetInitialSize(), EPSGDGetKrylovStart()
311: @*/
312: PetscErrorCode EPSGDSetInitialSize(EPS eps,PetscInt initialsize)
313: {
314:   PetscFunctionBegin;
317:   PetscTryMethod(eps,"EPSGDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
318:   PetscFunctionReturn(PETSC_SUCCESS);
319: }

321: /*@
322:    EPSGDGetInitialSize - Returns the initial size of the searching space.

324:    Not Collective

326:    Input Parameter:
327: .  eps - the eigenproblem solver context

329:    Output Parameter:
330: .  initialsize - number of vectors of the initial searching subspace

332:    Notes:
333:    If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
334:    EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
335:    provided vectors are not enough, the solver completes the subspace with
336:    random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
337:    gets the first vector provided by the user or, if not available, a random vector,
338:    and expands the Krylov basis up to initialsize vectors.

340:    Level: advanced

342: .seealso: EPSGDSetInitialSize(), EPSGDGetKrylovStart()
343: @*/
344: PetscErrorCode EPSGDGetInitialSize(EPS eps,PetscInt *initialsize)
345: {
346:   PetscFunctionBegin;
348:   PetscAssertPointer(initialsize,2);
349:   PetscUseMethod(eps,"EPSGDGetInitialSize_C",(EPS,PetscInt*),(eps,initialsize));
350:   PetscFunctionReturn(PETSC_SUCCESS);
351: }

353: /*@
354:    EPSGDSetBOrth - Selects the orthogonalization that will be used in the search
355:    subspace in case of generalized Hermitian problems.

357:    Logically Collective

359:    Input Parameters:
360: +  eps   - the eigenproblem solver context
361: -  borth - whether to B-orthogonalize the search subspace

363:    Options Database Key:
364: .  -eps_gd_borth - Set the orthogonalization used in the search subspace

366:    Level: advanced

368: .seealso: EPSGDGetBOrth()
369: @*/
370: PetscErrorCode EPSGDSetBOrth(EPS eps,PetscBool borth)
371: {
372:   PetscFunctionBegin;
375:   PetscTryMethod(eps,"EPSGDSetBOrth_C",(EPS,PetscBool),(eps,borth));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: /*@
380:    EPSGDGetBOrth - Returns the orthogonalization used in the search
381:    subspace in case of generalized Hermitian problems.

383:    Not Collective

385:    Input Parameter:
386: .  eps - the eigenproblem solver context

388:    Output Parameters:
389: .  borth - whether to B-orthogonalize the search subspace

391:    Level: advanced

393: .seealso: EPSGDSetBOrth()
394: @*/
395: PetscErrorCode EPSGDGetBOrth(EPS eps,PetscBool *borth)
396: {
397:   PetscFunctionBegin;
399:   PetscAssertPointer(borth,2);
400:   PetscUseMethod(eps,"EPSGDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: static PetscErrorCode EPSGDSetDoubleExpansion_GD(EPS eps,PetscBool doubleexp)
405: {
406:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

408:   PetscFunctionBegin;
409:   data->doubleexp = doubleexp;
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: /*@
414:    EPSGDSetDoubleExpansion - Activate the double expansion variant of GD.

416:    Logically Collective

418:    Input Parameters:
419: +  eps - the eigenproblem solver context
420: -  doubleexp - the boolean flag

422:    Options Database Keys:
423: .  -eps_gd_double_expansion - activate the double-expansion variant of GD

425:    Notes:
426:    In the double expansion variant the search subspace is expanded with K*[A*x B*x]
427:    instead of the classic K*r, where K is the preconditioner, x the selected
428:    approximate eigenvector and r its associated residual vector.

430:    Level: advanced

432: .seealso: EPSGDGetDoubleExpansion()
433: @*/
434: PetscErrorCode EPSGDSetDoubleExpansion(EPS eps,PetscBool doubleexp)
435: {
436:   PetscFunctionBegin;
439:   PetscTryMethod(eps,"EPSGDSetDoubleExpansion_C",(EPS,PetscBool),(eps,doubleexp));
440:   PetscFunctionReturn(PETSC_SUCCESS);
441: }

443: static PetscErrorCode EPSGDGetDoubleExpansion_GD(EPS eps,PetscBool *doubleexp)
444: {
445:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

447:   PetscFunctionBegin;
448:   *doubleexp = data->doubleexp;
449:   PetscFunctionReturn(PETSC_SUCCESS);
450: }

452: /*@
453:    EPSGDGetDoubleExpansion - Gets a flag indicating whether the double
454:    expansion variant has been activated or not.

456:    Not Collective

458:    Input Parameter:
459: .  eps - the eigenproblem solver context

461:    Output Parameter:
462: .  doubleexp - the flag

464:    Level: advanced

466: .seealso: EPSGDSetDoubleExpansion()
467: @*/
468: PetscErrorCode EPSGDGetDoubleExpansion(EPS eps,PetscBool *doubleexp)
469: {
470:   PetscFunctionBegin;
472:   PetscAssertPointer(doubleexp,2);
473:   PetscUseMethod(eps,"EPSGDGetDoubleExpansion_C",(EPS,PetscBool*),(eps,doubleexp));
474:   PetscFunctionReturn(PETSC_SUCCESS);
475: }

477: SLEPC_EXTERN PetscErrorCode EPSCreate_GD(EPS eps)
478: {
479:   EPS_DAVIDSON    *data;

481:   PetscFunctionBegin;
482:   PetscCall(PetscNew(&data));
483:   eps->data = (void*)data;

485:   data->blocksize   = 1;
486:   data->initialsize = 0;
487:   data->minv        = 0;
488:   data->plusk       = PETSC_DEFAULT;
489:   data->ipB         = PETSC_TRUE;
490:   data->fix         = 0.0;
491:   data->krylovstart = PETSC_FALSE;
492:   data->dynamic     = PETSC_FALSE;

494:   eps->useds = PETSC_TRUE;
495:   eps->categ = EPS_CATEGORY_PRECOND;

497:   eps->ops->solve          = EPSSolve_XD;
498:   eps->ops->setup          = EPSSetUp_GD;
499:   eps->ops->setupsort      = EPSSetUpSort_Default;
500:   eps->ops->setfromoptions = EPSSetFromOptions_GD;
501:   eps->ops->destroy        = EPSDestroy_GD;
502:   eps->ops->reset          = EPSReset_XD;
503:   eps->ops->view           = EPSView_GD;
504:   eps->ops->backtransform  = EPSBackTransform_Default;
505:   eps->ops->computevectors = EPSComputeVectors_XD;
506:   eps->ops->setdefaultst   = EPSSetDefaultST_Precond;

508:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",EPSXDSetKrylovStart_XD));
509:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",EPSXDGetKrylovStart_XD));
510:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",EPSXDSetBOrth_XD));
511:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",EPSXDGetBOrth_XD));
512:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",EPSXDSetBlockSize_XD));
513:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",EPSXDGetBlockSize_XD));
514:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",EPSXDSetRestart_XD));
515:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",EPSXDGetRestart_XD));
516:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",EPSXDSetInitialSize_XD));
517:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",EPSXDGetInitialSize_XD));
518:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",EPSGDSetDoubleExpansion_GD));
519:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",EPSGDGetDoubleExpansion_GD));
520:   PetscFunctionReturn(PETSC_SUCCESS);
521: }