Actual source code: gd.c
slepc-3.20.2 2024-03-15
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: }