Actual source code: gd.c
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 linear eigensolver context
138: - krylovstart - boolean flag
140: Options Database Key:
141: . -eps_gd_krylov_start - activate starting the searching subspace with a Krylov basis
143: Note:
144: See discussion at `EPSGDSetInitialSize()`.
146: Level: advanced
148: .seealso: [](ch:eps), `EPSGD`, `EPSGDGetKrylovStart()`, `EPSGDSetInitialSize()`
149: @*/
150: PetscErrorCode EPSGDSetKrylovStart(EPS eps,PetscBool krylovstart)
151: {
152: PetscFunctionBegin;
155: PetscTryMethod(eps,"EPSGDSetKrylovStart_C",(EPS,PetscBool),(eps,krylovstart));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@
160: EPSGDGetKrylovStart - Returns a flag indicating if the search subspace is started with a
161: Krylov basis.
163: Not Collective
165: Input Parameter:
166: . eps - the linear eigensolver context
168: Output Parameter:
169: . krylovstart - boolean flag indicating if the search subspace is started with a Krylov basis
171: Level: advanced
173: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetKrylovStart()`
174: @*/
175: PetscErrorCode EPSGDGetKrylovStart(EPS eps,PetscBool *krylovstart)
176: {
177: PetscFunctionBegin;
179: PetscAssertPointer(krylovstart,2);
180: PetscUseMethod(eps,"EPSGDGetKrylovStart_C",(EPS,PetscBool*),(eps,krylovstart));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@
185: EPSGDSetBlockSize - Sets the number of vectors to be added to the searching space
186: in every iteration.
188: Logically Collective
190: Input Parameters:
191: + eps - the linear eigensolver context
192: - blocksize - number of vectors added to the search space in every iteration
194: Options Database Key:
195: . -eps_gd_blocksize \<blocksize\> - number of vectors added to the search space in every iteration
197: Note:
198: Detailed information can be found at {cite:p}`Rom14`.
200: Level: advanced
202: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetKrylovStart()`
203: @*/
204: PetscErrorCode EPSGDSetBlockSize(EPS eps,PetscInt blocksize)
205: {
206: PetscFunctionBegin;
209: PetscTryMethod(eps,"EPSGDSetBlockSize_C",(EPS,PetscInt),(eps,blocksize));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@
214: EPSGDGetBlockSize - Returns the number of vectors to be added to the searching space
215: in every iteration.
217: Not Collective
219: Input Parameter:
220: . eps - the linear eigensolver context
222: Output Parameter:
223: . blocksize - number of vectors added to the search space in every iteration
225: Level: advanced
227: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetBlockSize()`
228: @*/
229: PetscErrorCode EPSGDGetBlockSize(EPS eps,PetscInt *blocksize)
230: {
231: PetscFunctionBegin;
233: PetscAssertPointer(blocksize,2);
234: PetscUseMethod(eps,"EPSGDGetBlockSize_C",(EPS,PetscInt*),(eps,blocksize));
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: EPSGDSetRestart - Sets the number of vectors of the searching space after
240: restarting and the number of vectors saved from the previous iteration.
242: Logically Collective
244: Input Parameters:
245: + eps - the linear eigensolver context
246: . minv - number of vectors of the searching subspace after restarting
247: - plusk - number of vectors saved from the previous iteration
249: Options Database Keys:
250: + -eps_gd_minv \<minv\> - number of vectors of the searching subspace after restarting
251: - -eps_gd_plusk \<plusk\> - number of vectors saved from the previous iteration
253: Notes:
254: `PETSC_CURRENT` can be used to preserve the current value of any of the
255: arguments, and `PETSC_DETERMINE` to set them to a default value.
257: Detailed information can be found at {cite:p}`Rom14`.
259: Level: advanced
261: .seealso: [](ch:eps), `EPSGD`, `EPSGDGetRestart()`
262: @*/
263: PetscErrorCode EPSGDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)
264: {
265: PetscFunctionBegin;
269: PetscTryMethod(eps,"EPSGDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: EPSGDGetRestart - Gets the number of vectors of the searching space after
275: restarting and the number of vectors saved from the previous iteration.
277: Not Collective
279: Input Parameter:
280: . eps - the linear eigensolver context
282: Output Parameters:
283: + minv - number of vectors of the searching subspace after restarting
284: - plusk - number of vectors saved from the previous iteration
286: Level: advanced
288: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetRestart()`
289: @*/
290: PetscErrorCode EPSGDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)
291: {
292: PetscFunctionBegin;
294: PetscUseMethod(eps,"EPSGDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*@
299: EPSGDSetInitialSize - Sets the initial size of the searching space.
301: Logically Collective
303: Input Parameters:
304: + eps - the linear eigensolver context
305: - initialsize - number of vectors of the initial searching subspace
307: Options Database Key:
308: . -eps_gd_initial_size \<initialsize\> - number of vectors of the initial searching subspace
310: Notes:
311: If the flag in `EPSGDSetKrylovStart()` is set to `PETSC_FALSE` and the user
312: provides vectors with `EPSSetInitialSpace()`, up to `initialsize` vectors will be used;
313: and if the provided vectors are not enough, the solver completes the subspace with
314: random vectors. In case the `EPSGDSetKrylovStart()` flag is `PETSC_TRUE`, the solver
315: gets the first vector provided by the user or, if not available, a random vector,
316: and expands the Krylov basis up to `initialsize` vectors.
318: Detailed information can be found at {cite:p}`Rom14`.
320: Level: advanced
322: .seealso: [](ch:eps), `EPSGD`, `EPSGDGetInitialSize()`, `EPSGDSetKrylovStart()`
323: @*/
324: PetscErrorCode EPSGDSetInitialSize(EPS eps,PetscInt initialsize)
325: {
326: PetscFunctionBegin;
329: PetscTryMethod(eps,"EPSGDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /*@
334: EPSGDGetInitialSize - Returns the initial size of the searching space.
336: Not Collective
338: Input Parameter:
339: . eps - the linear eigensolver context
341: Output Parameter:
342: . initialsize - number of vectors of the initial searching subspace
344: Level: advanced
346: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetInitialSize()`, `EPSGDGetKrylovStart()`
347: @*/
348: PetscErrorCode EPSGDGetInitialSize(EPS eps,PetscInt *initialsize)
349: {
350: PetscFunctionBegin;
352: PetscAssertPointer(initialsize,2);
353: PetscUseMethod(eps,"EPSGDGetInitialSize_C",(EPS,PetscInt*),(eps,initialsize));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@
358: EPSGDSetBOrth - Selects the orthogonalization that will be used in the search
359: subspace in case of generalized Hermitian problems.
361: Logically Collective
363: Input Parameters:
364: + eps - the linear eigensolver context
365: - borth - whether to $B$-orthogonalize the search subspace
367: Options Database Key:
368: . -eps_gd_borth - toggle the $B$-orthogonalization
370: Note:
371: Detailed information can be found at {cite:p}`Rom14`.
373: Level: advanced
375: .seealso: [](ch:eps), `EPSGD`, `EPSGDGetBOrth()`
376: @*/
377: PetscErrorCode EPSGDSetBOrth(EPS eps,PetscBool borth)
378: {
379: PetscFunctionBegin;
382: PetscTryMethod(eps,"EPSGDSetBOrth_C",(EPS,PetscBool),(eps,borth));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: EPSGDGetBOrth - Returns the orthogonalization used in the search
388: subspace in case of generalized Hermitian problems.
390: Not Collective
392: Input Parameter:
393: . eps - the linear eigensolver context
395: Output Parameter:
396: . borth - whether to $B$-orthogonalize the search subspace
398: Level: advanced
400: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetBOrth()`
401: @*/
402: PetscErrorCode EPSGDGetBOrth(EPS eps,PetscBool *borth)
403: {
404: PetscFunctionBegin;
406: PetscAssertPointer(borth,2);
407: PetscUseMethod(eps,"EPSGDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: static PetscErrorCode EPSGDSetDoubleExpansion_GD(EPS eps,PetscBool doubleexp)
412: {
413: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
415: PetscFunctionBegin;
416: data->doubleexp = doubleexp;
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: /*@
421: EPSGDSetDoubleExpansion - Activate the double expansion variant of GD.
423: Logically Collective
425: Input Parameters:
426: + eps - the linear eigensolver context
427: - doubleexp - the boolean flag
429: Options Database Key:
430: . -eps_gd_double_expansion - activate the double-expansion variant of GD
432: Note:
433: In the double expansion variant the search subspace is expanded with $K[Ax, Bx]$
434: instead of the classic $Kr$, where $K$ is the preconditioner, $x$ the selected
435: approximate eigenvector and $r$ its associated residual vector. More details
436: can be found in {cite:p}`Hoc12`.
438: Level: advanced
440: .seealso: [](ch:eps), `EPSGD`, `EPSGDGetDoubleExpansion()`
441: @*/
442: PetscErrorCode EPSGDSetDoubleExpansion(EPS eps,PetscBool doubleexp)
443: {
444: PetscFunctionBegin;
447: PetscTryMethod(eps,"EPSGDSetDoubleExpansion_C",(EPS,PetscBool),(eps,doubleexp));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode EPSGDGetDoubleExpansion_GD(EPS eps,PetscBool *doubleexp)
452: {
453: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
455: PetscFunctionBegin;
456: *doubleexp = data->doubleexp;
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@
461: EPSGDGetDoubleExpansion - Gets a flag indicating whether the double
462: expansion variant has been activated or not.
464: Not Collective
466: Input Parameter:
467: . eps - the linear eigensolver context
469: Output Parameter:
470: . doubleexp - the flag
472: Level: advanced
474: .seealso: [](ch:eps), `EPSGD`, `EPSGDSetDoubleExpansion()`
475: @*/
476: PetscErrorCode EPSGDGetDoubleExpansion(EPS eps,PetscBool *doubleexp)
477: {
478: PetscFunctionBegin;
480: PetscAssertPointer(doubleexp,2);
481: PetscUseMethod(eps,"EPSGDGetDoubleExpansion_C",(EPS,PetscBool*),(eps,doubleexp));
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: /*MC
486: EPSGD - EPSGD = "gd" - The generalized Davidson method.
488: Notes:
489: This is a preconditioned eigensolver, that is, it may be competitive
490: when computing interior eigenvalues in case the shift-and-invert spectral
491: transformation is too costly and a good preconditioner is available.
493: The implemented method is Generalized Davidson, where any `PC`
494: preconditioner can be used, in contrast to the original Davidson method
495: that used a diagonal preconditioner.
497: The preconditioner is specified via the internal `ST` object and its
498: associated `KSP`.
500: Details of the implementation are described in {cite:p}`Rom14`.
502: Level: beginner
504: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSGetST()`
505: M*/
506: SLEPC_EXTERN PetscErrorCode EPSCreate_GD(EPS eps)
507: {
508: EPS_DAVIDSON *data;
510: PetscFunctionBegin;
511: PetscCall(PetscNew(&data));
512: eps->data = (void*)data;
514: data->blocksize = 1;
515: data->initialsize = 0;
516: data->minv = 0;
517: data->plusk = PETSC_DETERMINE;
518: data->ipB = PETSC_TRUE;
519: data->fix = 0.0;
520: data->krylovstart = PETSC_FALSE;
521: data->dynamic = PETSC_FALSE;
523: eps->useds = PETSC_TRUE;
524: eps->categ = EPS_CATEGORY_PRECOND;
526: eps->ops->solve = EPSSolve_XD;
527: eps->ops->setup = EPSSetUp_GD;
528: eps->ops->setupsort = EPSSetUpSort_Default;
529: eps->ops->setfromoptions = EPSSetFromOptions_GD;
530: eps->ops->destroy = EPSDestroy_GD;
531: eps->ops->reset = EPSReset_XD;
532: eps->ops->view = EPSView_GD;
533: eps->ops->backtransform = EPSBackTransform_Default;
534: eps->ops->computevectors = EPSComputeVectors_XD;
535: eps->ops->setdefaultst = EPSSetDefaultST_Precond;
537: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",EPSXDSetKrylovStart_XD));
538: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",EPSXDGetKrylovStart_XD));
539: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",EPSXDSetBOrth_XD));
540: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",EPSXDGetBOrth_XD));
541: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",EPSXDSetBlockSize_XD));
542: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",EPSXDGetBlockSize_XD));
543: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",EPSXDSetRestart_XD));
544: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",EPSXDGetRestart_XD));
545: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",EPSXDSetInitialSize_XD));
546: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",EPSXDGetInitialSize_XD));
547: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",EPSGDSetDoubleExpansion_GD));
548: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",EPSGDGetDoubleExpansion_GD));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }