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: PEP routines related to problem setup
12: */
14: #include <slepc/private/pepimpl.h> /*I "slepcpep.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 PEPSetFromOptions (before STSetFromOptions)
20: and also at PEPSetUp (in case PEPSetFromOptions was not called).
21: */
22: PetscErrorCode PEPSetDefaultST(PEP pep) 23: {
27: if (pep->ops->setdefaultst) { (*pep->ops->setdefaultst)(pep); }
28: if (!((PetscObject)pep->st)->type_name) {
29: STSetType(pep->st,STSHIFT);
30: }
31: return(0);
32: }
34: /*
35: This is used in Q-Arnoldi and STOAR to set the transform flag by
36: default, otherwise the user has to explicitly run with -st_transform
37: */
38: PetscErrorCode PEPSetDefaultST_Transform(PEP pep) 39: {
43: STSetTransform(pep->st,PETSC_TRUE);
44: return(0);
45: }
47: /*@
48: PEPSetUp - Sets up all the internal data structures necessary for the
49: execution of the PEP solver.
51: Collective on PEP 53: Input Parameter:
54: . pep - solver context
56: Notes:
57: This function need not be called explicitly in most cases, since PEPSolve()
58: calls it. It can be useful when one wants to measure the set-up time
59: separately from the solve time.
61: Level: developer
63: .seealso: PEPCreate(), PEPSolve(), PEPDestroy()
64: @*/
65: PetscErrorCode PEPSetUp(PEP pep) 66: {
68: SlepcSC sc;
69: PetscBool istrivial,flg;
70: PetscInt k;
71: KSP ksp;
72: PC pc;
73: PetscMPIInt size;
74: MatSolverType stype;
78: if (pep->state) return(0);
79: PetscLogEventBegin(PEP_SetUp,pep,0,0,0);
81: /* reset the convergence flag from the previous solves */
82: pep->reason = PEP_CONVERGED_ITERATING;
84: /* set default solver type (PEPSetFromOptions was not called) */
85: if (!((PetscObject)pep)->type_name) {
86: PEPSetType(pep,PEPTOAR);
87: }
88: if (!pep->st) { PEPGetST(pep,&pep->st); }
89: PEPSetDefaultST(pep);
90: if (!pep->ds) { PEPGetDS(pep,&pep->ds); }
91: if (!pep->rg) { PEPGetRG(pep,&pep->rg); }
92: if (!((PetscObject)pep->rg)->type_name) {
93: RGSetType(pep->rg,RGINTERVAL);
94: }
96: /* check matrices, transfer them to ST */
97: if (!pep->A) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"PEPSetOperators must be called first");
98: STSetMatrices(pep->st,pep->nmat,pep->A);
100: /* set problem dimensions */
101: MatGetSize(pep->A[0],&pep->n,NULL);
102: MatGetLocalSize(pep->A[0],&pep->nloc,NULL);
104: /* set default problem type */
105: if (!pep->problem_type) {
106: PEPSetProblemType(pep,PEP_GENERAL);
107: }
108: if (pep->nev > (pep->nmat-1)*pep->n) pep->nev = (pep->nmat-1)*pep->n;
109: if (pep->ncv > (pep->nmat-1)*pep->n) pep->ncv = (pep->nmat-1)*pep->n;
111: /* check consistency of refinement options */
112: if (pep->refine) {
113: if (!pep->scheme) { /* set default scheme */
114: PEPRefineGetKSP(pep,&ksp);
115: KSPGetPC(ksp,&pc);
116: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
117: if (flg) {
118: PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
119: }
120: pep->scheme = flg? PEP_REFINE_SCHEME_MBE: PEP_REFINE_SCHEME_SCHUR;
121: }
122: if (pep->scheme==PEP_REFINE_SCHEME_MBE) {
123: PEPRefineGetKSP(pep,&ksp);
124: KSPGetPC(ksp,&pc);
125: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
126: if (flg) {
127: PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
128: }
129: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
130: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
131: if (size>1) { /* currently selected PC is a factorization */
132: PCFactorGetMatSolverType(pc,&stype);
133: PetscStrcmp(stype,MATSOLVERPETSC,&flg);
134: if (flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"For Newton refinement, you chose to solve linear systems with a factorization, but in parallel runs you need to select an external package");
135: }
136: }
137: if (pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
138: if (pep->npart>1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
139: }
140: }
141: /* call specific solver setup */
142: (*pep->ops->setup)(pep);
144: /* set tolerance if not yet set */
145: if (pep->tol==PETSC_DEFAULT) pep->tol = SLEPC_DEFAULT_TOL;
146: if (pep->refine) {
147: if (pep->rtol==PETSC_DEFAULT) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
148: if (pep->rits==PETSC_DEFAULT) pep->rits = (pep->refine==PEP_REFINE_SIMPLE)? 10: 1;
149: }
151: /* set default extraction */
152: if (!pep->extract) {
153: pep->extract = (pep->basis==PEP_BASIS_MONOMIAL)? PEP_EXTRACT_NORM: PEP_EXTRACT_NONE;
154: }
156: /* fill sorting criterion context */
157: switch (pep->which) {
158: case PEP_LARGEST_MAGNITUDE:
159: pep->sc->comparison = SlepcCompareLargestMagnitude;
160: pep->sc->comparisonctx = NULL;
161: break;
162: case PEP_SMALLEST_MAGNITUDE:
163: pep->sc->comparison = SlepcCompareSmallestMagnitude;
164: pep->sc->comparisonctx = NULL;
165: break;
166: case PEP_LARGEST_REAL:
167: pep->sc->comparison = SlepcCompareLargestReal;
168: pep->sc->comparisonctx = NULL;
169: break;
170: case PEP_SMALLEST_REAL:
171: pep->sc->comparison = SlepcCompareSmallestReal;
172: pep->sc->comparisonctx = NULL;
173: break;
174: case PEP_LARGEST_IMAGINARY:
175: pep->sc->comparison = SlepcCompareLargestImaginary;
176: pep->sc->comparisonctx = NULL;
177: break;
178: case PEP_SMALLEST_IMAGINARY:
179: pep->sc->comparison = SlepcCompareSmallestImaginary;
180: pep->sc->comparisonctx = NULL;
181: break;
182: case PEP_TARGET_MAGNITUDE:
183: pep->sc->comparison = SlepcCompareTargetMagnitude;
184: pep->sc->comparisonctx = &pep->target;
185: break;
186: case PEP_TARGET_REAL:
187: pep->sc->comparison = SlepcCompareTargetReal;
188: pep->sc->comparisonctx = &pep->target;
189: break;
190: case PEP_TARGET_IMAGINARY:
191: #if defined(PETSC_USE_COMPLEX)
192: pep->sc->comparison = SlepcCompareTargetImaginary;
193: pep->sc->comparisonctx = &pep->target;
194: #endif
195: break;
196: case PEP_ALL:
197: pep->sc->comparison = SlepcCompareSmallestReal;
198: pep->sc->comparisonctx = NULL;
199: break;
200: case PEP_WHICH_USER:
201: break;
202: }
203: pep->sc->map = NULL;
204: pep->sc->mapobj = NULL;
206: /* fill sorting criterion for DS */
207: if (pep->which!=PEP_ALL) {
208: DSGetSlepcSC(pep->ds,&sc);
209: RGIsTrivial(pep->rg,&istrivial);
210: sc->rg = istrivial? NULL: pep->rg;
211: sc->comparison = pep->sc->comparison;
212: sc->comparisonctx = pep->sc->comparisonctx;
213: sc->map = SlepcMap_ST;
214: sc->mapobj = (PetscObject)pep->st;
215: }
216: /* setup ST */
217: STSetUp(pep->st);
219: /* compute matrix coefficients */
220: STGetTransform(pep->st,&flg);
221: if (!flg) {
222: if (pep->solvematcoeffs) { STMatSetUp(pep->st,1.0,pep->solvematcoeffs); }
223: } else {
224: if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Cannot use ST-transform with non-monomial basis in PEP");
225: }
227: /* compute scale factor if no set by user */
228: PEPComputeScaleFactor(pep);
230: /* build balancing matrix if required */
231: if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) {
232: if (!pep->Dl) {
233: BVCreateVec(pep->V,&pep->Dl);
234: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->Dl);
235: }
236: if (!pep->Dr) {
237: BVCreateVec(pep->V,&pep->Dr);
238: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->Dr);
239: }
240: PEPBuildDiagonalScaling(pep);
241: }
243: /* process initial vectors */
244: if (pep->nini<0) {
245: k = -pep->nini;
246: if (k>pep->ncv) SETERRQ(PetscObjectComm((PetscObject)pep),1,"The number of initial vectors is larger than ncv");
247: BVInsertVecs(pep->V,0,&k,pep->IS,PETSC_TRUE);
248: SlepcBasisDestroy_Private(&pep->nini,&pep->IS);
249: pep->nini = k;
250: }
251: PetscLogEventEnd(PEP_SetUp,pep,0,0,0);
252: pep->state = PEP_STATE_SETUP;
253: return(0);
254: }
256: /*@
257: PEPSetOperators - Sets the coefficient matrices associated with the polynomial
258: eigenvalue problem.
260: Collective on PEP and Mat
262: Input Parameters:
263: + pep - the eigenproblem solver context
264: . nmat - number of matrices in array A
265: - A - the array of matrices associated with the eigenproblem
267: Notes:
268: The polynomial eigenproblem is defined as P(l)*x=0, where l is
269: the eigenvalue, x is the eigenvector, and P(l) is defined as
270: P(l) = A_0 + l*A_1 + ... + l^d*A_d, with d=nmat-1 (the degree of P).
271: For non-monomial bases, this expression is different.
273: Level: beginner
275: .seealso: PEPSolve(), PEPGetOperators(), PEPGetNumMatrices(), PEPSetBasis()
276: @*/
277: PetscErrorCode PEPSetOperators(PEP pep,PetscInt nmat,Mat A[])278: {
280: PetscInt i,n=0,m,m0=0;
285: if (nmat <= 0) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Non-positive value of nmat: %D",nmat);
286: if (nmat <= 2) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Cannot solve linear eigenproblems with PEP; use EPS instead");
289: for (i=0;i<nmat;i++) {
292: MatGetSize(A[i],&m,&n);
293: if (m!=n) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"A[%D] is a non-square matrix",i);
294: if (!i) m0 = m;
295: if (m!=m0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_INCOMP,"Dimensions of matrices do not match with each other");
296: PetscObjectReference((PetscObject)A[i]);
297: }
299: if (pep->state && n!=pep->n) { PEPReset(pep); }
300: else if (pep->nmat) {
301: MatDestroyMatrices(pep->nmat,&pep->A);
302: PetscFree2(pep->pbc,pep->nrma);
303: PetscFree(pep->solvematcoeffs);
304: }
306: PetscMalloc1(nmat,&pep->A);
307: PetscCalloc2(3*nmat,&pep->pbc,nmat,&pep->nrma);
308: PetscLogObjectMemory((PetscObject)pep,nmat*sizeof(Mat)+4*nmat*sizeof(PetscReal));
309: for (i=0;i<nmat;i++) {
310: pep->A[i] = A[i];
311: pep->pbc[i] = 1.0; /* default to monomial basis */
312: }
313: pep->nmat = nmat;
314: pep->state = PEP_STATE_INITIAL;
315: return(0);
316: }
318: /*@
319: PEPGetOperators - Gets the matrices associated with the polynomial eigensystem.
321: Not collective, though parallel Mats are returned if the PEP is parallel
323: Input Parameters:
324: + pep - the PEP context
325: - k - the index of the requested matrix (starting in 0)
327: Output Parameter:
328: . A - the requested matrix
330: Level: intermediate
332: .seealso: PEPSolve(), PEPSetOperators(), PEPGetNumMatrices()
333: @*/
334: PetscErrorCode PEPGetOperators(PEP pep,PetscInt k,Mat *A)335: {
339: if (k<0 || k>=pep->nmat) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",pep->nmat-1);
340: *A = pep->A[k];
341: return(0);
342: }
344: /*@
345: PEPGetNumMatrices - Returns the number of matrices stored in the PEP.
347: Not collective
349: Input Parameter:
350: . pep - the PEP context
352: Output Parameters:
353: . nmat - the number of matrices passed in PEPSetOperators()
355: Level: intermediate
357: .seealso: PEPSetOperators()
358: @*/
359: PetscErrorCode PEPGetNumMatrices(PEP pep,PetscInt *nmat)360: {
364: *nmat = pep->nmat;
365: return(0);
366: }
368: /*@C
369: PEPSetInitialSpace - Specify a basis of vectors that constitute the initial
370: space, that is, the subspace from which the solver starts to iterate.
372: Collective on PEP and Vec
374: Input Parameter:
375: + pep - the polynomial eigensolver context
376: . n - number of vectors
377: - is - set of basis vectors of the initial space
379: Notes:
380: Some solvers start to iterate on a single vector (initial vector). In that case,
381: the other vectors are ignored.
383: These vectors do not persist from one PEPSolve() call to the other, so the
384: initial space should be set every time.
386: The vectors do not need to be mutually orthonormal, since they are explicitly
387: orthonormalized internally.
389: Common usage of this function is when the user can provide a rough approximation
390: of the wanted eigenspace. Then, convergence may be faster.
392: Level: intermediate
393: @*/
394: PetscErrorCode PEPSetInitialSpace(PEP pep,PetscInt n,Vec *is)395: {
401: if (n<0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
402: SlepcBasisReference_Private(n,is,&pep->nini,&pep->IS);
403: if (n>0) pep->state = PEP_STATE_INITIAL;
404: return(0);
405: }
407: /*
408: PEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
409: by the user. This is called at setup.
410: */
411: PetscErrorCode PEPSetDimensions_Default(PEP pep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)412: {
414: PetscBool krylov;
415: PetscInt dim;
418: PetscObjectTypeCompareAny((PetscObject)pep,&krylov,PEPTOAR,PEPSTOAR,PEPQARNOLDI,"");
419: dim = (pep->nmat-1)*pep->n;
420: if (*ncv) { /* ncv set */
421: if (krylov) {
422: if (*ncv<nev+1 && !(*ncv==nev && *ncv==dim)) SETERRQ(PetscObjectComm((PetscObject)pep),1,"The value of ncv must be at least nev+1");
423: } else {
424: if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)pep),1,"The value of ncv must be at least nev");
425: }
426: } else if (*mpd) { /* mpd set */
427: *ncv = PetscMin(dim,nev+(*mpd));
428: } else { /* neither set: defaults depend on nev being small or large */
429: if (nev<500) *ncv = PetscMin(dim,PetscMax(2*nev,nev+15));
430: else {
431: *mpd = 500;
432: *ncv = PetscMin(dim,nev+(*mpd));
433: }
434: }
435: if (!*mpd) *mpd = *ncv;
436: return(0);
437: }
439: /*@
440: PEPAllocateSolution - Allocate memory storage for common variables such
441: as eigenvalues and eigenvectors.
443: Collective on PEP445: Input Parameters:
446: + pep - eigensolver context
447: - extra - number of additional positions, used for methods that require a
448: working basis slightly larger than ncv
450: Developers Note:
451: This is SLEPC_EXTERN because it may be required by user plugin PEP452: implementations.
454: Level: developer
455: @*/
456: PetscErrorCode PEPAllocateSolution(PEP pep,PetscInt extra)457: {
459: PetscInt oldsize,newc,requested,requestedbv;
460: PetscLogDouble cnt;
461: Vec t;
464: requested = (pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1)) + extra;
465: requestedbv = pep->ncv + extra;
467: /* oldsize is zero if this is the first time setup is called */
468: BVGetSizes(pep->V,NULL,NULL,&oldsize);
470: /* allocate space for eigenvalues and friends */
471: if (requested != oldsize || !pep->eigr) {
472: PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm);
473: PetscMalloc4(requested,&pep->eigr,requested,&pep->eigi,requested,&pep->errest,requested,&pep->perm);
474: newc = PetscMax(0,requested-oldsize);
475: cnt = 2*newc*sizeof(PetscScalar) + newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
476: PetscLogObjectMemory((PetscObject)pep,cnt);
477: }
479: /* allocate V */
480: if (!pep->V) { PEPGetBV(pep,&pep->V); }
481: if (!oldsize) {
482: if (!((PetscObject)(pep->V))->type_name) {
483: BVSetType(pep->V,BVSVEC);
484: }
485: STMatCreateVecsEmpty(pep->st,&t,NULL);
486: BVSetSizesFromVec(pep->V,t,requestedbv);
487: VecDestroy(&t);
488: } else {
489: BVResize(pep->V,requestedbv,PETSC_FALSE);
490: }
491: return(0);
492: }