Actual source code: epssolve.c
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: EPS routines related to the solution process
12: */
14: #include <slepc/private/epsimpl.h>
15: #include <slepc/private/bvimpl.h>
16: #include <petscdraw.h>
18: PetscErrorCode EPSComputeVectors(EPS eps)
19: {
20: PetscFunctionBegin;
21: EPSCheckSolved(eps,1);
22: if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
23: eps->state = EPS_STATE_EIGENVECTORS;
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode EPSComputeValues(EPS eps)
28: {
29: PetscBool injective,iscomp,isfilter;
30: PetscInt i,n,aux,nconv0;
31: Mat A,B=NULL,G,Z;
33: PetscFunctionBegin;
34: switch (eps->categ) {
35: case EPS_CATEGORY_KRYLOV:
36: case EPS_CATEGORY_OTHER:
37: PetscCall(STIsInjective(eps->st,&injective));
38: if (injective) {
39: /* one-to-one mapping: backtransform eigenvalues */
40: PetscUseTypeMethod(eps,backtransform);
41: } else {
42: /* compute eigenvalues from Rayleigh quotient */
43: PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
44: if (!n) break;
45: PetscCall(EPSGetOperators(eps,&A,&B));
46: PetscCall(BVSetActiveColumns(eps->V,0,n));
47: PetscCall(DSGetCompact(eps->ds,&iscomp));
48: PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
49: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
50: PetscCall(BVMatProject(eps->V,A,eps->V,G));
51: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
52: if (B) {
53: PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
54: PetscCall(BVMatProject(eps->V,B,eps->V,G));
55: PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
56: }
57: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
58: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
59: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
60: PetscCall(DSSetCompact(eps->ds,iscomp));
61: if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
62: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
63: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
64: PetscCall(BVMultInPlace(eps->V,Z,0,n));
65: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
66: }
67: /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
68: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
69: if (isfilter) {
70: nconv0 = eps->nconv;
71: for (i=0;i<eps->nconv;i++) {
72: if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
73: eps->nconv--;
74: if (i<eps->nconv) { SlepcSwap(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
75: }
76: }
77: if (nconv0>eps->nconv) PetscCall(PetscInfo(eps,"Discarded %" PetscInt_FMT " computed eigenvalues lying outside the interval\n",nconv0-eps->nconv));
78: }
79: }
80: break;
81: case EPS_CATEGORY_PRECOND:
82: case EPS_CATEGORY_CONTOUR:
83: /* eigenvalues already available as an output of the solver */
84: break;
85: }
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*@
90: EPSSolve - Solves the eigensystem.
92: Collective
94: Input Parameter:
95: . eps - the linear eigensolver context
97: Options Database Keys:
98: + -eps_view - print information about the solver once the solve is complete
99: . -eps_view_pre - print information about the solver before the solve starts
100: . -eps_view_mat0 - view the first matrix ($A$)
101: . -eps_view_mat1 - view the second matrix ($B$)
102: . -eps_view_vectors - view the computed eigenvectors
103: . -eps_view_values - view the computed eigenvalues
104: . -eps_converged_reason - print reason for convergence/divergence, and number of iterations
105: . -eps_error_absolute - print absolute errors of each eigenpair
106: . -eps_error_relative - print relative errors of each eigenpair
107: - -eps_error_backward - print backward errors of each eigenpair
109: Notes:
110: The problem matrices are specified with `EPSSetOperators()`.
112: `EPSSolve()` will return without generating an error regardless of whether
113: all requested solutions were computed or not. Call `EPSGetConverged()` to get the
114: actual number of computed solutions, and `EPSGetConvergedReason()` to determine if
115: the solver converged or failed and why.
117: All the command-line options listed above admit an optional argument specifying
118: the viewer type and options. For instance, use `-eps_view_mat0 binary:amatrix.bin`
119: to save the $A$ matrix to a binary file, `-eps_view_values draw` to draw the computed
120: eigenvalues graphically, or `-eps_error_relative :myerr.m:ascii_matlab` to save
121: the errors in a file that can be executed in Matlab.
122: See `PetscObjectViewFromOptions()` for more details.
124: Level: beginner
126: .seealso: [](ch:eps), `EPSCreate()`, `EPSSetUp()`, `EPSDestroy()`, `EPSSetOperators()`, `EPSGetConverged()`, `EPSGetConvergedReason()`
127: @*/
128: PetscErrorCode EPSSolve(EPS eps)
129: {
130: PetscInt i;
131: PetscBool hasname;
132: STMatMode matmode;
133: Mat A,B;
135: PetscFunctionBegin;
137: if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
138: PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
140: /* Call setup */
141: PetscCall(EPSSetUp(eps));
143: /* Safeguard for matrices of size 0 */
144: if (eps->n == 0) {
145: eps->nconv = 0;
146: eps->reason = EPS_CONVERGED_TOL;
147: eps->state = EPS_STATE_SOLVED;
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: eps->nconv = 0;
152: eps->its = 0;
153: for (i=0;i<eps->ncv;i++) {
154: eps->eigr[i] = 0.0;
155: eps->eigi[i] = 0.0;
156: eps->errest[i] = 0.0;
157: eps->perm[i] = i;
158: }
159: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
160: PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
162: /* Call solver */
163: PetscUseTypeMethod(eps,solve);
164: PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
165: eps->state = EPS_STATE_SOLVED;
167: /* Only the first nconv columns contain useful information (except in CISS) */
168: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
169: if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
171: /* If inplace, purify eigenvectors before reverting operator */
172: PetscCall(STGetMatMode(eps->st,&matmode));
173: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
174: PetscCall(STPostSolve(eps->st));
176: /* Map eigenvalues back to the original problem if appropriate */
177: PetscCall(EPSComputeValues(eps));
179: #if !defined(PETSC_USE_COMPLEX)
180: /* Reorder conjugate eigenvalues (positive imaginary first) */
181: for (i=0;i<eps->nconv-1;i++) {
182: if (eps->eigi[i] != 0 && (eps->problem_type!=EPS_HAMILT || eps->eigr[i]!=0)) {
183: /* conjugate eigenvalues */
184: if (eps->eigi[i] < 0) {
185: eps->eigi[i] = -eps->eigi[i];
186: eps->eigi[i+1] = -eps->eigi[i+1];
187: /* the next correction only works with eigenvectors */
188: PetscCall(EPSComputeVectors(eps));
189: PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
190: if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
191: }
192: i++;
193: }
194: }
195: #endif
197: /* Sort eigenvalues according to eps->which parameter */
198: if (eps->problem_type==EPS_HAMILT) PetscCall(SlepcSortEigenvaluesSpecial(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
199: else PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
200: PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
202: /* Various viewers */
203: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
204: PetscCall(EPSConvergedReasonViewFromOptions(eps));
205: PetscCall(EPSErrorViewFromOptions(eps));
206: PetscCall(EPSValuesViewFromOptions(eps));
207: PetscCall(EPSVectorsViewFromOptions(eps));
209: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
210: if (hasname) {
211: PetscCall(EPSGetOperators(eps,&A,NULL));
212: PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
213: }
214: if (eps->isgeneralized) {
215: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
216: if (hasname) {
217: PetscCall(EPSGetOperators(eps,NULL,&B));
218: PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
219: }
220: }
222: /* Remove deflation and initial subspaces */
223: if (eps->nds) {
224: PetscCall(BVSetNumConstraints(eps->V,0));
225: eps->nds = 0;
226: }
227: eps->nini = 0;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: EPSGetIterationNumber - Gets the current iteration number. If the
233: call to `EPSSolve()` is complete, then it returns the number of iterations
234: carried out by the solution method.
236: Not Collective
238: Input Parameter:
239: . eps - the linear eigensolver context
241: Output Parameter:
242: . its - number of iterations
244: Note:
245: During the $i$-th iteration this call returns $i-1$. If `EPSSolve()` is
246: complete, then parameter `its` contains either the iteration number at
247: which convergence was successfully reached, or failure was detected.
248: Call `EPSGetConvergedReason()` to determine if the solver converged or
249: failed and why.
251: Level: intermediate
253: .seealso: [](ch:eps), `EPSGetConvergedReason()`, `EPSSetTolerances()`
254: @*/
255: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
256: {
257: PetscFunctionBegin;
259: PetscAssertPointer(its,2);
260: *its = eps->its;
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@
265: EPSGetConverged - Gets the number of converged eigenpairs.
267: Not Collective
269: Input Parameter:
270: . eps - the linear eigensolver context
272: Output Parameter:
273: . nconv - number of converged eigenpairs
275: Notes:
276: This function should be called after `EPSSolve()` has finished.
278: The value `nconv` may be different from the number of requested solutions
279: `nev`, but not larger than `ncv`, see `EPSSetDimensions()`.
281: Level: beginner
283: .seealso: [](ch:eps), `EPSSetDimensions()`, `EPSSolve()`, `EPSGetEigenpair()`
284: @*/
285: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
286: {
287: PetscFunctionBegin;
289: PetscAssertPointer(nconv,2);
290: EPSCheckSolved(eps,1);
291: PetscCall(EPS_GetActualConverged(eps,nconv));
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@
296: EPSGetConvergedReason - Gets the reason why the `EPSSolve()` iteration was
297: stopped.
299: Not Collective
301: Input Parameter:
302: . eps - the linear eigensolver context
304: Output Parameter:
305: . reason - negative value indicates diverged, positive value converged, see
306: `EPSConvergedReason` for the possible values
308: Options Database Key:
309: . -eps_converged_reason - print reason for convergence/divergence, and number of iterations
311: Note:
312: If this routine is called before or doing the `EPSSolve()` the value of
313: `EPS_CONVERGED_ITERATING` is returned.
315: Level: intermediate
317: .seealso: [](ch:eps), `EPSSetTolerances()`, `EPSSolve()`, `EPSConvergedReason`
318: @*/
319: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
320: {
321: PetscFunctionBegin;
323: PetscAssertPointer(reason,2);
324: EPSCheckSolved(eps,1);
325: *reason = eps->reason;
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: /*@
330: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
331: subspace.
333: Collective
335: Input Parameter:
336: . eps - the linear eigensolver context
338: Output Parameter:
339: . v - an array of vectors
341: Notes:
342: This function should be called after `EPSSolve()` has finished.
344: The user should provide in `v` an array of `nconv` vectors, where `nconv`
345: is the value returned by `EPSGetConverged()`.
347: The first $k$ vectors returned in `v` span an invariant subspace associated
348: with the first $k$ computed eigenvalues (note that this is not true if the
349: $k$-th eigenvalue is complex and matrix $A$ is real; in this case the first
350: $k+1$ vectors should be used). An invariant subspace $X$ of $A$ satisfies
351: $Ax\in X, \forall x\in X$ (a similar definition applies for generalized
352: eigenproblems).
354: Level: intermediate
356: .seealso: [](ch:eps), `EPSGetEigenpair()`, `EPSGetConverged()`, `EPSSolve()`
357: @*/
358: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
359: {
360: PetscInt i;
361: BV V=eps->V;
362: Vec w;
364: PetscFunctionBegin;
366: PetscAssertPointer(v,2);
368: EPSCheckSolved(eps,1);
369: PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
370: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
371: PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
372: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
373: PetscCall(BVCopy(eps->V,V));
374: for (i=0;i<eps->nconv;i++) {
375: PetscCall(BVGetColumn(V,i,&w));
376: PetscCall(VecPointwiseDivide(w,w,eps->D));
377: PetscCall(BVRestoreColumn(V,i,&w));
378: }
379: PetscCall(BVOrthogonalize(V,NULL));
380: }
381: for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
382: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: EPSGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
388: `EPSSolve()`. The solution consists in both the eigenvalue and the eigenvector.
390: Collective
392: Input Parameters:
393: + eps - the linear eigensolver context
394: - i - index of the solution
396: Output Parameters:
397: + eigr - real part of eigenvalue
398: . eigi - imaginary part of eigenvalue
399: . Vr - real part of eigenvector
400: - Vi - imaginary part of eigenvector
402: Notes:
403: It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
404: required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
405: they must be created by the calling program with e.g. `MatCreateVecs()`.
407: If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
408: configured with complex scalars the eigenvalue is stored
409: directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi` is
410: set to zero). In both cases, the user can pass `NULL` in `eigi` and `Vi`.
412: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
413: Eigenpairs are indexed according to the ordering criterion established
414: with `EPSSetWhichEigenpairs()`.
416: The 2-norm of the eigenvector is one unless the problem is generalized
417: Hermitian. In this case the eigenvector is normalized with respect to the
418: norm defined by the $B$ matrix.
420: In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
421: [](#sec:structured-vectors).
423: Level: beginner
425: .seealso: [](ch:eps), `EPSGetEigenvalue()`, `EPSGetEigenvector()`, `EPSGetLeftEigenvector()`, `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetInvariantSubspace()`
426: @*/
427: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
428: {
429: PetscInt nconv;
431: PetscFunctionBegin;
434: EPSCheckSolved(eps,1);
435: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
436: PetscCall(EPS_GetActualConverged(eps,&nconv));
437: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
438: PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
439: if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@
444: EPSGetEigenvalue - Gets the `i`-th eigenvalue as computed by `EPSSolve()`.
446: Not Collective
448: Input Parameters:
449: + eps - the linear eigensolver context
450: - i - index of the solution
452: Output Parameters:
453: + eigr - real part of eigenvalue
454: - eigi - imaginary part of eigenvalue
456: Notes:
457: If the eigenvalue is real, then `eigi` is set to zero. If PETSc is
458: configured with complex scalars the eigenvalue is stored
459: directly in `eigr` (`eigi` is set to zero).
461: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
462: Eigenpairs are indexed according to the ordering criterion established
463: with `EPSSetWhichEigenpairs()`.
465: Level: beginner
467: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`
468: @*/
469: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
470: {
471: PetscInt k,nconv;
472: #if !defined(PETSC_USE_COMPLEX)
473: PetscInt k2, iquad;
474: #endif
476: PetscFunctionBegin;
478: EPSCheckSolved(eps,1);
479: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
480: PetscCall(EPS_GetActualConverged(eps,&nconv));
481: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
482: if (nconv==eps->nconv) {
483: k = eps->perm[i];
484: #if defined(PETSC_USE_COMPLEX)
485: if (eigr) *eigr = eps->eigr[k];
486: if (eigi) *eigi = 0;
487: #else
488: if (eigr) *eigr = eps->eigr[k];
489: if (eigi) *eigi = eps->eigi[k];
490: #endif
491: } else {
492: PetscCheck(eps->problem_type==EPS_BSE || eps->problem_type==EPS_HAMILT || eps->problem_type==EPS_LREP,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE, Hamiltonian, or LREP");
493: if (eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP) {
494: /* BSE problem, even index is +lambda, odd index is -lambda */
495: k = eps->perm[i/2];
496: #if defined(PETSC_USE_COMPLEX)
497: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
498: if (eigi) *eigi = 0;
499: #else
500: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
501: if (eigi) *eigi = eps->eigi[k];
502: #endif
503: } else if (eps->problem_type==EPS_HAMILT) {
504: /* Hamiltonian eigenproblem */
505: k = eps->perm[i/2];
506: #if defined(PETSC_USE_COMPLEX)
507: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
508: if (eigi) *eigi = 0;
509: #else
510: if (eps->eigi[k]==0.0) { /* real eigenvalue */
511: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
512: if (eigi) *eigi = 0.0;
513: } else if (eps->eigr[k]==0.0) { /* purely imaginary eigenvalue */
514: if (eigr) *eigr = 0.0;
515: if (eigi) *eigi = (i%2)? -eps->eigi[k]: eps->eigi[k];
516: } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
517: iquad = i%2; /* index within the 4 values */
518: if (i>1) {
519: k2 = eps->perm[(i-2)/2];
520: if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
521: }
522: if (eigr) *eigr = (iquad<2)? -eps->eigr[k]: eps->eigr[k];
523: if (eigi) *eigi = (iquad%3)? -eps->eigi[k]: eps->eigi[k];
524: }
525: #endif
526: }
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: /*@
532: EPSGetEigenvector - Gets the `i`-th right eigenvector as computed by `EPSSolve()`.
534: Collective
536: Input Parameters:
537: + eps - the linear eigensolver context
538: - i - index of the solution
540: Output Parameters:
541: + Vr - real part of eigenvector
542: - Vi - imaginary part of eigenvector
544: Notes:
545: The caller must provide valid `Vec` objects, i.e., they must be created
546: by the calling program with e.g. `MatCreateVecs()`.
548: If the corresponding eigenvalue is real, then `Vi` is set to zero. If PETSc is
549: configured with complex scalars the eigenvector is stored
550: directly in `Vr` (`Vi` is set to zero). In any case, the user can pass `NULL` in `Vr`
551: or `Vi` if one of them is not required.
553: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
554: Eigenpairs are indexed according to the ordering criterion established
555: with `EPSSetWhichEigenpairs()`.
557: The 2-norm of the eigenvector is one unless the problem is generalized
558: Hermitian. In this case the eigenvector is normalized with respect to the
559: norm defined by the $B$ matrix.
561: In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
562: [](#sec:structured-vectors).
564: Level: beginner
566: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`, `EPSGetLeftEigenvector()`
567: @*/
568: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
569: {
570: PetscInt nconv;
572: PetscFunctionBegin;
577: EPSCheckSolved(eps,1);
578: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
579: PetscCall(EPS_GetActualConverged(eps,&nconv));
580: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
581: PetscCall(EPSComputeVectors(eps));
582: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: /*@
587: EPSGetLeftEigenvector - Gets the `i`-th left eigenvector as computed by `EPSSolve()`.
589: Collective
591: Input Parameters:
592: + eps - the linear eigensolver context
593: - i - index of the solution
595: Output Parameters:
596: + Wr - real part of left eigenvector
597: - Wi - imaginary part of left eigenvector
599: Notes:
600: The caller must provide valid `Vec` objects, i.e., they must be created
601: by the calling program with e.g. `MatCreateVecs()`.
603: If the corresponding eigenvalue is real, then `Wi` is set to zero. If PETSc is
604: configured with complex scalars the eigenvector is stored directly in `Wr`
605: (`Wi` is set to zero). In any case, the user can pass `NULL` in `Wr` or `Wi` if
606: one of them is not required.
608: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
609: Eigensolutions are indexed according to the ordering criterion established
610: with `EPSSetWhichEigenpairs()`.
612: Left eigenvectors are available only if the `twosided` flag was set, see
613: `EPSSetTwoSided()`.
615: In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
616: [](#sec:structured-vectors).
618: Level: intermediate
620: .seealso: [](ch:eps), `EPSGetEigenvector()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSSetTwoSided()`
621: @*/
622: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
623: {
624: PetscInt nconv;
625: PetscBool trivial,lrepreduced=PETSC_FALSE;
626: Mat H;
627: IS is[2];
628: Vec v0,v1,z,z0,z1;
630: PetscFunctionBegin;
635: EPSCheckSolved(eps,1);
636: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
637: PetscCall(EPS_GetActualConverged(eps,&nconv));
638: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
640: trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP)? PETSC_TRUE: PETSC_FALSE;
641: if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
643: PetscCall(EPSComputeVectors(eps));
644: if (trivial) {
645: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
646: PetscCall(STGetMatrix(eps->st,0,&H));
647: if (eps->problem_type==EPS_LREP) PetscCall(SlepcCheckMatLREPReduced(H,&lrepreduced));
648: if (eps->problem_type==EPS_BSE || !lrepreduced) { /* change sign of bottom part of the vector */
649: PetscCall(MatNestGetISs(H,is,NULL));
650: if (Wr) {
651: PetscCall(VecGetSubVector(Wr,is[1],&v1));
652: PetscCall(VecScale(v1,-1.0));
653: PetscCall(VecRestoreSubVector(Wr,is[1],&v1));
654: }
655: #if !defined(PETSC_USE_COMPLEX)
656: if (Wi) {
657: PetscCall(VecGetSubVector(Wi,is[1],&v1));
658: PetscCall(VecScale(v1,-1.0));
659: PetscCall(VecRestoreSubVector(Wi,is[1],&v1));
660: }
661: #endif
662: } else if (eps->problem_type==EPS_LREP) { /* swap the two parts of the vector */
663: if (Wr) {
664: PetscCall(VecDuplicate(Wr,&z));
665: PetscCall(VecCopy(Wr,z));
666: PetscCall(MatNestGetISs(H,is,NULL));
667: PetscCall(VecGetSubVector(Wr,is[0],&v0));
668: PetscCall(VecGetSubVector(Wr,is[1],&v1));
669: PetscCall(VecGetSubVector(z,is[0],&z0));
670: PetscCall(VecGetSubVector(z,is[1],&z1));
671: PetscCall(VecCopy(z0,v1));
672: PetscCall(VecCopy(z1,v0));
673: PetscCall(VecRestoreSubVector(Wr,is[0],&v0));
674: PetscCall(VecRestoreSubVector(Wr,is[1],&v1));
675: PetscCall(VecRestoreSubVector(z,is[0],&z0));
676: PetscCall(VecRestoreSubVector(z,is[1],&z1));
677: PetscCall(VecDestroy(&z));
678: }
679: }
680: } else {
681: PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
682: }
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: /*@
687: EPSGetErrorEstimate - Returns the error estimate associated to the `i`-th
688: computed eigenpair.
690: Not Collective
692: Input Parameters:
693: + eps - the linear eigensolver context
694: - i - index of eigenpair
696: Output Parameter:
697: . errest - the error estimate
699: Note:
700: This is the error estimate used internally by the eigensolver. The actual
701: error bound can be computed with `EPSComputeError()`. See discussion at
702: section [](#sec:errbnd).
704: Level: advanced
706: .seealso: [](ch:eps), [](#sec:errbnd), `EPSComputeError()`
707: @*/
708: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
709: {
710: PetscInt nconv;
712: PetscFunctionBegin;
714: PetscAssertPointer(errest,3);
715: EPSCheckSolved(eps,1);
716: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
717: PetscCall(EPS_GetActualConverged(eps,&nconv));
718: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
719: if (nconv==eps->nconv) {
720: *errest = eps->errest[eps->perm[i]];
721: } else {
722: PetscCheck(eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP || eps->problem_type==EPS_HAMILT,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Wrong problem type");
723: /* even index is +lambda, odd index is -lambda, assume both have same error */
724: *errest = eps->errest[eps->perm[i/2]];
725: }
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /*
730: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
731: associated with an eigenpair.
733: Input Parameters:
734: trans - whether A' must be used instead of A
735: kr,ki - eigenvalue
736: xr,xi - eigenvector
737: z - three work vectors (the second one not referenced in complex scalars)
738: */
739: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
740: {
741: PetscInt nmat;
742: Mat A,B;
743: Vec u,w;
744: PetscScalar alpha;
745: #if !defined(PETSC_USE_COMPLEX)
746: Vec v;
747: PetscReal ni,nr;
748: #endif
749: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
751: PetscFunctionBegin;
752: u = z[0]; w = z[2];
753: PetscCall(STGetNumMatrices(eps->st,&nmat));
754: PetscCall(STGetMatrix(eps->st,0,&A));
755: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
757: #if !defined(PETSC_USE_COMPLEX)
758: v = z[1];
759: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
760: #endif
761: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
762: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
763: if (nmat>1) PetscCall((*matmult)(B,xr,w));
764: else PetscCall(VecCopy(xr,w)); /* w=B*x */
765: alpha = trans? -PetscConj(kr): -kr;
766: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
767: }
768: PetscCall(VecNorm(u,NORM_2,norm));
769: #if !defined(PETSC_USE_COMPLEX)
770: } else {
771: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
772: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
773: if (nmat>1) PetscCall((*matmult)(B,xr,v));
774: else PetscCall(VecCopy(xr,v)); /* v=B*xr */
775: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
776: if (nmat>1) PetscCall((*matmult)(B,xi,w));
777: else PetscCall(VecCopy(xi,w)); /* w=B*xi */
778: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
779: }
780: PetscCall(VecNorm(u,NORM_2,&nr));
781: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
782: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
783: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
784: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
785: }
786: PetscCall(VecNorm(u,NORM_2,&ni));
787: *norm = SlepcAbsEigenvalue(nr,ni);
788: }
789: #endif
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: EPSComputeError - Computes the error (based on the residual norm) associated
795: with the `i`-th computed eigenpair.
797: Collective
799: Input Parameters:
800: + eps - the linear eigensolver context
801: . i - the solution index
802: - type - the type of error to compute, see `EPSErrorType`
804: Output Parameter:
805: . error - the error
807: Notes:
808: The error can be computed in various ways, all of them based on the residual
809: norm $\|Ax-\lambda Bx\|_2$ where $(\lambda,x)$ is the approximate eigenpair.
811: If the computation of left eigenvectors was enabled with `EPSSetTwoSided()`,
812: then the error will be computed using the maximum of the value above and
813: the left residual norm $\|y^*A-\lambda y^*B\|_2$, where $y$ is the approximate left
814: eigenvector.
816: Level: beginner
818: .seealso: [](ch:eps), `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`, `EPSSetTwoSided()`
819: @*/
820: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
821: {
822: Mat A,B;
823: Vec xr,xi,w[3];
824: PetscReal t,vecnorm=1.0,errorl;
825: PetscScalar kr,ki;
826: PetscBool flg;
828: PetscFunctionBegin;
832: PetscAssertPointer(error,4);
833: EPSCheckSolved(eps,1);
835: /* allocate work vectors */
836: #if defined(PETSC_USE_COMPLEX)
837: PetscCall(EPSSetWorkVecs(eps,3));
838: xi = NULL;
839: w[1] = NULL;
840: #else
841: PetscCall(EPSSetWorkVecs(eps,5));
842: xi = eps->work[3];
843: w[1] = eps->work[4];
844: #endif
845: xr = eps->work[0];
846: w[0] = eps->work[1];
847: w[2] = eps->work[2];
849: /* compute residual norm */
850: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
851: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
853: /* compute 2-norm of eigenvector */
854: if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
856: /* if two-sided, compute left residual norm and take the maximum */
857: if (eps->twosided) {
858: PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
859: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
860: *error = PetscMax(*error,errorl);
861: }
863: /* compute error */
864: switch (type) {
865: case EPS_ERROR_ABSOLUTE:
866: break;
867: case EPS_ERROR_RELATIVE:
868: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
869: break;
870: case EPS_ERROR_BACKWARD:
871: /* initialization of matrix norms */
872: if (!eps->nrma) {
873: PetscCall(STGetMatrix(eps->st,0,&A));
874: PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
875: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
876: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
877: }
878: if (eps->isgeneralized) {
879: if (!eps->nrmb) {
880: PetscCall(STGetMatrix(eps->st,1,&B));
881: PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
882: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
883: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
884: }
885: } else eps->nrmb = 1.0;
886: t = SlepcAbsEigenvalue(kr,ki);
887: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
888: break;
889: default:
890: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
891: }
892: PetscFunctionReturn(PETSC_SUCCESS);
893: }
895: /*
896: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
897: for the recurrence that builds the right subspace.
899: Collective
901: Input Parameters:
902: + eps - the linear eigensolver context
903: - i - iteration number
905: Output Parameter:
906: . breakdown - flag indicating that a breakdown has occurred
908: Notes:
909: The start vector is computed from another vector: for the first step (i=0),
910: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
911: vector is created. Then this vector is forced to be in the range of OP (only
912: for generalized definite problems) and orthonormalized with respect to all
913: V-vectors up to i-1. The resulting vector is placed in V[i].
915: The flag breakdown is set to true if either i=0 and the vector belongs to the
916: deflation space, or i>0 and the vector is linearly dependent with respect
917: to the V-vectors.
918: */
919: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
920: {
921: PetscReal norm;
922: PetscBool lindep;
923: Vec w,z;
925: PetscFunctionBegin;
929: /* For the first step, use the first initial vector, otherwise a random one */
930: if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
932: /* Force the vector to be in the range of OP for generalized problems with B-inner product */
933: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
934: PetscCall(BVCreateVec(eps->V,&w));
935: PetscCall(BVCopyVec(eps->V,i,w));
936: PetscCall(BVGetColumn(eps->V,i,&z));
937: PetscCall(STApply(eps->st,w,z));
938: PetscCall(BVRestoreColumn(eps->V,i,&z));
939: PetscCall(VecDestroy(&w));
940: }
942: /* Orthonormalize the vector with respect to previous vectors */
943: PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
944: if (breakdown) *breakdown = lindep;
945: else if (lindep || norm == 0.0) {
946: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
947: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
948: }
949: PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: /*
954: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
955: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
956: */
957: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
958: {
959: PetscReal norm;
960: PetscBool lindep;
961: Vec w,z;
963: PetscFunctionBegin;
967: /* For the first step, use the first initial vector, otherwise a random one */
968: if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
970: /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
971: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
972: PetscCall(BVCreateVec(eps->W,&w));
973: PetscCall(BVCopyVec(eps->W,i,w));
974: PetscCall(BVGetColumn(eps->W,i,&z));
975: PetscCall(STApplyHermitianTranspose(eps->st,w,z));
976: PetscCall(BVRestoreColumn(eps->W,i,&z));
977: PetscCall(VecDestroy(&w));
978: }
980: /* Orthonormalize the vector with respect to previous vectors */
981: PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
982: if (breakdown) *breakdown = lindep;
983: else if (lindep || norm == 0.0) {
984: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
985: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
986: }
987: PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
988: PetscFunctionReturn(PETSC_SUCCESS);
989: }