1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepceps.h" I*/
32: #include <../src/eps/impls/davidson/davidson.h>
34: PetscErrorCode EPSSetFromOptions_JD(PetscOptionItems *PetscOptionsObject,EPS eps) 35: {
37: PetscBool flg,flg2,op,orth;
38: PetscInt opi,opi0;
39: PetscReal opf;
42: PetscOptionsHead(PetscOptionsObject,"EPS Jacobi-Davidson (JD) Options");
44: EPSJDGetKrylovStart(eps,&op);
45: PetscOptionsBool("-eps_jd_krylov_start","Start the search subspace with a Krylov basis","EPSJDSetKrylovStart",op,&op,&flg);
46: if (flg) { EPSJDSetKrylovStart(eps,op); }
48: EPSJDGetBOrth(eps,&orth);
49: PetscOptionsBool("-eps_jd_borth","Use B-orthogonalization in the search subspace","EPSJDSetBOrth",op,&op,&flg);
50: if (flg) { EPSJDSetBOrth(eps,op); }
52: EPSJDGetBlockSize(eps,&opi);
53: PetscOptionsInt("-eps_jd_blocksize","Number of vectors to add to the search subspace","EPSJDSetBlockSize",opi,&opi,&flg);
54: if (flg) { EPSJDSetBlockSize(eps,opi); }
56: EPSJDGetRestart(eps,&opi,&opi0);
57: PetscOptionsInt("-eps_jd_minv","Size of the search subspace after restarting","EPSJDSetRestart",opi,&opi,&flg);
58: PetscOptionsInt("-eps_jd_plusk","Number of eigenvectors saved from the previous iteration when restarting","EPSJDSetRestart",opi0,&opi0,&flg2);
59: if (flg || flg2) { EPSJDSetRestart(eps,opi,opi0); }
61: EPSJDGetInitialSize(eps,&opi);
62: PetscOptionsInt("-eps_jd_initial_size","Initial size of the search subspace","EPSJDSetInitialSize",opi,&opi,&flg);
63: if (flg) { EPSJDSetInitialSize(eps,opi); }
65: EPSJDGetFix(eps,&opf);
66: PetscOptionsReal("-eps_jd_fix","Tolerance for changing the target in the correction equation","EPSJDSetFix",opf,&opf,&flg);
67: if (flg) { EPSJDSetFix(eps,opf); }
69: EPSJDGetConstCorrectionTol(eps,&op);
70: PetscOptionsBool("-eps_jd_const_correction_tol","Disable the dynamic stopping criterion when solving the correction equation","EPSJDSetConstCorrectionTol",op,&op,&flg);
71: if (flg) { EPSJDSetConstCorrectionTol(eps,op); }
73: PetscOptionsTail();
74: return(0);
75: }
77: PetscErrorCode EPSSetDefaultST_JD(EPS eps) 78: {
80: KSP ksp;
83: if (!((PetscObject)eps->st)->type_name) {
84: STSetType(eps->st,STPRECOND);
85: STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
86: }
87: STGetKSP(eps->st,&ksp);
88: if (!((PetscObject)ksp)->type_name) {
89: KSPSetType(ksp,KSPBCGSL);
90: KSPSetTolerances(ksp,1e-4,PETSC_DEFAULT,PETSC_DEFAULT,90);
91: }
92: return(0);
93: }
95: PetscErrorCode EPSSetUp_JD(EPS eps) 96: {
98: PetscBool t;
99: KSP ksp;
102: /* Setup common for all davidson solvers */
103: EPSSetUp_XD(eps);
105: /* Check some constraints */
106: STGetKSP(eps->st,&ksp);
107: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&t);
108: if (t) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"EPSJD does not work with KSPPREONLY");
109: return(0);
110: }
112: PetscErrorCode EPSView_JD(EPS eps,PetscViewer viewer)113: {
115: PetscBool isascii,opb;
116: PetscReal opf;
117: PetscInt opi,opi0;
120: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
121: if (isascii) {
122: EPSXDGetBOrth_XD(eps,&opb);
123: if (opb) {
124: PetscViewerASCIIPrintf(viewer," search subspace is B-orthogonalized\n");
125: } else {
126: PetscViewerASCIIPrintf(viewer," search subspace is orthogonalized\n");
127: }
128: EPSXDGetBlockSize_XD(eps,&opi);
129: PetscViewerASCIIPrintf(viewer," block size=%D\n",opi);
130: EPSXDGetKrylovStart_XD(eps,&opb);
131: if (!opb) {
132: PetscViewerASCIIPrintf(viewer," type of the initial subspace: non-Krylov\n");
133: } else {
134: PetscViewerASCIIPrintf(viewer," type of the initial subspace: Krylov\n");
135: }
136: EPSXDGetRestart_XD(eps,&opi,&opi0);
137: PetscViewerASCIIPrintf(viewer," size of the subspace after restarting: %D\n",opi);
138: PetscViewerASCIIPrintf(viewer," number of vectors after restarting from the previous iteration: %D\n",opi0);
140: EPSJDGetFix_JD(eps,&opf);
141: PetscViewerASCIIPrintf(viewer," threshold for changing the target in the correction equation (fix): %g\n",(double)opf);
143: EPSJDGetConstCorrectionTol_JD(eps,&opb);
144: if (!opb) { PetscViewerASCIIPrintf(viewer," using dynamic tolerance for the correction equation\n"); }
145: }
146: return(0);
147: }
149: PetscErrorCode EPSDestroy_JD(EPS eps)150: {
151: PetscErrorCode ierr;
154: PetscFree(eps->data);
155: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetKrylovStart_C",NULL);
156: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetKrylovStart_C",NULL);
157: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBlockSize_C",NULL);
158: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBlockSize_C",NULL);
159: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetRestart_C",NULL);
160: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetRestart_C",NULL);
161: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetInitialSize_C",NULL);
162: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetInitialSize_C",NULL);
163: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetFix_C",NULL);
164: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetFix_C",NULL);
165: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetConstCorrectionTol_C",NULL);
166: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetConstCorrectionTol_C",NULL);
167: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBOrth_C",NULL);
168: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBOrth_C",NULL);
169: return(0);
170: }
172: /*@
173: EPSJDSetKrylovStart - Activates or deactivates starting the searching
174: subspace with a Krylov basis.
176: Logically Collective on EPS178: Input Parameters:
179: + eps - the eigenproblem solver context
180: - krylovstart - boolean flag
182: Options Database Key:
183: . -eps_jd_krylov_start - Activates starting the searching subspace with a
184: Krylov basis
186: Level: advanced
188: .seealso: EPSJDGetKrylovStart()
189: @*/
190: PetscErrorCode EPSJDSetKrylovStart(EPS eps,PetscBool krylovstart)191: {
197: PetscTryMethod(eps,"EPSJDSetKrylovStart_C",(EPS,PetscBool),(eps,krylovstart));
198: return(0);
199: }
201: /*@
202: EPSJDGetKrylovStart - Returns a flag indicating if the searching subspace is started with a
203: Krylov basis.
205: Not Collective
207: Input Parameter:
208: . eps - the eigenproblem solver context
210: Output Parameters:
211: . krylovstart - boolean flag indicating if the searching subspace is started
212: with a Krylov basis
214: Level: advanced
216: .seealso: EPSJDGetKrylovStart()
217: @*/
218: PetscErrorCode EPSJDGetKrylovStart(EPS eps,PetscBool *krylovstart)219: {
225: PetscUseMethod(eps,"EPSJDGetKrylovStart_C",(EPS,PetscBool*),(eps,krylovstart));
226: return(0);
227: }
229: /*@
230: EPSJDSetBlockSize - Sets the number of vectors to be added to the searching space
231: in every iteration.
233: Logically Collective on EPS235: Input Parameters:
236: + eps - the eigenproblem solver context
237: - blocksize - number of vectors added to the search space in every iteration
239: Options Database Key:
240: . -eps_jd_blocksize - number of vectors added to the searching space every iteration
242: Level: advanced
244: .seealso: EPSJDSetKrylovStart()
245: @*/
246: PetscErrorCode EPSJDSetBlockSize(EPS eps,PetscInt blocksize)247: {
253: PetscTryMethod(eps,"EPSJDSetBlockSize_C",(EPS,PetscInt),(eps,blocksize));
254: return(0);
255: }
257: /*@
258: EPSJDGetBlockSize - Returns the number of vectors to be added to the searching space
259: in every iteration.
261: Not Collective
263: Input Parameter:
264: . eps - the eigenproblem solver context
266: Output Parameter:
267: . blocksize - number of vectors added to the search space in every iteration
269: Level: advanced
271: .seealso: EPSJDSetBlockSize()
272: @*/
273: PetscErrorCode EPSJDGetBlockSize(EPS eps,PetscInt *blocksize)274: {
280: PetscUseMethod(eps,"EPSJDGetBlockSize_C",(EPS,PetscInt*),(eps,blocksize));
281: return(0);
282: }
284: /*@
285: EPSJDSetRestart - Sets the number of vectors of the searching space after
286: restarting and the number of vectors saved from the previous iteration.
288: Logically Collective on EPS290: Input Parameters:
291: + eps - the eigenproblem solver context
292: . minv - number of vectors of the searching subspace after restarting
293: - plusk - number of vectors saved from the previous iteration
295: Options Database Keys:
296: + -eps_jd_minv - number of vectors of the searching subspace after restarting
297: - -eps_jd_plusk - number of vectors saved from the previous iteration
299: Level: advanced
301: .seealso: EPSJDGetRestart()
302: @*/
303: PetscErrorCode EPSJDSetRestart(EPS eps,PetscInt minv,PetscInt plusk)304: {
311: PetscTryMethod(eps,"EPSJDSetRestart_C",(EPS,PetscInt,PetscInt),(eps,minv,plusk));
312: return(0);
313: }
315: /*@
316: EPSJDGetRestart - Gets the number of vectors of the searching space after
317: restarting and the number of vectors saved from the previous iteration.
319: Not Collective
321: Input Parameter:
322: . eps - the eigenproblem solver context
324: Output Parameter:
325: + minv - number of vectors of the searching subspace after restarting
326: - plusk - number of vectors saved from the previous iteration
328: Level: advanced
330: .seealso: EPSJDSetRestart()
331: @*/
332: PetscErrorCode EPSJDGetRestart(EPS eps,PetscInt *minv,PetscInt *plusk)333: {
338: PetscUseMethod(eps,"EPSJDGetRestart_C",(EPS,PetscInt*,PetscInt*),(eps,minv,plusk));
339: return(0);
340: }
342: /*@
343: EPSJDSetInitialSize - Sets the initial size of the searching space.
345: Logically Collective on EPS347: Input Parameters:
348: + eps - the eigenproblem solver context
349: - initialsize - number of vectors of the initial searching subspace
351: Options Database Key:
352: . -eps_jd_initial_size - number of vectors of the initial searching subspace
354: Notes:
355: If EPSJDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
356: EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
357: provided vectors are not enough, the solver completes the subspace with
358: random vectors. In the case of EPSJDGetKrylovStart() being PETSC_TRUE, the solver
359: gets the first vector provided by the user or, if not available, a random vector,
360: and expands the Krylov basis up to initialsize vectors.
362: Level: advanced
364: .seealso: EPSJDGetInitialSize(), EPSJDGetKrylovStart()
365: @*/
366: PetscErrorCode EPSJDSetInitialSize(EPS eps,PetscInt initialsize)367: {
373: PetscTryMethod(eps,"EPSJDSetInitialSize_C",(EPS,PetscInt),(eps,initialsize));
374: return(0);
375: }
377: /*@
378: EPSJDGetInitialSize - Returns the initial size of the searching space.
380: Not Collective
382: Input Parameter:
383: . eps - the eigenproblem solver context
385: Output Parameter:
386: . initialsize - number of vectors of the initial searching subspace
388: Notes:
389: If EPSJDGetKrylovStart() is PETSC_FALSE and the user provides vectors with
390: EPSSetInitialSpace(), up to initialsize vectors will be used; and if the
391: provided vectors are not enough, the solver completes the subspace with
392: random vectors. In the case of EPSJDGetKrylovStart() being PETSC_TRUE, the solver
393: gets the first vector provided by the user or, if not available, a random vector,
394: and expands the Krylov basis up to initialsize vectors.
396: Level: advanced
398: .seealso: EPSJDSetInitialSize(), EPSJDGetKrylovStart()
399: @*/
400: PetscErrorCode EPSJDGetInitialSize(EPS eps,PetscInt *initialsize)401: {
407: PetscUseMethod(eps,"EPSJDGetInitialSize_C",(EPS,PetscInt*),(eps,initialsize));
408: return(0);
409: }
411: PetscErrorCode EPSJDSetFix_JD(EPS eps,PetscReal fix)412: {
413: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
416: if (fix == PETSC_DEFAULT || fix == PETSC_DECIDE) fix = 0.01;
417: if (fix < 0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid fix value");
418: data->fix = fix;
419: return(0);
420: }
422: /*@
423: EPSJDSetFix - Sets the threshold for changing the target in the correction
424: equation.
426: Logically Collective on EPS428: Input Parameters:
429: + eps - the eigenproblem solver context
430: - fix - threshold for changing the target
432: Options Database Key:
433: . -eps_jd_fix - the fix value
435: Note:
436: The target in the correction equation is fixed at the first iterations.
437: When the norm of the residual vector is lower than the fix value,
438: the target is set to the corresponding eigenvalue.
440: Level: advanced
442: .seealso: EPSJDGetFix()
443: @*/
444: PetscErrorCode EPSJDSetFix(EPS eps,PetscReal fix)445: {
451: PetscTryMethod(eps,"EPSJDSetFix_C",(EPS,PetscReal),(eps,fix));
452: return(0);
453: }
455: PetscErrorCode EPSJDGetFix_JD(EPS eps,PetscReal *fix)456: {
457: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
460: *fix = data->fix;
461: return(0);
462: }
464: /*@
465: EPSJDGetFix - Returns the threshold for changing the target in the correction
466: equation.
468: Not Collective
470: Input Parameter:
471: . eps - the eigenproblem solver context
473: Output Parameter:
474: . fix - threshold for changing the target
476: Note:
477: The target in the correction equation is fixed at the first iterations.
478: When the norm of the residual vector is lower than the fix value,
479: the target is set to the corresponding eigenvalue.
481: Level: advanced
483: .seealso: EPSJDSetFix()
484: @*/
485: PetscErrorCode EPSJDGetFix(EPS eps,PetscReal *fix)486: {
492: PetscUseMethod(eps,"EPSJDGetFix_C",(EPS,PetscReal*),(eps,fix));
493: return(0);
494: }
496: PetscErrorCode EPSJDSetConstCorrectionTol_JD(EPS eps,PetscBool constant)497: {
498: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
501: data->dynamic = PetscNot(constant);
502: return(0);
503: }
505: /*@
506: EPSJDSetConstCorrectionTol - If true, deactivates the dynamic stopping criterion
507: (also called Newton) that sets the KSP relative tolerance
508: to 0.5**i, where i is the number of EPS iterations from the last converged value.
510: Logically Collective on EPS512: Input Parameters:
513: + eps - the eigenproblem solver context
514: - constant - if false, the KSP relative tolerance is set to 0.5**i.
516: Options Database Key:
517: . -eps_jd_const_correction_tol - Deactivates the dynamic stopping criterion.
519: Level: advanced
521: .seealso: EPSJDGetConstCorrectionTol()
522: @*/
523: PetscErrorCode EPSJDSetConstCorrectionTol(EPS eps,PetscBool constant)524: {
530: PetscTryMethod(eps,"EPSJDSetConstCorrectionTol_C",(EPS,PetscBool),(eps,constant));
531: return(0);
532: }
534: PetscErrorCode EPSJDGetConstCorrectionTol_JD(EPS eps,PetscBool *constant)535: {
536: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
539: *constant = PetscNot(data->dynamic);
540: return(0);
541: }
543: /*@
544: EPSJDGetConstCorrectionTol - Returns a flag indicating if the dynamic stopping is being used for
545: solving the correction equation. If the flag is false the KSP relative tolerance is set
546: to 0.5**i, where i is the number of EPS iterations from the last converged value.
548: Not Collective
550: Input Parameter:
551: . eps - the eigenproblem solver context
553: Output Parameters:
554: . constant - boolean flag indicating if the dynamic stopping criterion is not being used.
556: Level: advanced
558: .seealso: EPSJDGetConstCorrectionTol()
559: @*/
560: PetscErrorCode EPSJDGetConstCorrectionTol(EPS eps,PetscBool *constant)561: {
567: PetscUseMethod(eps,"EPSJDGetConstCorrectionTol_C",(EPS,PetscBool*),(eps,constant));
568: return(0);
569: }
571: /*@
572: EPSJDSetBOrth - Selects the orthogonalization that will be used in the search
573: subspace in case of generalized Hermitian problems.
575: Logically Collective on EPS577: Input Parameters:
578: + eps - the eigenproblem solver context
579: - borth - whether to B-orthogonalize the search subspace
581: Options Database Key:
582: . -eps_jd_borth - Set the orthogonalization used in the search subspace
584: Level: advanced
586: .seealso: EPSJDGetBOrth()
587: @*/
588: PetscErrorCode EPSJDSetBOrth(EPS eps,PetscBool borth)589: {
595: PetscTryMethod(eps,"EPSJDSetBOrth_C",(EPS,PetscBool),(eps,borth));
596: return(0);
597: }
599: /*@
600: EPSJDGetBOrth - Returns the orthogonalization used in the search
601: subspace in case of generalized Hermitian problems.
603: Not Collective
605: Input Parameter:
606: . eps - the eigenproblem solver context
608: Output Parameters:
609: . borth - whether to B-orthogonalize the search subspace
611: Level: advanced
613: .seealso: EPSJDSetBOrth()
614: @*/
615: PetscErrorCode EPSJDGetBOrth(EPS eps,PetscBool *borth)616: {
622: PetscUseMethod(eps,"EPSJDGetBOrth_C",(EPS,PetscBool*),(eps,borth));
623: return(0);
624: }
626: SLEPC_EXTERN PetscErrorCode EPSCreate_JD(EPS eps)627: {
629: EPS_DAVIDSON *data;
632: PetscNewLog(eps,&data);
633: eps->data = (void*)data;
635: data->blocksize = 1;
636: data->initialsize = 6;
637: data->minv = 6;
638: data->plusk = PETSC_DEFAULT;
639: data->ipB = PETSC_TRUE;
640: data->fix = 0.01;
641: data->krylovstart = PETSC_FALSE;
642: data->dynamic = PETSC_FALSE;
644: eps->useds = PETSC_TRUE;
645: eps->categ = EPS_CATEGORY_PRECOND;
647: eps->ops->solve = EPSSolve_XD;
648: eps->ops->setup = EPSSetUp_JD;
649: eps->ops->setfromoptions = EPSSetFromOptions_JD;
650: eps->ops->destroy = EPSDestroy_JD;
651: eps->ops->reset = EPSReset_XD;
652: eps->ops->view = EPSView_JD;
653: eps->ops->backtransform = EPSBackTransform_Default;
654: eps->ops->computevectors = EPSComputeVectors_XD;
655: eps->ops->setdefaultst = EPSSetDefaultST_JD;
657: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetKrylovStart_C",EPSXDSetKrylovStart_XD);
658: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetKrylovStart_C",EPSXDGetKrylovStart_XD);
659: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBlockSize_C",EPSXDSetBlockSize_XD);
660: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBlockSize_C",EPSXDGetBlockSize_XD);
661: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetRestart_C",EPSXDSetRestart_XD);
662: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetRestart_C",EPSXDGetRestart_XD);
663: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetInitialSize_C",EPSXDSetInitialSize_XD);
664: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetInitialSize_C",EPSXDGetInitialSize_XD);
665: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetFix_C",EPSJDSetFix_JD);
666: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetFix_C",EPSJDGetFix_JD);
667: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetConstCorrectionTol_C",EPSJDSetConstCorrectionTol_JD);
668: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetConstCorrectionTol_C",EPSJDGetConstCorrectionTol_JD);
669: PetscObjectComposeFunction((PetscObject)eps,"EPSJDSetBOrth_C",EPSXDSetBOrth_XD);
670: PetscObjectComposeFunction((PetscObject)eps,"EPSJDGetBOrth_C",EPSXDGetBOrth_XD);
671: return(0);
672: }