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: EPS routines related to problem setup
12: */
14: #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
16: /*
17: Let the solver choose the ST type that should be used by default,
18: otherwise set it to SHIFT.
19: This is called at EPSSetFromOptions (before STSetFromOptions)
20: and also at EPSSetUp (in case EPSSetFromOptions was not called).
21: */
22: PetscErrorCode EPSSetDefaultST(EPS eps) 23: {
27: if (eps->ops->setdefaultst) { (*eps->ops->setdefaultst)(eps); }
28: if (!((PetscObject)eps->st)->type_name) {
29: STSetType(eps->st,STSHIFT);
30: }
31: return(0);
32: }
34: /*
35: This is done by preconditioned eigensolvers that use the PC only.
36: It sets STPRECOND with KSPPREONLY.
37: */
38: PetscErrorCode EPSSetDefaultST_Precond(EPS eps) 39: {
41: KSP ksp;
44: if (!((PetscObject)eps->st)->type_name) {
45: STSetType(eps->st,STPRECOND);
46: STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
47: }
48: STGetKSP(eps->st,&ksp);
49: if (!((PetscObject)ksp)->type_name) {
50: KSPSetType(ksp,KSPPREONLY);
51: }
52: return(0);
53: }
55: /*
56: This is done by preconditioned eigensolvers that can also use the KSP.
57: It sets STPRECOND with the default KSP (GMRES) and maxit=5.
58: */
59: PetscErrorCode EPSSetDefaultST_GMRES(EPS eps) 60: {
62: KSP ksp;
65: if (!((PetscObject)eps->st)->type_name) {
66: STSetType(eps->st,STPRECOND);
67: STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
68: }
69: STGetKSP(eps->st,&ksp);
70: if (!((PetscObject)ksp)->type_name) {
71: KSPSetType(ksp,KSPGMRES);
72: KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
73: }
74: return(0);
75: }
77: /*@
78: EPSSetUp - Sets up all the internal data structures necessary for the
79: execution of the eigensolver. Then calls STSetUp() for any set-up
80: operations associated to the ST object.
82: Collective on EPS 84: Input Parameter:
85: . eps - eigenproblem solver context
87: Notes:
88: This function need not be called explicitly in most cases, since EPSSolve()
89: calls it. It can be useful when one wants to measure the set-up time
90: separately from the solve time.
92: Level: developer
94: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
95: @*/
96: PetscErrorCode EPSSetUp(EPS eps) 97: {
99: Mat A,B;
100: SlepcSC sc;
101: PetscInt k,nmat;
102: PetscBool flg,istrivial,precond;
103: #if defined(PETSC_USE_COMPLEX)
104: PetscScalar sigma;
105: #endif
109: if (eps->state) return(0);
110: PetscLogEventBegin(EPS_SetUp,eps,0,0,0);
112: /* reset the convergence flag from the previous solves */
113: eps->reason = EPS_CONVERGED_ITERATING;
115: /* Set default solver type (EPSSetFromOptions was not called) */
116: if (!((PetscObject)eps)->type_name) {
117: EPSSetType(eps,EPSKRYLOVSCHUR);
118: }
119: if (!eps->st) { EPSGetST(eps,&eps->st); }
120: EPSSetDefaultST(eps);
122: PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond);
123: if (eps->categ==EPS_CATEGORY_PRECOND && !precond) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
124: if (eps->categ!=EPS_CATEGORY_PRECOND && precond) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"STPRECOND is intended for preconditioned eigensolvers only");
126: STSetTransform(eps->st,PETSC_TRUE);
127: if (eps->useds && !eps->ds) { EPSGetDS(eps,&eps->ds); }
128: if (eps->twosided) {
129: if (!eps->hasts) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing left eigenvectors (no two-sided variant)");
130: if (eps->ishermitian) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided methods are not intended for symmetric problems");
131: if (!eps->dsts) { DSDuplicate(eps->ds,&eps->dsts); }
132: }
133: if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
134: if (!((PetscObject)eps->rg)->type_name) {
135: RGSetType(eps->rg,RGINTERVAL);
136: }
138: /* Set problem dimensions */
139: STGetNumMatrices(eps->st,&nmat);
140: if (!nmat) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
141: STMatGetSize(eps->st,&eps->n,NULL);
142: STMatGetLocalSize(eps->st,&eps->nloc,NULL);
144: /* Set default problem type */
145: if (!eps->problem_type) {
146: if (nmat==1) {
147: EPSSetProblemType(eps,EPS_NHEP);
148: } else {
149: EPSSetProblemType(eps,EPS_GNHEP);
150: }
151: } else if (nmat==1 && eps->isgeneralized) {
152: PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n");
153: eps->isgeneralized = PETSC_FALSE;
154: eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
155: } else if (nmat>1 && !eps->isgeneralized) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Inconsistent EPS state");
157: if (eps->nev > eps->n) eps->nev = eps->n;
158: if (eps->ncv > eps->n) eps->ncv = eps->n;
160: /* initialization of matrix norms */
161: if (eps->conv==EPS_CONV_NORM) {
162: if (!eps->nrma) {
163: STGetMatrix(eps->st,0,&A);
164: MatNorm(A,NORM_INFINITY,&eps->nrma);
165: }
166: if (nmat>1 && !eps->nrmb) {
167: STGetMatrix(eps->st,1,&B);
168: MatNorm(B,NORM_INFINITY,&eps->nrmb);
169: }
170: }
172: /* call specific solver setup */
173: (*eps->ops->setup)(eps);
175: /* if purification is set, check that it really makes sense */
176: if (eps->purify) {
177: if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
178: else {
179: if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
180: else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
181: else {
182: PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
183: if (flg) eps->purify = PETSC_FALSE;
184: }
185: }
186: }
188: /* check extraction */
189: PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STPRECOND,STSHIFT,"");
190: if (!flg && eps->extraction && eps->extraction!=EPS_RITZ) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot use a spectral transformation combined with harmonic extraction");
192: /* set tolerance if not yet set */
193: if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL;
195: /* fill sorting criterion context */
196: switch (eps->which) {
197: case EPS_LARGEST_MAGNITUDE:
198: eps->sc->comparison = SlepcCompareLargestMagnitude;
199: eps->sc->comparisonctx = NULL;
200: break;
201: case EPS_SMALLEST_MAGNITUDE:
202: eps->sc->comparison = SlepcCompareSmallestMagnitude;
203: eps->sc->comparisonctx = NULL;
204: break;
205: case EPS_LARGEST_REAL:
206: eps->sc->comparison = SlepcCompareLargestReal;
207: eps->sc->comparisonctx = NULL;
208: break;
209: case EPS_SMALLEST_REAL:
210: eps->sc->comparison = SlepcCompareSmallestReal;
211: eps->sc->comparisonctx = NULL;
212: break;
213: case EPS_LARGEST_IMAGINARY:
214: eps->sc->comparison = SlepcCompareLargestImaginary;
215: eps->sc->comparisonctx = NULL;
216: break;
217: case EPS_SMALLEST_IMAGINARY:
218: eps->sc->comparison = SlepcCompareSmallestImaginary;
219: eps->sc->comparisonctx = NULL;
220: break;
221: case EPS_TARGET_MAGNITUDE:
222: eps->sc->comparison = SlepcCompareTargetMagnitude;
223: eps->sc->comparisonctx = &eps->target;
224: break;
225: case EPS_TARGET_REAL:
226: eps->sc->comparison = SlepcCompareTargetReal;
227: eps->sc->comparisonctx = &eps->target;
228: break;
229: case EPS_TARGET_IMAGINARY:
230: #if defined(PETSC_USE_COMPLEX)
231: eps->sc->comparison = SlepcCompareTargetImaginary;
232: eps->sc->comparisonctx = &eps->target;
233: #endif
234: break;
235: case EPS_ALL:
236: eps->sc->comparison = SlepcCompareSmallestReal;
237: eps->sc->comparisonctx = NULL;
238: break;
239: case EPS_WHICH_USER:
240: break;
241: }
242: eps->sc->map = NULL;
243: eps->sc->mapobj = NULL;
245: /* fill sorting criterion for DS */
246: if (eps->useds && eps->which!=EPS_ALL) {
247: DSGetSlepcSC(eps->ds,&sc);
248: RGIsTrivial(eps->rg,&istrivial);
249: sc->rg = istrivial? NULL: eps->rg;
250: sc->comparison = eps->sc->comparison;
251: sc->comparisonctx = eps->sc->comparisonctx;
252: sc->map = SlepcMap_ST;
253: sc->mapobj = (PetscObject)eps->st;
254: if (eps->twosided) {
255: DSSetSlepcSC(eps->dsts,sc);
256: }
257: }
259: /* Build balancing matrix if required */
260: if (eps->balance!=EPS_BALANCE_USER) {
261: STSetBalanceMatrix(eps->st,NULL);
262: if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
263: if (!eps->D) {
264: BVCreateVec(eps->V,&eps->D);
265: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->D);
266: }
267: EPSBuildBalance_Krylov(eps);
268: STSetBalanceMatrix(eps->st,eps->D);
269: }
270: }
272: /* Setup ST */
273: STSetUp(eps->st);
275: #if defined(PETSC_USE_COMPLEX)
276: STGetShift(eps->st,&sigma);
277: if (eps->ishermitian && PetscImaginaryPart(sigma) != 0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Hermitian problems are not compatible with complex shifts");
278: #endif
279: PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
280: if (flg && eps->problem_type == EPS_PGNHEP) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cayley spectral transformation is not compatible with PGNHEP");
281: PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STCAYLEY,"");
282: if (flg && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER) && (eps->which!=EPS_TARGET_MAGNITUDE && eps->which!=EPS_TARGET_REAL && eps->which!=EPS_TARGET_IMAGINARY && eps->which!=EPS_ALL && eps->which!=EPS_WHICH_USER)) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Shift-and-invert requires a target 'which' (see EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 -eps_target_magnitude");
284: /* process deflation and initial vectors */
285: if (eps->nds<0) {
286: k = -eps->nds;
287: BVInsertConstraints(eps->V,&k,eps->defl);
288: SlepcBasisDestroy_Private(&eps->nds,&eps->defl);
289: eps->nds = k;
290: STCheckNullSpace(eps->st,eps->V);
291: }
292: if (eps->nini<0) {
293: k = -eps->nini;
294: if (k>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The number of initial vectors is larger than ncv");
295: BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE);
296: SlepcBasisDestroy_Private(&eps->nini,&eps->IS);
297: eps->nini = k;
298: }
300: PetscLogEventEnd(EPS_SetUp,eps,0,0,0);
301: eps->state = EPS_STATE_SETUP;
302: return(0);
303: }
305: /*@
306: EPSSetOperators - Sets the matrices associated with the eigenvalue problem.
308: Collective on EPS and Mat
310: Input Parameters:
311: + eps - the eigenproblem solver context
312: . A - the matrix associated with the eigensystem
313: - B - the second matrix in the case of generalized eigenproblems
315: Notes:
316: To specify a standard eigenproblem, use NULL for parameter B.
318: It must be called before EPSSetUp(). If it is called again after EPSSetUp() then
319: the EPS object is reset.
321: Level: beginner
323: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
324: @*/
325: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)326: {
328: PetscInt m,n,m0,nmat;
329: Mat mat[2];
338: /* Check for square matrices */
339: MatGetSize(A,&m,&n);
340: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
341: if (B) {
342: MatGetSize(B,&m0,&n);
343: if (m0!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix");
344: if (m!=m0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match");
345: }
346: if (eps->state && n!=eps->n) { EPSReset(eps); }
347: eps->nrma = 0.0;
348: eps->nrmb = 0.0;
349: if (!eps->st) { EPSGetST(eps,&eps->st); }
350: mat[0] = A;
351: if (B) {
352: mat[1] = B;
353: nmat = 2;
354: } else nmat = 1;
355: STSetMatrices(eps->st,nmat,mat);
356: eps->state = EPS_STATE_INITIAL;
357: return(0);
358: }
360: /*@
361: EPSGetOperators - Gets the matrices associated with the eigensystem.
363: Collective on EPS and Mat
365: Input Parameter:
366: . eps - the EPS context
368: Output Parameters:
369: + A - the matrix associated with the eigensystem
370: - B - the second matrix in the case of generalized eigenproblems
372: Level: intermediate
374: .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
375: @*/
376: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)377: {
379: ST st;
380: PetscInt k;
384: EPSGetST(eps,&st);
385: if (A) { STGetMatrix(st,0,A); }
386: if (B) {
387: STGetNumMatrices(st,&k);
388: if (k==1) B = NULL;
389: else {
390: STGetMatrix(st,1,B);
391: }
392: }
393: return(0);
394: }
396: /*@
397: EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
398: space.
400: Collective on EPS and Vec
402: Input Parameter:
403: + eps - the eigenproblem solver context
404: . n - number of vectors
405: - v - set of basis vectors of the deflation space
407: Notes:
408: When a deflation space is given, the eigensolver seeks the eigensolution
409: in the restriction of the problem to the orthogonal complement of this
410: space. This can be used for instance in the case that an invariant
411: subspace is known beforehand (such as the nullspace of the matrix).
413: These vectors do not persist from one EPSSolve() call to the other, so the
414: deflation space should be set every time.
416: The vectors do not need to be mutually orthonormal, since they are explicitly
417: orthonormalized internally.
419: Level: intermediate
421: .seealso: EPSSetInitialSpace()
422: @*/
423: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec *v)424: {
430: if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n out of range");
431: SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
432: if (n>0) eps->state = EPS_STATE_INITIAL;
433: return(0);
434: }
436: /*@C
437: EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
438: space, that is, the subspace from which the solver starts to iterate.
440: Collective on EPS and Vec
442: Input Parameter:
443: + eps - the eigenproblem solver context
444: . n - number of vectors
445: - is - set of basis vectors of the initial space
447: Notes:
448: Some solvers start to iterate on a single vector (initial vector). In that case,
449: the other vectors are ignored.
451: These vectors do not persist from one EPSSolve() call to the other, so the
452: initial space should be set every time.
454: The vectors do not need to be mutually orthonormal, since they are explicitly
455: orthonormalized internally.
457: Common usage of this function is when the user can provide a rough approximation
458: of the wanted eigenspace. Then, convergence may be faster.
460: Level: intermediate
462: .seealso: EPSSetDeflationSpace()
463: @*/
464: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec *is)465: {
471: if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
472: SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
473: if (n>0) eps->state = EPS_STATE_INITIAL;
474: return(0);
475: }
477: /*
478: EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
479: by the user. This is called at setup.
480: */
481: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)482: {
484: PetscBool krylov;
487: if (*ncv) { /* ncv set */
488: PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
489: if (krylov) {
490: if (*ncv<nev+1 && !(*ncv==nev && *ncv==eps->n)) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev+1");
491: } else {
492: if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
493: }
494: } else if (*mpd) { /* mpd set */
495: *ncv = PetscMin(eps->n,nev+(*mpd));
496: } else { /* neither set: defaults depend on nev being small or large */
497: if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
498: else {
499: *mpd = 500;
500: *ncv = PetscMin(eps->n,nev+(*mpd));
501: }
502: }
503: if (!*mpd) *mpd = *ncv;
504: return(0);
505: }
507: /*@
508: EPSAllocateSolution - Allocate memory storage for common variables such
509: as eigenvalues and eigenvectors.
511: Collective on EPS513: Input Parameters:
514: + eps - eigensolver context
515: - extra - number of additional positions, used for methods that require a
516: working basis slightly larger than ncv
518: Developers Note:
519: This is SLEPC_EXTERN because it may be required by user plugin EPS520: implementations.
522: Level: developer
523: @*/
524: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)525: {
527: PetscInt oldsize,newc,requested;
528: PetscLogDouble cnt;
529: Vec t;
532: requested = eps->ncv + extra;
534: /* oldsize is zero if this is the first time setup is called */
535: BVGetSizes(eps->V,NULL,NULL,&oldsize);
536: newc = PetscMax(0,requested-oldsize);
538: /* allocate space for eigenvalues and friends */
539: if (requested != oldsize || !eps->eigr) {
540: PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
541: PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
542: cnt = 2*newc*sizeof(PetscScalar) + 2*newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
543: PetscLogObjectMemory((PetscObject)eps,cnt);
544: }
546: /* workspace for the case of arbitrary selection */
547: if (eps->arbitrary) {
548: if (eps->rr) {
549: PetscFree2(eps->rr,eps->ri);
550: }
551: PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
552: PetscLogObjectMemory((PetscObject)eps,2*newc*sizeof(PetscScalar));
553: }
555: /* allocate V */
556: if (!eps->V) { EPSGetBV(eps,&eps->V); }
557: if (!oldsize) {
558: if (!((PetscObject)(eps->V))->type_name) {
559: BVSetType(eps->V,BVSVEC);
560: }
561: STMatCreateVecsEmpty(eps->st,&t,NULL);
562: BVSetSizesFromVec(eps->V,t,requested);
563: VecDestroy(&t);
564: } else {
565: BVResize(eps->V,requested,PETSC_FALSE);
566: }
568: /* allocate W */
569: if (eps->twosided) {
570: BVDestroy(&eps->W);
571: BVDuplicate(eps->V,&eps->W);
572: }
573: return(0);
574: }