Actual source code: krylovschur.c
slepc-3.17.2 2022-08-09
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: "krylovschur"
13: Method: Krylov-Schur
15: Algorithm:
17: Single-vector Krylov-Schur method for non-symmetric problems,
18: including harmonic extraction.
20: References:
22: [1] "Krylov-Schur Methods in SLEPc", SLEPc Technical Report STR-7,
23: available at https://slepc.upv.es.
25: [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",
26: SIAM J. Matrix Anal. App. 23(3):601-614, 2001.
28: [3] "Practical Implementation of Harmonic Krylov-Schur", SLEPc Technical
29: Report STR-9, available at https://slepc.upv.es.
30: */
32: #include <slepc/private/epsimpl.h>
33: #include "krylovschur.h"
35: PetscErrorCode EPSGetArbitraryValues(EPS eps,PetscScalar *rr,PetscScalar *ri)
36: {
37: PetscInt i,newi,ld,n,l;
38: Vec xr=eps->work[0],xi=eps->work[1];
39: PetscScalar re,im,*Zr,*Zi,*X;
41: DSGetLeadingDimension(eps->ds,&ld);
42: DSGetDimensions(eps->ds,&n,&l,NULL,NULL);
43: for (i=l;i<n;i++) {
44: re = eps->eigr[i];
45: im = eps->eigi[i];
46: STBackTransform(eps->st,1,&re,&im);
47: newi = i;
48: DSVectors(eps->ds,DS_MAT_X,&newi,NULL);
49: DSGetArray(eps->ds,DS_MAT_X,&X);
50: Zr = X+i*ld;
51: if (newi==i+1) Zi = X+newi*ld;
52: else Zi = NULL;
53: EPSComputeRitzVector(eps,Zr,Zi,eps->V,xr,xi);
54: DSRestoreArray(eps->ds,DS_MAT_X,&X);
55: (*eps->arbitrary)(re,im,xr,xi,rr+i,ri+i,eps->arbitraryctx);
56: }
57: PetscFunctionReturn(0);
58: }
60: static PetscErrorCode EPSSetUp_KrylovSchur_Filter(EPS eps)
61: {
62: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
63: PetscBool estimaterange=PETSC_TRUE;
64: PetscReal rleft,rright;
65: Mat A;
67: EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
68: EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
70: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
71: if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL*1e-2; /* use tighter tolerance */
72: STFilterSetInterval(eps->st,eps->inta,eps->intb);
73: if (!ctx->estimatedrange) {
74: STFilterGetRange(eps->st,&rleft,&rright);
75: estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
76: }
77: if (estimaterange) { /* user did not set a range */
78: STGetMatrix(eps->st,0,&A);
79: MatEstimateSpectralRange_EPS(A,&rleft,&rright);
80: PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright);
81: STFilterSetRange(eps->st,rleft,rright);
82: ctx->estimatedrange = PETSC_TRUE;
83: }
84: if (eps->ncv==PETSC_DEFAULT && eps->nev==1) eps->nev = 40; /* user did not provide nev estimation */
85: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
87: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
88: PetscFunctionReturn(0);
89: }
91: PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
92: {
93: PetscReal eta;
94: PetscBool isfilt=PETSC_FALSE;
95: BVOrthogType otype;
96: BVOrthogBlockType obtype;
97: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
98: enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;
100: if (eps->which==EPS_ALL) { /* default values in case of spectrum slicing or polynomial filter */
101: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
102: if (isfilt) EPSSetUp_KrylovSchur_Filter(eps);
103: else EPSSetUp_KrylovSchur_Slice(eps);
104: } else {
105: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
107: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
108: if (!eps->which) EPSSetWhichEigenpairs_Default(eps);
109: }
112: EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
116: if (!ctx->keep) ctx->keep = 0.5;
118: EPSAllocateSolution(eps,1);
119: EPS_SetInnerProduct(eps);
120: if (eps->arbitrary) EPSSetWorkVecs(eps,2);
121: else if (eps->ishermitian && !eps->ispositive) EPSSetWorkVecs(eps,1);
123: /* dispatch solve method */
124: if (eps->ishermitian) {
125: if (eps->which==EPS_ALL) {
126: EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
127: variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
128: } else if (eps->isgeneralized && !eps->ispositive) {
129: variant = EPS_KS_INDEF;
130: } else {
131: switch (eps->extraction) {
132: case EPS_RITZ: variant = EPS_KS_SYMM; break;
133: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
134: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
135: }
136: }
137: } else if (eps->twosided) {
138: variant = EPS_KS_TWOSIDED;
139: } else {
140: switch (eps->extraction) {
141: case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
142: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
143: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
144: }
145: }
146: switch (variant) {
147: case EPS_KS_DEFAULT:
148: eps->ops->solve = EPSSolve_KrylovSchur_Default;
149: eps->ops->computevectors = EPSComputeVectors_Schur;
150: DSSetType(eps->ds,DSNHEP);
151: DSSetExtraRow(eps->ds,PETSC_TRUE);
152: DSAllocate(eps->ds,eps->ncv+1);
153: break;
154: case EPS_KS_SYMM:
155: case EPS_KS_FILTER:
156: eps->ops->solve = EPSSolve_KrylovSchur_Default;
157: eps->ops->computevectors = EPSComputeVectors_Hermitian;
158: DSSetType(eps->ds,DSHEP);
159: DSSetCompact(eps->ds,PETSC_TRUE);
160: DSSetExtraRow(eps->ds,PETSC_TRUE);
161: DSAllocate(eps->ds,eps->ncv+1);
162: break;
163: case EPS_KS_SLICE:
164: eps->ops->solve = EPSSolve_KrylovSchur_Slice;
165: eps->ops->computevectors = EPSComputeVectors_Slice;
166: break;
167: case EPS_KS_INDEF:
168: eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
169: eps->ops->computevectors = EPSComputeVectors_Indefinite;
170: DSSetType(eps->ds,DSGHIEP);
171: DSSetCompact(eps->ds,PETSC_TRUE);
172: DSSetExtraRow(eps->ds,PETSC_TRUE);
173: DSAllocate(eps->ds,eps->ncv+1);
174: /* force reorthogonalization for pseudo-Lanczos */
175: BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype);
176: BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
177: break;
178: case EPS_KS_TWOSIDED:
179: eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
180: eps->ops->computevectors = EPSComputeVectors_Schur;
181: DSSetType(eps->ds,DSNHEPTS);
182: DSAllocate(eps->ds,eps->ncv+1);
183: DSSetExtraRow(eps->ds,PETSC_TRUE);
184: break;
185: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
186: }
187: PetscFunctionReturn(0);
188: }
190: PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
191: {
192: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
193: SlepcSC sc;
194: PetscBool isfilt;
196: EPSSetUpSort_Default(eps);
197: if (eps->which==EPS_ALL) {
198: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
199: if (isfilt) {
200: DSGetSlepcSC(eps->ds,&sc);
201: sc->rg = NULL;
202: sc->comparison = SlepcCompareLargestReal;
203: sc->comparisonctx = NULL;
204: sc->map = NULL;
205: sc->mapobj = NULL;
206: } else {
207: if (!ctx->global && ctx->sr->numEigs>0) {
208: DSGetSlepcSC(eps->ds,&sc);
209: sc->rg = NULL;
210: sc->comparison = SlepcCompareLargestMagnitude;
211: sc->comparisonctx = NULL;
212: sc->map = NULL;
213: sc->mapobj = NULL;
214: }
215: }
216: }
217: PetscFunctionReturn(0);
218: }
220: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
221: {
222: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
223: PetscInt j,*pj,k,l,nv,ld,nconv;
224: Mat U,Op,H;
225: PetscScalar *g;
226: PetscReal beta,gamma=1.0,*a,*b;
227: PetscBool breakdown,harmonic,hermitian;
229: DSGetLeadingDimension(eps->ds,&ld);
230: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
231: hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
232: if (harmonic) PetscMalloc1(ld,&g);
233: if (eps->arbitrary) pj = &j;
234: else pj = NULL;
236: /* Get the starting Arnoldi vector */
237: EPSGetStartVector(eps,0,NULL);
238: l = 0;
240: /* Restart loop */
241: while (eps->reason == EPS_CONVERGED_ITERATING) {
242: eps->its++;
244: /* Compute an nv-step Arnoldi factorization */
245: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
246: DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
247: STGetOperator(eps->st,&Op);
248: if (hermitian) {
249: DSGetArrayReal(eps->ds,DS_MAT_T,&a);
250: b = a + ld;
251: BVMatLanczos(eps->V,Op,a,b,eps->nconv+l,&nv,&breakdown);
252: beta = b[nv-1];
253: DSRestoreArrayReal(eps->ds,DS_MAT_T,&a);
254: } else {
255: DSGetMat(eps->ds,DS_MAT_A,&H);
256: BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown);
257: DSRestoreMat(eps->ds,DS_MAT_A,&H);
258: }
259: STRestoreOperator(eps->st,&Op);
260: DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l);
261: DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE);
262: BVSetActiveColumns(eps->V,eps->nconv,nv);
264: /* Compute translation of Krylov decomposition if harmonic extraction used */
265: if (PetscUnlikely(harmonic)) DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma);
267: /* Solve projected problem */
268: DSSolve(eps->ds,eps->eigr,eps->eigi);
269: if (PetscUnlikely(eps->arbitrary)) {
270: EPSGetArbitraryValues(eps,eps->rr,eps->ri);
271: j=1;
272: }
273: DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj);
274: DSUpdateExtraRow(eps->ds);
275: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
277: /* Check convergence */
278: EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k);
279: (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
280: nconv = k;
282: /* Update l */
283: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
284: else {
285: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
286: if (!hermitian) DSGetTruncateSize(eps->ds,k,nv,&l);
287: }
288: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
289: if (l) PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l);
291: if (eps->reason == EPS_CONVERGED_ITERATING) {
292: if (PetscUnlikely(breakdown || k==nv)) {
293: /* Start a new Arnoldi factorization */
294: PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta);
295: if (k<eps->nev) {
296: EPSGetStartVector(eps,k,&breakdown);
297: if (breakdown) {
298: eps->reason = EPS_DIVERGED_BREAKDOWN;
299: PetscInfo(eps,"Unable to generate more start vectors\n");
300: }
301: }
302: } else {
303: /* Undo translation of Krylov decomposition */
304: if (PetscUnlikely(harmonic)) {
305: DSSetDimensions(eps->ds,nv,k,l);
306: DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma);
307: /* gamma u^ = u - U*g~ */
308: BVSetActiveColumns(eps->V,0,nv);
309: BVMultColumn(eps->V,-1.0,1.0,nv,g);
310: BVScaleColumn(eps->V,nv,1.0/gamma);
311: BVSetActiveColumns(eps->V,eps->nconv,nv);
312: DSSetDimensions(eps->ds,nv,k,nv);
313: }
314: /* Prepare the Rayleigh quotient for restart */
315: DSTruncate(eps->ds,k+l,PETSC_FALSE);
316: }
317: }
318: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
319: DSGetMat(eps->ds,DS_MAT_Q,&U);
320: BVMultInPlace(eps->V,U,eps->nconv,k+l);
321: MatDestroy(&U);
323: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) BVCopyColumn(eps->V,nv,k+l);
324: eps->nconv = k;
325: EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
326: }
328: if (harmonic) PetscFree(g);
329: DSTruncate(eps->ds,eps->nconv,PETSC_TRUE);
330: PetscFunctionReturn(0);
331: }
333: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
334: {
335: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
337: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
338: else {
340: ctx->keep = keep;
341: }
342: PetscFunctionReturn(0);
343: }
345: /*@
346: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
347: method, in particular the proportion of basis vectors that must be kept
348: after restart.
350: Logically Collective on eps
352: Input Parameters:
353: + eps - the eigenproblem solver context
354: - keep - the number of vectors to be kept at restart
356: Options Database Key:
357: . -eps_krylovschur_restart - Sets the restart parameter
359: Notes:
360: Allowed values are in the range [0.1,0.9]. The default is 0.5.
362: Level: advanced
364: .seealso: EPSKrylovSchurGetRestart()
365: @*/
366: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
367: {
370: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
371: PetscFunctionReturn(0);
372: }
374: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
375: {
376: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
378: *keep = ctx->keep;
379: PetscFunctionReturn(0);
380: }
382: /*@
383: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
384: Krylov-Schur method.
386: Not Collective
388: Input Parameter:
389: . eps - the eigenproblem solver context
391: Output Parameter:
392: . keep - the restart parameter
394: Level: advanced
396: .seealso: EPSKrylovSchurSetRestart()
397: @*/
398: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
399: {
402: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
403: PetscFunctionReturn(0);
404: }
406: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
407: {
408: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
410: ctx->lock = lock;
411: PetscFunctionReturn(0);
412: }
414: /*@
415: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
416: the Krylov-Schur method.
418: Logically Collective on eps
420: Input Parameters:
421: + eps - the eigenproblem solver context
422: - lock - true if the locking variant must be selected
424: Options Database Key:
425: . -eps_krylovschur_locking - Sets the locking flag
427: Notes:
428: The default is to lock converged eigenpairs when the method restarts.
429: This behaviour can be changed so that all directions are kept in the
430: working subspace even if already converged to working accuracy (the
431: non-locking variant).
433: Level: advanced
435: .seealso: EPSKrylovSchurGetLocking()
436: @*/
437: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
438: {
441: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
442: PetscFunctionReturn(0);
443: }
445: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
446: {
447: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
449: *lock = ctx->lock;
450: PetscFunctionReturn(0);
451: }
453: /*@
454: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
455: method.
457: Not Collective
459: Input Parameter:
460: . eps - the eigenproblem solver context
462: Output Parameter:
463: . lock - the locking flag
465: Level: advanced
467: .seealso: EPSKrylovSchurSetLocking()
468: @*/
469: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
470: {
473: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
474: PetscFunctionReturn(0);
475: }
477: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
478: {
479: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
480: PetscMPIInt size;
481: PetscInt newnpart;
483: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
484: newnpart = 1;
485: } else {
486: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
488: newnpart = npart;
489: }
490: if (ctx->npart!=newnpart) {
491: if (ctx->npart>1) {
492: PetscSubcommDestroy(&ctx->subc);
493: if (ctx->commset) {
494: MPI_Comm_free(&ctx->commrank);
495: ctx->commset = PETSC_FALSE;
496: }
497: }
498: EPSDestroy(&ctx->eps);
499: ctx->npart = newnpart;
500: eps->state = EPS_STATE_INITIAL;
501: }
502: PetscFunctionReturn(0);
503: }
505: /*@
506: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
507: case of doing spectrum slicing for a computational interval with the
508: communicator split in several sub-communicators.
510: Logically Collective on eps
512: Input Parameters:
513: + eps - the eigenproblem solver context
514: - npart - number of partitions
516: Options Database Key:
517: . -eps_krylovschur_partitions <npart> - Sets the number of partitions
519: Notes:
520: By default, npart=1 so all processes in the communicator participate in
521: the processing of the whole interval. If npart>1 then the interval is
522: divided into npart subintervals, each of them being processed by a
523: subset of processes.
525: The interval is split proportionally unless the separation points are
526: specified with EPSKrylovSchurSetSubintervals().
528: Level: advanced
530: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
531: @*/
532: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
533: {
536: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
537: PetscFunctionReturn(0);
538: }
540: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
541: {
542: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
544: *npart = ctx->npart;
545: PetscFunctionReturn(0);
546: }
548: /*@
549: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
550: communicator in case of spectrum slicing.
552: Not Collective
554: Input Parameter:
555: . eps - the eigenproblem solver context
557: Output Parameter:
558: . npart - number of partitions
560: Level: advanced
562: .seealso: EPSKrylovSchurSetPartitions()
563: @*/
564: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
565: {
568: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
569: PetscFunctionReturn(0);
570: }
572: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
573: {
574: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
576: ctx->detect = detect;
577: eps->state = EPS_STATE_INITIAL;
578: PetscFunctionReturn(0);
579: }
581: /*@
582: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
583: zeros during the factorizations throughout the spectrum slicing computation.
585: Logically Collective on eps
587: Input Parameters:
588: + eps - the eigenproblem solver context
589: - detect - check for zeros
591: Options Database Key:
592: . -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
593: bool value (0/1/no/yes/true/false)
595: Notes:
596: A zero in the factorization indicates that a shift coincides with an eigenvalue.
598: This flag is turned off by default, and may be necessary in some cases,
599: especially when several partitions are being used. This feature currently
600: requires an external package for factorizations with support for zero
601: detection, e.g. MUMPS.
603: Level: advanced
605: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
606: @*/
607: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
608: {
611: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
612: PetscFunctionReturn(0);
613: }
615: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
616: {
617: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
619: *detect = ctx->detect;
620: PetscFunctionReturn(0);
621: }
623: /*@
624: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
625: in spectrum slicing.
627: Not Collective
629: Input Parameter:
630: . eps - the eigenproblem solver context
632: Output Parameter:
633: . detect - whether zeros detection is enforced during factorizations
635: Level: advanced
637: .seealso: EPSKrylovSchurSetDetectZeros()
638: @*/
639: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
640: {
643: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
644: PetscFunctionReturn(0);
645: }
647: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
648: {
649: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
652: ctx->nev = nev;
653: if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
654: ctx->ncv = PETSC_DEFAULT;
655: } else {
657: ctx->ncv = ncv;
658: }
659: if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
660: ctx->mpd = PETSC_DEFAULT;
661: } else {
663: ctx->mpd = mpd;
664: }
665: eps->state = EPS_STATE_INITIAL;
666: PetscFunctionReturn(0);
667: }
669: /*@
670: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
671: step in case of doing spectrum slicing for a computational interval.
672: The meaning of the parameters is the same as in EPSSetDimensions().
674: Logically Collective on eps
676: Input Parameters:
677: + eps - the eigenproblem solver context
678: . nev - number of eigenvalues to compute
679: . ncv - the maximum dimension of the subspace to be used by the subsolve
680: - mpd - the maximum dimension allowed for the projected problem
682: Options Database Key:
683: + -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
684: . -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
685: - -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension
687: Level: advanced
689: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
690: @*/
691: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
692: {
697: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
698: PetscFunctionReturn(0);
699: }
701: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
702: {
703: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
705: if (nev) *nev = ctx->nev;
706: if (ncv) *ncv = ctx->ncv;
707: if (mpd) *mpd = ctx->mpd;
708: PetscFunctionReturn(0);
709: }
711: /*@
712: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
713: step in case of doing spectrum slicing for a computational interval.
715: Not Collective
717: Input Parameter:
718: . eps - the eigenproblem solver context
720: Output Parameters:
721: + nev - number of eigenvalues to compute
722: . ncv - the maximum dimension of the subspace to be used by the subsolve
723: - mpd - the maximum dimension allowed for the projected problem
725: Level: advanced
727: .seealso: EPSKrylovSchurSetDimensions()
728: @*/
729: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
730: {
732: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
733: PetscFunctionReturn(0);
734: }
736: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
737: {
738: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
739: PetscInt i;
743: if (ctx->subintervals) PetscFree(ctx->subintervals);
744: PetscMalloc1(ctx->npart+1,&ctx->subintervals);
745: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
746: ctx->subintset = PETSC_TRUE;
747: eps->state = EPS_STATE_INITIAL;
748: PetscFunctionReturn(0);
749: }
751: /*@C
752: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
753: subintervals to be used in spectrum slicing with several partitions.
755: Logically Collective on eps
757: Input Parameters:
758: + eps - the eigenproblem solver context
759: - subint - array of real values specifying subintervals
761: Notes:
762: This function must be called after EPSKrylovSchurSetPartitions(). For npart
763: partitions, the argument subint must contain npart+1 real values sorted in
764: ascending order, subint_0, subint_1, ..., subint_npart, where the first
765: and last values must coincide with the interval endpoints set with
766: EPSSetInterval().
768: The subintervals are then defined by two consecutive points [subint_0,subint_1],
769: [subint_1,subint_2], and so on.
771: Level: advanced
773: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
774: @*/
775: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint)
776: {
778: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
779: PetscFunctionReturn(0);
780: }
782: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
783: {
784: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
785: PetscInt i;
787: if (!ctx->subintset) {
790: }
791: PetscMalloc1(ctx->npart+1,subint);
792: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
793: PetscFunctionReturn(0);
794: }
796: /*@C
797: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
798: subintervals used in spectrum slicing with several partitions.
800: Logically Collective on eps
802: Input Parameter:
803: . eps - the eigenproblem solver context
805: Output Parameter:
806: . subint - array of real values specifying subintervals
808: Notes:
809: If the user passed values with EPSKrylovSchurSetSubintervals(), then the
810: same values are returned. Otherwise, the values computed internally are
811: obtained.
813: This function is only available for spectrum slicing runs.
815: The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
816: and should be freed by the user.
818: Fortran Notes:
819: The calling sequence from Fortran is
820: .vb
821: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
822: double precision subint(npart+1) output
823: .ve
825: Level: advanced
827: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
828: @*/
829: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
830: {
833: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
834: PetscFunctionReturn(0);
835: }
837: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
838: {
839: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
840: PetscInt i,numsh;
841: EPS_SR sr = ctx->sr;
845: switch (eps->state) {
846: case EPS_STATE_INITIAL:
847: break;
848: case EPS_STATE_SETUP:
849: numsh = ctx->npart+1;
850: if (n) *n = numsh;
851: if (shifts) {
852: PetscMalloc1(numsh,shifts);
853: (*shifts)[0] = eps->inta;
854: if (ctx->npart==1) (*shifts)[1] = eps->intb;
855: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
856: }
857: if (inertias) {
858: PetscMalloc1(numsh,inertias);
859: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
860: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
861: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
862: }
863: break;
864: case EPS_STATE_SOLVED:
865: case EPS_STATE_EIGENVECTORS:
866: numsh = ctx->nshifts;
867: if (n) *n = numsh;
868: if (shifts) {
869: PetscMalloc1(numsh,shifts);
870: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
871: }
872: if (inertias) {
873: PetscMalloc1(numsh,inertias);
874: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
875: }
876: break;
877: }
878: PetscFunctionReturn(0);
879: }
881: /*@C
882: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
883: corresponding inertias in case of doing spectrum slicing for a
884: computational interval.
886: Not Collective
888: Input Parameter:
889: . eps - the eigenproblem solver context
891: Output Parameters:
892: + n - number of shifts, including the endpoints of the interval
893: . shifts - the values of the shifts used internally in the solver
894: - inertias - the values of the inertia in each shift
896: Notes:
897: If called after EPSSolve(), all shifts used internally by the solver are
898: returned (including both endpoints and any intermediate ones). If called
899: before EPSSolve() and after EPSSetUp() then only the information of the
900: endpoints of subintervals is available.
902: This function is only available for spectrum slicing runs.
904: The returned arrays should be freed by the user. Can pass NULL in any of
905: the two arrays if not required.
907: Fortran Notes:
908: The calling sequence from Fortran is
909: .vb
910: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
911: integer n
912: double precision shifts(*)
913: integer inertias(*)
914: .ve
915: The arrays should be at least of length n. The value of n can be determined
916: by an initial call
917: .vb
918: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
919: .ve
921: Level: advanced
923: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
924: @*/
925: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
926: {
929: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
930: PetscFunctionReturn(0);
931: }
933: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
934: {
935: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
936: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
940: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
941: if (n) *n = sr->numEigs;
942: if (v) BVCreateVec(sr->V,v);
943: PetscFunctionReturn(0);
944: }
946: /*@C
947: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
948: doing spectrum slicing for a computational interval with multiple
949: communicators.
951: Collective on the subcommunicator (if v is given)
953: Input Parameter:
954: . eps - the eigenproblem solver context
956: Output Parameters:
957: + k - index of the subinterval for the calling process
958: . n - number of eigenvalues found in the k-th subinterval
959: - v - a vector owned by processes in the subcommunicator with dimensions
960: compatible for locally computed eigenvectors (or NULL)
962: Notes:
963: This function is only available for spectrum slicing runs.
965: The returned Vec should be destroyed by the user.
967: Level: advanced
969: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
970: @*/
971: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
972: {
974: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
975: PetscFunctionReturn(0);
976: }
978: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
979: {
980: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
981: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
983: EPSCheckSolved(eps,1);
986: if (eig) *eig = sr->eigr[sr->perm[i]];
987: if (v) BVCopyVec(sr->V,sr->perm[i],v);
988: PetscFunctionReturn(0);
989: }
991: /*@
992: EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
993: internally in the subcommunicator to which the calling process belongs.
995: Collective on the subcommunicator (if v is given)
997: Input Parameters:
998: + eps - the eigenproblem solver context
999: - i - index of the solution
1001: Output Parameters:
1002: + eig - the eigenvalue
1003: - v - the eigenvector
1005: Notes:
1006: It is allowed to pass NULL for v if the eigenvector is not required.
1007: Otherwise, the caller must provide a valid Vec objects, i.e.,
1008: it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().
1010: The index i should be a value between 0 and n-1, where n is the number of
1011: vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().
1013: Level: advanced
1015: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1016: @*/
1017: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1018: {
1021: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1022: PetscFunctionReturn(0);
1023: }
1025: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1026: {
1027: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1031: EPSGetOperators(ctx->eps,A,B);
1032: PetscFunctionReturn(0);
1033: }
1035: /*@C
1036: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1037: internally in the subcommunicator to which the calling process belongs.
1039: Collective on the subcommunicator
1041: Input Parameter:
1042: . eps - the eigenproblem solver context
1044: Output Parameters:
1045: + A - the matrix associated with the eigensystem
1046: - B - the second matrix in the case of generalized eigenproblems
1048: Notes:
1049: This is the analog of EPSGetOperators(), but returns the matrices distributed
1050: differently (in the subcommunicator rather than in the parent communicator).
1052: These matrices should not be modified by the user.
1054: Level: advanced
1056: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1057: @*/
1058: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1059: {
1061: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1062: PetscFunctionReturn(0);
1063: }
1065: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1066: {
1067: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1068: Mat A,B=NULL,Ag,Bg=NULL;
1069: PetscBool reuse=PETSC_TRUE;
1073: EPSGetOperators(eps,&Ag,&Bg);
1074: EPSGetOperators(ctx->eps,&A,&B);
1076: MatScale(A,a);
1077: if (Au) MatAXPY(A,ap,Au,str);
1078: if (B) MatScale(B,b);
1079: if (Bu) MatAXPY(B,bp,Bu,str);
1080: EPSSetOperators(ctx->eps,A,B);
1082: /* Update stored matrix state */
1083: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1084: PetscObjectStateGet((PetscObject)A,&subctx->Astate);
1085: if (B) PetscObjectStateGet((PetscObject)B,&subctx->Bstate);
1087: /* Update matrices in the parent communicator if requested by user */
1088: if (globalup) {
1089: if (ctx->npart>1) {
1090: if (!ctx->isrow) {
1091: MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol);
1092: reuse = PETSC_FALSE;
1093: }
1094: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1095: if (ctx->submata && !reuse) MatDestroyMatrices(1,&ctx->submata);
1096: MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata);
1097: MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag);
1098: if (B) {
1099: if (ctx->submatb && !reuse) MatDestroyMatrices(1,&ctx->submatb);
1100: MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb);
1101: MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg);
1102: }
1103: }
1104: PetscObjectStateGet((PetscObject)Ag,&ctx->Astate);
1105: if (Bg) PetscObjectStateGet((PetscObject)Bg,&ctx->Bstate);
1106: }
1107: EPSSetOperators(eps,Ag,Bg);
1108: PetscFunctionReturn(0);
1109: }
1111: /*@
1112: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1113: internally in the subcommunicator to which the calling process belongs.
1115: Collective on eps
1117: Input Parameters:
1118: + eps - the eigenproblem solver context
1119: . s - scalar that multiplies the existing A matrix
1120: . a - scalar used in the axpy operation on A
1121: . Au - matrix used in the axpy operation on A
1122: . t - scalar that multiplies the existing B matrix
1123: . b - scalar used in the axpy operation on B
1124: . Bu - matrix used in the axpy operation on B
1125: . str - structure flag
1126: - globalup - flag indicating if global matrices must be updated
1128: Notes:
1129: This function modifies the eigenproblem matrices at the subcommunicator level,
1130: and optionally updates the global matrices in the parent communicator. The updates
1131: are expressed as A <-- s*A + a*Au, B <-- t*B + b*Bu.
1133: It is possible to update one of the matrices, or both.
1135: The matrices Au and Bu must be equal in all subcommunicators.
1137: The str flag is passed to the MatAXPY() operations to perform the updates.
1139: If globalup is true, communication is carried out to reconstruct the updated
1140: matrices in the parent communicator. The user must be warned that if global
1141: matrices are not in sync with subcommunicator matrices, the errors computed
1142: by EPSComputeError() will be wrong even if the computed solution is correct
1143: (the synchronization may be done only once at the end).
1145: Level: advanced
1147: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1148: @*/
1149: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1150: {
1160: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1161: PetscFunctionReturn(0);
1162: }
1164: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1165: {
1166: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1167: Mat A,B=NULL,Ar=NULL,Br=NULL;
1168: PetscMPIInt rank;
1169: PetscObjectState Astate,Bstate=0;
1170: PetscObjectId Aid,Bid=0;
1171: STType sttype;
1172: PetscInt nmat;
1173: const char *prefix;
1174: MPI_Comm child;
1176: EPSGetOperators(eps,&A,&B);
1177: if (ctx->npart==1) {
1178: if (!ctx->eps) EPSCreate(((PetscObject)eps)->comm,&ctx->eps);
1179: EPSGetOptionsPrefix(eps,&prefix);
1180: EPSSetOptionsPrefix(ctx->eps,prefix);
1181: EPSSetOperators(ctx->eps,A,B);
1182: } else {
1183: PetscObjectStateGet((PetscObject)A,&Astate);
1184: PetscObjectGetId((PetscObject)A,&Aid);
1185: if (B) {
1186: PetscObjectStateGet((PetscObject)B,&Bstate);
1187: PetscObjectGetId((PetscObject)B,&Bid);
1188: }
1189: if (!ctx->subc) {
1190: /* Create context for subcommunicators */
1191: PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc);
1192: PetscSubcommSetNumber(ctx->subc,ctx->npart);
1193: PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS);
1194: PetscLogObjectMemory((PetscObject)eps,sizeof(PetscSubcomm));
1195: PetscSubcommGetChild(ctx->subc,&child);
1197: /* Duplicate matrices */
1198: MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar);
1199: PetscLogObjectParent((PetscObject)eps,(PetscObject)Ar);
1200: ctx->Astate = Astate;
1201: ctx->Aid = Aid;
1202: MatPropagateSymmetryOptions(A,Ar);
1203: if (B) {
1204: MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br);
1205: PetscLogObjectParent((PetscObject)eps,(PetscObject)Br);
1206: ctx->Bstate = Bstate;
1207: ctx->Bid = Bid;
1208: MatPropagateSymmetryOptions(B,Br);
1209: }
1210: } else {
1211: PetscSubcommGetChild(ctx->subc,&child);
1212: if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1213: STGetNumMatrices(ctx->eps->st,&nmat);
1214: if (nmat) EPSGetOperators(ctx->eps,&Ar,&Br);
1215: MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar);
1216: ctx->Astate = Astate;
1217: ctx->Aid = Aid;
1218: MatPropagateSymmetryOptions(A,Ar);
1219: if (B) {
1220: MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br);
1221: ctx->Bstate = Bstate;
1222: ctx->Bid = Bid;
1223: MatPropagateSymmetryOptions(B,Br);
1224: }
1225: EPSSetOperators(ctx->eps,Ar,Br);
1226: MatDestroy(&Ar);
1227: MatDestroy(&Br);
1228: }
1229: }
1231: /* Create auxiliary EPS */
1232: if (!ctx->eps) {
1233: EPSCreate(child,&ctx->eps);
1234: EPSGetOptionsPrefix(eps,&prefix);
1235: EPSSetOptionsPrefix(ctx->eps,prefix);
1236: EPSSetOperators(ctx->eps,Ar,Br);
1237: MatDestroy(&Ar);
1238: MatDestroy(&Br);
1239: }
1240: /* Create subcommunicator grouping processes with same rank */
1241: if (!ctx->commset) {
1242: MPI_Comm_rank(child,&rank);
1243: MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank);
1244: ctx->commset = PETSC_TRUE;
1245: }
1246: }
1247: EPSSetType(ctx->eps,((PetscObject)eps)->type_name);
1248: STGetType(eps->st,&sttype);
1249: STSetType(ctx->eps->st,sttype);
1251: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1252: ctx_local->npart = ctx->npart;
1253: ctx_local->global = PETSC_FALSE;
1254: ctx_local->eps = eps;
1255: ctx_local->subc = ctx->subc;
1256: ctx_local->commrank = ctx->commrank;
1257: *childeps = ctx->eps;
1258: PetscFunctionReturn(0);
1259: }
1261: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1262: {
1263: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1264: ST st;
1265: PetscBool isfilt;
1267: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1269: EPSKrylovSchurGetChildEPS(eps,&ctx->eps);
1270: EPSGetST(ctx->eps,&st);
1271: STGetOperator(st,NULL);
1272: STGetKSP(st,ksp);
1273: PetscFunctionReturn(0);
1274: }
1276: /*@
1277: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1278: internal EPS object in case of doing spectrum slicing for a computational interval.
1280: Collective on eps
1282: Input Parameter:
1283: . eps - the eigenproblem solver context
1285: Output Parameter:
1286: . ksp - the internal KSP object
1288: Notes:
1289: When invoked to compute all eigenvalues in an interval with spectrum
1290: slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1291: used to compute eigenvalues by chunks near selected shifts. This function
1292: allows access to the KSP object associated to this internal EPS object.
1294: This function is only available for spectrum slicing runs. In case of
1295: having more than one partition, the returned KSP will be different
1296: in MPI processes belonging to different partitions. Hence, if required,
1297: EPSKrylovSchurSetPartitions() must be called BEFORE this function.
1299: Level: advanced
1301: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1302: @*/
1303: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1304: {
1306: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1307: PetscFunctionReturn(0);
1308: }
1310: PetscErrorCode EPSSetFromOptions_KrylovSchur(PetscOptionItems *PetscOptionsObject,EPS eps)
1311: {
1312: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1313: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1314: PetscReal keep;
1315: PetscInt i,j,k;
1316: KSP ksp;
1318: PetscOptionsHead(PetscOptionsObject,"EPS Krylov-Schur Options");
1320: PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg);
1321: if (flg) EPSKrylovSchurSetRestart(eps,keep);
1323: PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg);
1324: if (flg) EPSKrylovSchurSetLocking(eps,lock);
1326: i = ctx->npart;
1327: PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg);
1328: if (flg) EPSKrylovSchurSetPartitions(eps,i);
1330: b = ctx->detect;
1331: PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg);
1332: if (flg) EPSKrylovSchurSetDetectZeros(eps,b);
1334: i = 1;
1335: j = k = PETSC_DECIDE;
1336: PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1);
1337: PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2);
1338: PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3);
1339: if (f1 || f2 || f3) EPSKrylovSchurSetDimensions(eps,i,j,k);
1341: PetscOptionsTail();
1343: /* set options of child KSP in spectrum slicing */
1344: if (eps->which==EPS_ALL) {
1345: if (!eps->st) EPSGetST(eps,&eps->st);
1346: EPSSetDefaultST(eps);
1347: STSetFromOptions(eps->st); /* need to advance this to check ST type */
1348: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1349: if (!isfilt) {
1350: EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp);
1351: KSPSetFromOptions(ksp);
1352: }
1353: }
1354: PetscFunctionReturn(0);
1355: }
1357: PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1358: {
1359: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1360: PetscBool isascii,isfilt;
1361: KSP ksp;
1362: PetscViewer sviewer;
1364: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1365: if (isascii) {
1366: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1367: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
1368: if (eps->which==EPS_ALL) {
1369: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1370: if (isfilt) PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n");
1371: else {
1372: PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd);
1373: if (ctx->npart>1) {
1374: PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart);
1375: if (ctx->detect) PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n");
1376: }
1377: /* view child KSP */
1378: EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp);
1379: PetscViewerASCIIPushTab(viewer);
1380: if (ctx->npart>1 && ctx->subc) {
1381: PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer);
1382: if (!ctx->subc->color) KSPView(ksp,sviewer);
1383: PetscViewerFlush(sviewer);
1384: PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer);
1385: PetscViewerFlush(viewer);
1386: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1387: PetscViewerASCIIPopSynchronized(viewer);
1388: } else KSPView(ksp,viewer);
1389: PetscViewerASCIIPopTab(viewer);
1390: }
1391: }
1392: }
1393: PetscFunctionReturn(0);
1394: }
1396: PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1397: {
1398: PetscBool isfilt;
1400: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1401: if (eps->which==EPS_ALL && !isfilt) EPSDestroy_KrylovSchur_Slice(eps);
1402: PetscFree(eps->data);
1403: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL);
1404: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL);
1405: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL);
1406: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL);
1407: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL);
1408: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL);
1409: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL);
1410: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL);
1411: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL);
1412: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL);
1413: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL);
1414: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL);
1415: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL);
1416: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL);
1417: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL);
1418: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL);
1419: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL);
1420: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL);
1421: PetscFunctionReturn(0);
1422: }
1424: PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1425: {
1426: PetscBool isfilt;
1428: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1429: if (eps->which==EPS_ALL && !isfilt) EPSReset_KrylovSchur_Slice(eps);
1430: PetscFunctionReturn(0);
1431: }
1433: PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1434: {
1435: if (eps->which==EPS_ALL) {
1436: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSINVERT);
1437: }
1438: PetscFunctionReturn(0);
1439: }
1441: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1442: {
1443: EPS_KRYLOVSCHUR *ctx;
1445: PetscNewLog(eps,&ctx);
1446: eps->data = (void*)ctx;
1447: ctx->lock = PETSC_TRUE;
1448: ctx->nev = 1;
1449: ctx->ncv = PETSC_DEFAULT;
1450: ctx->mpd = PETSC_DEFAULT;
1451: ctx->npart = 1;
1452: ctx->detect = PETSC_FALSE;
1453: ctx->global = PETSC_TRUE;
1455: eps->useds = PETSC_TRUE;
1457: /* solve and computevectors determined at setup */
1458: eps->ops->setup = EPSSetUp_KrylovSchur;
1459: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1460: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1461: eps->ops->destroy = EPSDestroy_KrylovSchur;
1462: eps->ops->reset = EPSReset_KrylovSchur;
1463: eps->ops->view = EPSView_KrylovSchur;
1464: eps->ops->backtransform = EPSBackTransform_Default;
1465: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1467: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur);
1468: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur);
1469: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur);
1470: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur);
1471: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur);
1472: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur);
1473: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur);
1474: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur);
1475: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur);
1476: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur);
1477: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur);
1478: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur);
1479: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur);
1480: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur);
1481: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur);
1482: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur);
1483: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur);
1484: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur);
1485: PetscFunctionReturn(0);
1486: }