Actual source code: gd.c
slepc-main 2024-12-17
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: Note:
250: PETSC_CURRENT can be used to preserve the current value of any of the
251: arguments, and PETSC_DETERMINE to set them to a default value.
253: Level: advanced
255: .seealso: EPSGDGetRestart()
256: @*/
257: PetscErrorCode EPSGDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)
258: {
259: PetscFunctionBegin;
263: PetscTryMethod(eps,"EPSGDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: /*@
268: EPSGDGetRestart - Gets the number of vectors of the searching space after
269: restarting and the number of vectors saved from the previous iteration.
271: Not Collective
273: Input Parameter:
274: . eps - the eigenproblem solver context
276: Output Parameters:
277: + minv - number of vectors of the searching subspace after restarting
278: - plusk - number of vectors saved from the previous iteration
280: Level: advanced
282: .seealso: EPSGDSetRestart()
283: @*/
284: PetscErrorCode EPSGDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)
285: {
286: PetscFunctionBegin;
288: PetscUseMethod(eps,"EPSGDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /*@
293: EPSGDSetInitialSize - Sets the initial size of the searching space.
295: Logically Collective
297: Input Parameters:
298: + eps - the eigenproblem solver context
299: - initialsize - number of vectors of the initial searching subspace
301: Options Database Key:
302: . -eps_gd_initial_size - number of vectors of the initial searching subspace
304: Notes:
305: If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
306: EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
307: provided vectors are not enough, the solver completes the subspace with
308: random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
309: gets the first vector provided by the user or, if not available, a random vector,
310: and expands the Krylov basis up to initialsize vectors.
312: Level: advanced
314: .seealso: EPSGDGetInitialSize(), EPSGDGetKrylovStart()
315: @*/
316: PetscErrorCode EPSGDSetInitialSize(EPS eps,PetscInt initialsize)
317: {
318: PetscFunctionBegin;
321: PetscTryMethod(eps,"EPSGDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*@
326: EPSGDGetInitialSize - Returns the initial size of the searching space.
328: Not Collective
330: Input Parameter:
331: . eps - the eigenproblem solver context
333: Output Parameter:
334: . initialsize - number of vectors of the initial searching subspace
336: Notes:
337: If EPSGDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
338: EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
339: provided vectors are not enough, the solver completes the subspace with
340: random vectors. In the case of EPSGDGetKrylovStart() being PETSC_TRUE, the solver
341: gets the first vector provided by the user or, if not available, a random vector,
342: and expands the Krylov basis up to initialsize vectors.
344: Level: advanced
346: .seealso: 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 eigenproblem solver context
365: - borth - whether to B-orthogonalize the search subspace
367: Options Database Key:
368: . -eps_gd_borth - Set the orthogonalization used in the search subspace
370: Level: advanced
372: .seealso: EPSGDGetBOrth()
373: @*/
374: PetscErrorCode EPSGDSetBOrth(EPS eps,PetscBool borth)
375: {
376: PetscFunctionBegin;
379: PetscTryMethod(eps,"EPSGDSetBOrth_C",(EPS,PetscBool),(eps,borth));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: /*@
384: EPSGDGetBOrth - Returns the orthogonalization used in the search
385: subspace in case of generalized Hermitian problems.
387: Not Collective
389: Input Parameter:
390: . eps - the eigenproblem solver context
392: Output Parameters:
393: . borth - whether to B-orthogonalize the search subspace
395: Level: advanced
397: .seealso: EPSGDSetBOrth()
398: @*/
399: PetscErrorCode EPSGDGetBOrth(EPS eps,PetscBool *borth)
400: {
401: PetscFunctionBegin;
403: PetscAssertPointer(borth,2);
404: PetscUseMethod(eps,"EPSGDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode EPSGDSetDoubleExpansion_GD(EPS eps,PetscBool doubleexp)
409: {
410: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
412: PetscFunctionBegin;
413: data->doubleexp = doubleexp;
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: /*@
418: EPSGDSetDoubleExpansion - Activate the double expansion variant of GD.
420: Logically Collective
422: Input Parameters:
423: + eps - the eigenproblem solver context
424: - doubleexp - the boolean flag
426: Options Database Keys:
427: . -eps_gd_double_expansion - activate the double-expansion variant of GD
429: Notes:
430: In the double expansion variant the search subspace is expanded with K*[A*x B*x]
431: instead of the classic K*r, where K is the preconditioner, x the selected
432: approximate eigenvector and r its associated residual vector.
434: Level: advanced
436: .seealso: EPSGDGetDoubleExpansion()
437: @*/
438: PetscErrorCode EPSGDSetDoubleExpansion(EPS eps,PetscBool doubleexp)
439: {
440: PetscFunctionBegin;
443: PetscTryMethod(eps,"EPSGDSetDoubleExpansion_C",(EPS,PetscBool),(eps,doubleexp));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode EPSGDGetDoubleExpansion_GD(EPS eps,PetscBool *doubleexp)
448: {
449: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
451: PetscFunctionBegin;
452: *doubleexp = data->doubleexp;
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: /*@
457: EPSGDGetDoubleExpansion - Gets a flag indicating whether the double
458: expansion variant has been activated or not.
460: Not Collective
462: Input Parameter:
463: . eps - the eigenproblem solver context
465: Output Parameter:
466: . doubleexp - the flag
468: Level: advanced
470: .seealso: EPSGDSetDoubleExpansion()
471: @*/
472: PetscErrorCode EPSGDGetDoubleExpansion(EPS eps,PetscBool *doubleexp)
473: {
474: PetscFunctionBegin;
476: PetscAssertPointer(doubleexp,2);
477: PetscUseMethod(eps,"EPSGDGetDoubleExpansion_C",(EPS,PetscBool*),(eps,doubleexp));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: SLEPC_EXTERN PetscErrorCode EPSCreate_GD(EPS eps)
482: {
483: EPS_DAVIDSON *data;
485: PetscFunctionBegin;
486: PetscCall(PetscNew(&data));
487: eps->data = (void*)data;
489: data->blocksize = 1;
490: data->initialsize = 0;
491: data->minv = 0;
492: data->plusk = PETSC_DETERMINE;
493: data->ipB = PETSC_TRUE;
494: data->fix = 0.0;
495: data->krylovstart = PETSC_FALSE;
496: data->dynamic = PETSC_FALSE;
498: eps->useds = PETSC_TRUE;
499: eps->categ = EPS_CATEGORY_PRECOND;
501: eps->ops->solve = EPSSolve_XD;
502: eps->ops->setup = EPSSetUp_GD;
503: eps->ops->setupsort = EPSSetUpSort_Default;
504: eps->ops->setfromoptions = EPSSetFromOptions_GD;
505: eps->ops->destroy = EPSDestroy_GD;
506: eps->ops->reset = EPSReset_XD;
507: eps->ops->view = EPSView_GD;
508: eps->ops->backtransform = EPSBackTransform_Default;
509: eps->ops->computevectors = EPSComputeVectors_XD;
510: eps->ops->setdefaultst = EPSSetDefaultST_Precond;
512: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetKrylovStart_C",EPSXDSetKrylovStart_XD));
513: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetKrylovStart_C",EPSXDGetKrylovStart_XD));
514: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBOrth_C",EPSXDSetBOrth_XD));
515: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBOrth_C",EPSXDGetBOrth_XD));
516: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetBlockSize_C",EPSXDSetBlockSize_XD));
517: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetBlockSize_C",EPSXDGetBlockSize_XD));
518: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetRestart_C",EPSXDSetRestart_XD));
519: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetRestart_C",EPSXDGetRestart_XD));
520: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetInitialSize_C",EPSXDSetInitialSize_XD));
521: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetInitialSize_C",EPSXDGetInitialSize_XD));
522: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDSetDoubleExpansion_C",EPSGDSetDoubleExpansion_GD));
523: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSGDGetDoubleExpansion_C",EPSGDGetDoubleExpansion_GD));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }