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: "jd"
13: Method: Jacobi-Davidson
15: Algorithm:
17: Jacobi-Davidson with various subspace extraction and
18: restart techniques.
20: References:
22: [1] G.L.G. Sleijpen and H.A. van der Vorst, "A Jacobi-Davidson
23: iteration method for linear eigenvalue problems", SIAM J.
24: Matrix Anal. Appl. 17(2):401-425, 1996.
26: [2] E. Romero and J.E. Roman, "A parallel implementation of
27: Davidson methods for large-scale eigenvalue problems in
28: SLEPc", ACM Trans. Math. Software 40(2), Article 13, 2014.
29: */
31: #include <slepc/private/epsimpl.h> 32: #include <../src/eps/impls/davidson/davidson.h> 34: PetscErrorCode EPSSetFromOptions_JD(PetscOptionItems *PetscOptionsObject,EPS eps) 35: {
36: PetscBool flg,flg2,op,orth;
37: PetscInt opi,opi0;
38: PetscReal opf;
40: PetscOptionsHead(PetscOptionsObject,"EPS Jacobi-Davidson (JD) Options");
42: EPSJDGetKrylovStart(eps,&op);
43: PetscOptionsBool("-eps_jd_krylov_start","Start the search subspace with a Krylov basis","EPSJDSetKrylovStart",op,&op,&flg);
44: if (flg) EPSJDSetKrylovStart(eps,op);
46: EPSJDGetBOrth(eps,&orth);
47: PetscOptionsBool("-eps_jd_borth","Use B-orthogonalization in the search subspace","EPSJDSetBOrth",op,&op,&flg);
48: if (flg) EPSJDSetBOrth(eps,op);
50: EPSJDGetBlockSize(eps,&opi);
51: PetscOptionsInt("-eps_jd_blocksize","Number of vectors to add to the search subspace","EPSJDSetBlockSize",opi,&opi,&flg);
52: if (flg) EPSJDSetBlockSize(eps,opi);
54: EPSJDGetRestart(eps,&opi,&opi0);
55: PetscOptionsInt("-eps_jd_minv","Size of the search subspace after restarting","EPSJDSetRestart",opi,&opi,&flg);
56: PetscOptionsInt("-eps_jd_plusk","Number of eigenvectors saved from the previous iteration when restarting","EPSJDSetRestart",opi0,&opi0,&flg2);
57: if (flg || flg2) EPSJDSetRestart(eps,opi,opi0);
59: EPSJDGetInitialSize(eps,&opi);
60: PetscOptionsInt("-eps_jd_initial_size","Initial size of the search subspace","EPSJDSetInitialSize",opi,&opi,&flg);
61: if (flg) EPSJDSetInitialSize(eps,opi);
63: EPSJDGetFix(eps,&opf);
64: PetscOptionsReal("-eps_jd_fix","Tolerance for changing the target in the correction equation","EPSJDSetFix",opf,&opf,&flg);
65: if (flg) EPSJDSetFix(eps,opf);
67: EPSJDGetConstCorrectionTol(eps,&op);
68: PetscOptionsBool("-eps_jd_const_correction_tol","Disable the dynamic stopping criterion when solving the correction equation","EPSJDSetConstCorrectionTol",op,&op,&flg);
69: if (flg) EPSJDSetConstCorrectionTol(eps,op);
71: PetscOptionsTail();
72: PetscFunctionReturn(0);
73: }
75: PetscErrorCode EPSSetDefaultST_JD(EPS eps) 76: {
77: KSP ksp;
79: if (!((PetscObject)eps->st)->type_name) {
80: STSetType(eps->st,STPRECOND);
81: STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
82: }
83: STGetKSP(eps->st,&ksp);
84: if (!((PetscObject)ksp)->type_name) {
85: KSPSetType(ksp,KSPBCGSL);
86: KSPSetTolerances(ksp,1e-4,PETSC_DEFAULT,PETSC_DEFAULT,90);
87: }
88: PetscFunctionReturn(0);
89: }
91: PetscErrorCode EPSSetUp_JD(EPS eps) 92: {
93: PetscBool t;
94: KSP ksp;
96: /* Setup common for all davidson solvers */
97: EPSSetUp_XD(eps);
99: /* Check some constraints */
100: STGetKSP(eps->st,&ksp);
101: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&t);
103: PetscFunctionReturn(0);
104: }
106: PetscErrorCode EPSView_JD(EPS eps,PetscViewer viewer)107: {
108: PetscBool isascii,opb;
109: PetscReal opf;
110: PetscInt opi,opi0;
112: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
113: if (isascii) {
114: EPSXDGetBOrth_XD(eps,&opb);
115: if (opb) PetscViewerASCIIPrintf(viewer," search subspace is B-orthogonalized\n");
116: else PetscViewerASCIIPrintf(viewer," search subspace is orthogonalized\n");
117: EPSXDGetBlockSize_XD(eps,&opi);
118: PetscViewerASCIIPrintf(viewer," block size=%" PetscInt_FMT "\n",opi);
119: EPSXDGetKrylovStart_XD(eps,&opb);
120: if (!opb) PetscViewerASCIIPrintf(viewer," type of the initial subspace: non-Krylov\n");
121: else PetscViewerASCIIPrintf(viewer," type of the initial subspace: Krylov\n");
122: EPSXDGetRestart_XD(eps,&opi,&opi0);
123: PetscViewerASCIIPrintf(viewer," size of the subspace after restarting: %" PetscInt_FMT "\n",opi);
124: PetscViewerASCIIPrintf(viewer," number of vectors after restarting from the previous iteration: %" PetscInt_FMT "\n",opi0);
126: EPSJDGetFix_JD(eps,&opf);
127: PetscViewerASCIIPrintf(viewer," threshold for changing the target in the correction equation (fix): %g\n",(double)opf);
129: EPSJDGetConstCorrectionTol_JD(eps,&opb);
130: if (!opb) PetscViewerASCIIPrintf(viewer," using dynamic tolerance for the correction equation\n");
131: }
132: PetscFunctionReturn(0);
133: }
135: PetscErrorCode EPSDestroy_JD(EPS eps)136: {
137: PetscFree(eps->data);
138: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetKrylovStart_C",NULL);
139: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetKrylovStart_C",NULL);
140: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBlockSize_C",NULL);
141: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBlockSize_C",NULL);
142: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetRestart_C",NULL);
143: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetRestart_C",NULL);
144: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetInitialSize_C",NULL);
145: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetInitialSize_C",NULL);
146: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetFix_C",NULL);
147: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetFix_C",NULL);
148: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetConstCorrectionTol_C",NULL);
149: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetConstCorrectionTol_C",NULL);
150: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBOrth_C",NULL);
151: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBOrth_C",NULL);
152: PetscFunctionReturn(0);
153: }
155: /*@
156: EPSJDSetKrylovStart - Activates or deactivates starting the searching
157: subspace with a Krylov basis.
159: Logically Collective on eps
161: Input Parameters:
162: + eps - the eigenproblem solver context
163: - krylovstart - boolean flag
165: Options Database Key:
166: . -eps_jd_krylov_start - Activates starting the searching subspace with a
167: Krylov basis
169: Level: advanced
171: .seealso: EPSJDGetKrylovStart()
172: @*/
173: PetscErrorCode EPSJDSetKrylovStart(EPS eps,PetscBool krylovstart)174: {
177: PetscTryMethod(eps,"EPSJDSetKrylovStart_C",(EPS,PetscBool),(eps,krylovstart));
178: PetscFunctionReturn(0);
179: }
181: /*@
182: EPSJDGetKrylovStart - Returns a flag indicating if the searching subspace is started with a
183: Krylov basis.
185: Not Collective
187: Input Parameter:
188: . eps - the eigenproblem solver context
190: Output Parameters:
191: . krylovstart - boolean flag indicating if the searching subspace is started
192: with a Krylov basis
194: Level: advanced
196: .seealso: EPSJDSetKrylovStart()
197: @*/
198: PetscErrorCode EPSJDGetKrylovStart(EPS eps,PetscBool *krylovstart)199: {
202: PetscUseMethod(eps,"EPSJDGetKrylovStart_C",(EPS,PetscBool*),(eps,krylovstart));
203: PetscFunctionReturn(0);
204: }
206: /*@
207: EPSJDSetBlockSize - Sets the number of vectors to be added to the searching space
208: in every iteration.
210: Logically Collective on eps
212: Input Parameters:
213: + eps - the eigenproblem solver context
214: - blocksize - number of vectors added to the search space in every iteration
216: Options Database Key:
217: . -eps_jd_blocksize - number of vectors added to the searching space every iteration
219: Level: advanced
221: .seealso: EPSJDSetKrylovStart()
222: @*/
223: PetscErrorCode EPSJDSetBlockSize(EPS eps,PetscInt blocksize)224: {
227: PetscTryMethod(eps,"EPSJDSetBlockSize_C",(EPS,PetscInt),(eps,blocksize));
228: PetscFunctionReturn(0);
229: }
231: /*@
232: EPSJDGetBlockSize - Returns the number of vectors to be added to the searching space
233: in every iteration.
235: Not Collective
237: Input Parameter:
238: . eps - the eigenproblem solver context
240: Output Parameter:
241: . blocksize - number of vectors added to the search space in every iteration
243: Level: advanced
245: .seealso: EPSJDSetBlockSize()
246: @*/
247: PetscErrorCode EPSJDGetBlockSize(EPS eps,PetscInt *blocksize)248: {
251: PetscUseMethod(eps,"EPSJDGetBlockSize_C",(EPS,PetscInt*),(eps,blocksize));
252: PetscFunctionReturn(0);
253: }
255: /*@
256: EPSJDSetRestart - Sets the number of vectors of the searching space after
257: restarting and the number of vectors saved from the previous iteration.
259: Logically Collective on eps
261: Input Parameters:
262: + eps - the eigenproblem solver context
263: . minv - number of vectors of the searching subspace after restarting
264: - plusk - number of vectors saved from the previous iteration
266: Options Database Keys:
267: + -eps_jd_minv - number of vectors of the searching subspace after restarting
268: - -eps_jd_plusk - number of vectors saved from the previous iteration
270: Level: advanced
272: .seealso: EPSJDGetRestart()
273: @*/
274: PetscErrorCode EPSJDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)275: {
279: PetscTryMethod(eps,"EPSJDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
280: PetscFunctionReturn(0);
281: }
283: /*@
284: EPSJDGetRestart - Gets the number of vectors of the searching space after
285: restarting and the number of vectors saved from the previous iteration.
287: Not Collective
289: Input Parameter:
290: . eps - the eigenproblem solver context
292: Output Parameters:
293: + minv - number of vectors of the searching subspace after restarting
294: - plusk - number of vectors saved from the previous iteration
296: Level: advanced
298: .seealso: EPSJDSetRestart()
299: @*/
300: PetscErrorCode EPSJDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)301: {
303: PetscUseMethod(eps,"EPSJDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
304: PetscFunctionReturn(0);
305: }
307: /*@
308: EPSJDSetInitialSize - Sets the initial size of the searching space.
310: Logically Collective on eps
312: Input Parameters:
313: + eps - the eigenproblem solver context
314: - initialsize - number of vectors of the initial searching subspace
316: Options Database Key:
317: . -eps_jd_initial_size - number of vectors of the initial searching subspace
319: Notes:
320: If EPSJDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
321: EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
322: provided vectors are not enough, the solver completes the subspace with
323: random vectors. In the case of EPSJDGetKrylovStart() being PETSC_TRUE, the solver
324: gets the first vector provided by the user or, if not available, a random vector,
325: and expands the Krylov basis up to initialsize vectors.
327: Level: advanced
329: .seealso: EPSJDGetInitialSize(), EPSJDGetKrylovStart()
330: @*/
331: PetscErrorCode EPSJDSetInitialSize(EPS eps,PetscInt initialsize)332: {
335: PetscTryMethod(eps,"EPSJDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
336: PetscFunctionReturn(0);
337: }
339: /*@
340: EPSJDGetInitialSize - Returns the initial size of the searching space.
342: Not Collective
344: Input Parameter:
345: . eps - the eigenproblem solver context
347: Output Parameter:
348: . initialsize - number of vectors of the initial searching subspace
350: Notes:
351: If EPSJDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
352: EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
353: provided vectors are not enough, the solver completes the subspace with
354: random vectors. In the case of EPSJDGetKrylovStart() being PETSC_TRUE, the solver
355: gets the first vector provided by the user or, if not available, a random vector,
356: and expands the Krylov basis up to initialsize vectors.
358: Level: advanced
360: .seealso: EPSJDSetInitialSize(), EPSJDGetKrylovStart()
361: @*/
362: PetscErrorCode EPSJDGetInitialSize(EPS eps,PetscInt *initialsize)363: {
366: PetscUseMethod(eps,"EPSJDGetInitialSize_C",(EPS,PetscInt*),(eps,initialsize));
367: PetscFunctionReturn(0);
368: }
370: PetscErrorCode EPSJDSetFix_JD(EPS eps,PetscReal fix)371: {
372: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
374: if (fix == PETSC_DEFAULT || fix == PETSC_DECIDE) fix = 0.01;
376: data->fix = fix;
377: PetscFunctionReturn(0);
378: }
380: /*@
381: EPSJDSetFix - Sets the threshold for changing the target in the correction
382: equation.
384: Logically Collective on eps
386: Input Parameters:
387: + eps - the eigenproblem solver context
388: - fix - threshold for changing the target
390: Options Database Key:
391: . -eps_jd_fix - the fix value
393: Note:
394: The target in the correction equation is fixed at the first iterations.
395: When the norm of the residual vector is lower than the fix value,
396: the target is set to the corresponding eigenvalue.
398: Level: advanced
400: .seealso: EPSJDGetFix()
401: @*/
402: PetscErrorCode EPSJDSetFix(EPS eps,PetscReal fix)403: {
406: PetscTryMethod(eps,"EPSJDSetFix_C",(EPS,PetscReal),(eps,fix));
407: PetscFunctionReturn(0);
408: }
410: PetscErrorCode EPSJDGetFix_JD(EPS eps,PetscReal *fix)411: {
412: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
414: *fix = data->fix;
415: PetscFunctionReturn(0);
416: }
418: /*@
419: EPSJDGetFix - Returns the threshold for changing the target in the correction
420: equation.
422: Not Collective
424: Input Parameter:
425: . eps - the eigenproblem solver context
427: Output Parameter:
428: . fix - threshold for changing the target
430: Note:
431: The target in the correction equation is fixed at the first iterations.
432: When the norm of the residual vector is lower than the fix value,
433: the target is set to the corresponding eigenvalue.
435: Level: advanced
437: .seealso: EPSJDSetFix()
438: @*/
439: PetscErrorCode EPSJDGetFix(EPS eps,PetscReal *fix)440: {
443: PetscUseMethod(eps,"EPSJDGetFix_C",(EPS,PetscReal*),(eps,fix));
444: PetscFunctionReturn(0);
445: }
447: PetscErrorCode EPSJDSetConstCorrectionTol_JD(EPS eps,PetscBool constant)448: {
449: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
451: data->dynamic = PetscNot(constant);
452: PetscFunctionReturn(0);
453: }
455: /*@
456: EPSJDSetConstCorrectionTol - If true, deactivates the dynamic stopping criterion
457: (also called Newton) that sets the KSP relative tolerance
458: to 0.5**i, where i is the number of EPS iterations from the last converged value.
460: Logically Collective on eps
462: Input Parameters:
463: + eps - the eigenproblem solver context
464: - constant - if false, the KSP relative tolerance is set to 0.5**i.
466: Options Database Key:
467: . -eps_jd_const_correction_tol - Deactivates the dynamic stopping criterion.
469: Level: advanced
471: .seealso: EPSJDGetConstCorrectionTol()
472: @*/
473: PetscErrorCode EPSJDSetConstCorrectionTol(EPS eps,PetscBool constant)474: {
477: PetscTryMethod(eps,"EPSJDSetConstCorrectionTol_C",(EPS,PetscBool),(eps,constant));
478: PetscFunctionReturn(0);
479: }
481: PetscErrorCode EPSJDGetConstCorrectionTol_JD(EPS eps,PetscBool *constant)482: {
483: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
485: *constant = PetscNot(data->dynamic);
486: PetscFunctionReturn(0);
487: }
489: /*@
490: EPSJDGetConstCorrectionTol - Returns a flag indicating if the dynamic stopping is being used for
491: solving the correction equation.
493: Not Collective
495: Input Parameter:
496: . eps - the eigenproblem solver context
498: Output Parameter:
499: . constant - boolean flag indicating if the dynamic stopping criterion is not being used.
501: Notes:
502: If the flag is false the KSP relative tolerance is set to 0.5**i, where i is the number
503: of EPS iterations from the last converged value.
505: Level: advanced
507: .seealso: EPSJDSetConstCorrectionTol()
508: @*/
509: PetscErrorCode EPSJDGetConstCorrectionTol(EPS eps,PetscBool *constant)510: {
513: PetscUseMethod(eps,"EPSJDGetConstCorrectionTol_C",(EPS,PetscBool*),(eps,constant));
514: PetscFunctionReturn(0);
515: }
517: /*@
518: EPSJDSetBOrth - Selects the orthogonalization that will be used in the search
519: subspace in case of generalized Hermitian problems.
521: Logically Collective on eps
523: Input Parameters:
524: + eps - the eigenproblem solver context
525: - borth - whether to B-orthogonalize the search subspace
527: Options Database Key:
528: . -eps_jd_borth - Set the orthogonalization used in the search subspace
530: Level: advanced
532: .seealso: EPSJDGetBOrth()
533: @*/
534: PetscErrorCode EPSJDSetBOrth(EPS eps,PetscBool borth)535: {
538: PetscTryMethod(eps,"EPSJDSetBOrth_C",(EPS,PetscBool),(eps,borth));
539: PetscFunctionReturn(0);
540: }
542: /*@
543: EPSJDGetBOrth - Returns the orthogonalization used in the search
544: subspace in case of generalized Hermitian problems.
546: Not Collective
548: Input Parameter:
549: . eps - the eigenproblem solver context
551: Output Parameters:
552: . borth - whether to B-orthogonalize the search subspace
554: Level: advanced
556: .seealso: EPSJDSetBOrth()
557: @*/
558: PetscErrorCode EPSJDGetBOrth(EPS eps,PetscBool *borth)559: {
562: PetscUseMethod(eps,"EPSJDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
563: PetscFunctionReturn(0);
564: }
566: SLEPC_EXTERN PetscErrorCode EPSCreate_JD(EPS eps)567: {
568: EPS_DAVIDSON *data;
570: PetscNewLog(eps,&data);
571: eps->data = (void*)data;
573: data->blocksize = 1;
574: data->initialsize = 0;
575: data->minv = 0;
576: data->plusk = PETSC_DEFAULT;
577: data->ipB = PETSC_TRUE;
578: data->fix = 0.01;
579: data->krylovstart = PETSC_FALSE;
580: data->dynamic = PETSC_FALSE;
582: eps->useds = PETSC_TRUE;
583: eps->categ = EPS_CATEGORY_PRECOND;
585: eps->ops->solve = EPSSolve_XD;
586: eps->ops->setup = EPSSetUp_JD;
587: eps->ops->setupsort = EPSSetUpSort_Default;
588: eps->ops->setfromoptions = EPSSetFromOptions_JD;
589: eps->ops->destroy = EPSDestroy_JD;
590: eps->ops->reset = EPSReset_XD;
591: eps->ops->view = EPSView_JD;
592: eps->ops->backtransform = EPSBackTransform_Default;
593: eps->ops->computevectors = EPSComputeVectors_XD;
594: eps->ops->setdefaultst = EPSSetDefaultST_JD;
596: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetKrylovStart_C",EPSXDSetKrylovStart_XD);
597: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetKrylovStart_C",EPSXDGetKrylovStart_XD);
598: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBlockSize_C",EPSXDSetBlockSize_XD);
599: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBlockSize_C",EPSXDGetBlockSize_XD);
600: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetRestart_C",EPSXDSetRestart_XD);
601: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetRestart_C",EPSXDGetRestart_XD);
602: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetInitialSize_C",EPSXDSetInitialSize_XD);
603: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetInitialSize_C",EPSXDGetInitialSize_XD);
604: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetFix_C",EPSJDSetFix_JD);
605: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetFix_C",EPSJDGetFix_JD);
606: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetConstCorrectionTol_C",EPSJDSetConstCorrectionTol_JD);
607: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetConstCorrectionTol_C",EPSJDGetConstCorrectionTol_JD);
608: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBOrth_C",EPSXDSetBOrth_XD);
609: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBOrth_C",EPSXDGetBOrth_XD);
610: PetscFunctionReturn(0);
611: }