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.
123: Level: beginner
125: .seealso: [](ch:eps), `EPSCreate()`, `EPSSetUp()`, `EPSDestroy()`, `EPSSetOperators()`, `EPSGetConverged()`, `EPSGetConvergedReason()`
126: @*/
127: PetscErrorCode EPSSolve(EPS eps)
128: {
129: PetscInt i;
130: PetscBool hasname;
131: STMatMode matmode;
132: Mat A,B;
134: PetscFunctionBegin;
136: if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
137: PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
139: /* Call setup */
140: PetscCall(EPSSetUp(eps));
142: /* Safeguard for matrices of size 0 */
143: if (eps->n == 0) {
144: eps->nconv = 0;
145: eps->reason = EPS_CONVERGED_TOL;
146: eps->state = EPS_STATE_SOLVED;
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: eps->nconv = 0;
151: eps->its = 0;
152: for (i=0;i<eps->ncv;i++) {
153: eps->eigr[i] = 0.0;
154: eps->eigi[i] = 0.0;
155: eps->errest[i] = 0.0;
156: eps->perm[i] = i;
157: }
158: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
159: PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
161: /* Call solver */
162: PetscUseTypeMethod(eps,solve);
163: PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
164: eps->state = EPS_STATE_SOLVED;
166: /* Only the first nconv columns contain useful information (except in CISS) */
167: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
168: if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
170: /* If inplace, purify eigenvectors before reverting operator */
171: PetscCall(STGetMatMode(eps->st,&matmode));
172: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
173: PetscCall(STPostSolve(eps->st));
175: /* Map eigenvalues back to the original problem if appropriate */
176: PetscCall(EPSComputeValues(eps));
178: #if !defined(PETSC_USE_COMPLEX)
179: /* Reorder conjugate eigenvalues (positive imaginary first) */
180: for (i=0;i<eps->nconv-1;i++) {
181: if (eps->eigi[i] != 0 && (eps->problem_type!=EPS_HAMILT || eps->eigr[i]!=0)) {
182: /* conjugate eigenvalues */
183: if (eps->eigi[i] < 0) {
184: eps->eigi[i] = -eps->eigi[i];
185: eps->eigi[i+1] = -eps->eigi[i+1];
186: /* the next correction only works with eigenvectors */
187: PetscCall(EPSComputeVectors(eps));
188: PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
189: if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
190: }
191: i++;
192: }
193: }
194: #endif
196: /* Sort eigenvalues according to eps->which parameter */
197: if (eps->problem_type==EPS_HAMILT) PetscCall(SlepcSortEigenvaluesSpecial(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
198: else PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
199: PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
201: /* Various viewers */
202: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
203: PetscCall(EPSConvergedReasonViewFromOptions(eps));
204: PetscCall(EPSErrorViewFromOptions(eps));
205: PetscCall(EPSValuesViewFromOptions(eps));
206: PetscCall(EPSVectorsViewFromOptions(eps));
208: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
209: if (hasname) {
210: PetscCall(EPSGetOperators(eps,&A,NULL));
211: PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
212: }
213: if (eps->isgeneralized) {
214: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
215: if (hasname) {
216: PetscCall(EPSGetOperators(eps,NULL,&B));
217: PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
218: }
219: }
221: /* Remove deflation and initial subspaces */
222: if (eps->nds) {
223: PetscCall(BVSetNumConstraints(eps->V,0));
224: eps->nds = 0;
225: }
226: eps->nini = 0;
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: EPSGetIterationNumber - Gets the current iteration number. If the
232: call to `EPSSolve()` is complete, then it returns the number of iterations
233: carried out by the solution method.
235: Not Collective
237: Input Parameter:
238: . eps - the linear eigensolver context
240: Output Parameter:
241: . its - number of iterations
243: Note:
244: During the $i$-th iteration this call returns $i-1$. If `EPSSolve()` is
245: complete, then parameter `its` contains either the iteration number at
246: which convergence was successfully reached, or failure was detected.
247: Call `EPSGetConvergedReason()` to determine if the solver converged or
248: failed and why.
250: Level: intermediate
252: .seealso: [](ch:eps), `EPSGetConvergedReason()`, `EPSSetTolerances()`
253: @*/
254: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
255: {
256: PetscFunctionBegin;
258: PetscAssertPointer(its,2);
259: *its = eps->its;
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: /*@
264: EPSGetConverged - Gets the number of converged eigenpairs.
266: Not Collective
268: Input Parameter:
269: . eps - the linear eigensolver context
271: Output Parameter:
272: . nconv - number of converged eigenpairs
274: Notes:
275: This function should be called after `EPSSolve()` has finished.
277: The value `nconv` may be different from the number of requested solutions
278: `nev`, but not larger than `ncv`, see `EPSSetDimensions()`.
280: Level: beginner
282: .seealso: [](ch:eps), `EPSSetDimensions()`, `EPSSolve()`, `EPSGetEigenpair()`
283: @*/
284: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
285: {
286: PetscFunctionBegin;
288: PetscAssertPointer(nconv,2);
289: EPSCheckSolved(eps,1);
290: PetscCall(EPS_GetActualConverged(eps,nconv));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: EPSGetConvergedReason - Gets the reason why the `EPSSolve()` iteration was
296: stopped.
298: Not Collective
300: Input Parameter:
301: . eps - the linear eigensolver context
303: Output Parameter:
304: . reason - negative value indicates diverged, positive value converged, see
305: `EPSConvergedReason` for the possible values
307: Options Database Key:
308: . -eps_converged_reason - print reason for convergence/divergence, and number of iterations
310: Note:
311: If this routine is called before or doing the `EPSSolve()` the value of
312: `EPS_CONVERGED_ITERATING` is returned.
314: Level: intermediate
316: .seealso: [](ch:eps), `EPSSetTolerances()`, `EPSSolve()`, `EPSConvergedReason`
317: @*/
318: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
319: {
320: PetscFunctionBegin;
322: PetscAssertPointer(reason,2);
323: EPSCheckSolved(eps,1);
324: *reason = eps->reason;
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*@
329: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
330: subspace.
332: Collective
334: Input Parameter:
335: . eps - the linear eigensolver context
337: Output Parameter:
338: . v - an array of vectors
340: Notes:
341: This function should be called after `EPSSolve()` has finished.
343: The user should provide in `v` an array of `nconv` vectors, where `nconv`
344: is the value returned by `EPSGetConverged()`.
346: The first $k$ vectors returned in `v` span an invariant subspace associated
347: with the first $k$ computed eigenvalues (note that this is not true if the
348: $k$-th eigenvalue is complex and matrix $A$ is real; in this case the first
349: $k+1$ vectors should be used). An invariant subspace $X$ of $A$ satisfies
350: $Ax\in X, \forall x\in X$ (a similar definition applies for generalized
351: eigenproblems).
353: Level: intermediate
355: .seealso: [](ch:eps), `EPSGetEigenpair()`, `EPSGetConverged()`, `EPSSolve()`
356: @*/
357: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
358: {
359: PetscInt i;
360: BV V=eps->V;
361: Vec w;
363: PetscFunctionBegin;
365: PetscAssertPointer(v,2);
367: EPSCheckSolved(eps,1);
368: PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
369: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
370: PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
371: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
372: PetscCall(BVCopy(eps->V,V));
373: for (i=0;i<eps->nconv;i++) {
374: PetscCall(BVGetColumn(V,i,&w));
375: PetscCall(VecPointwiseDivide(w,w,eps->D));
376: PetscCall(BVRestoreColumn(V,i,&w));
377: }
378: PetscCall(BVOrthogonalize(V,NULL));
379: }
380: for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
381: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: EPSGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
387: `EPSSolve()`. The solution consists in both the eigenvalue and the eigenvector.
389: Collective
391: Input Parameters:
392: + eps - the linear eigensolver context
393: - i - index of the solution
395: Output Parameters:
396: + eigr - real part of eigenvalue
397: . eigi - imaginary part of eigenvalue
398: . Vr - real part of eigenvector
399: - Vi - imaginary part of eigenvector
401: Notes:
402: It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
403: required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
404: they must be created by the calling program with e.g. `MatCreateVecs()`.
406: If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
407: configured with complex scalars the eigenvalue is stored
408: directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi` is
409: set to zero). In both cases, the user can pass `NULL` in `eigi` and `Vi`.
411: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
412: Eigenpairs are indexed according to the ordering criterion established
413: with `EPSSetWhichEigenpairs()`.
415: The 2-norm of the eigenvector is one unless the problem is generalized
416: Hermitian. In this case the eigenvector is normalized with respect to the
417: norm defined by the $B$ matrix.
419: In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
420: [](#sec:structured-vectors).
422: Level: beginner
424: .seealso: [](ch:eps), `EPSGetEigenvalue()`, `EPSGetEigenvector()`, `EPSGetLeftEigenvector()`, `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetInvariantSubspace()`
425: @*/
426: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
427: {
428: PetscInt nconv;
430: PetscFunctionBegin;
433: EPSCheckSolved(eps,1);
434: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
435: PetscCall(EPS_GetActualConverged(eps,&nconv));
436: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
437: PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
438: if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: EPSGetEigenvalue - Gets the `i`-th eigenvalue as computed by `EPSSolve()`.
445: Not Collective
447: Input Parameters:
448: + eps - the linear eigensolver context
449: - i - index of the solution
451: Output Parameters:
452: + eigr - real part of eigenvalue
453: - eigi - imaginary part of eigenvalue
455: Notes:
456: If the eigenvalue is real, then `eigi` is set to zero. If PETSc is
457: configured with complex scalars the eigenvalue is stored
458: directly in `eigr` (`eigi` is set to zero).
460: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
461: Eigenpairs are indexed according to the ordering criterion established
462: with `EPSSetWhichEigenpairs()`.
464: Level: beginner
466: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`
467: @*/
468: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
469: {
470: PetscInt k,nconv;
471: #if !defined(PETSC_USE_COMPLEX)
472: PetscInt k2, iquad;
473: #endif
475: PetscFunctionBegin;
477: EPSCheckSolved(eps,1);
478: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
479: PetscCall(EPS_GetActualConverged(eps,&nconv));
480: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
481: if (nconv==eps->nconv) {
482: k = eps->perm[i];
483: #if defined(PETSC_USE_COMPLEX)
484: if (eigr) *eigr = eps->eigr[k];
485: if (eigi) *eigi = 0;
486: #else
487: if (eigr) *eigr = eps->eigr[k];
488: if (eigi) *eigi = eps->eigi[k];
489: #endif
490: } else {
491: 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");
492: if (eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP) {
493: /* BSE problem, even index is +lambda, odd index is -lambda */
494: k = eps->perm[i/2];
495: #if defined(PETSC_USE_COMPLEX)
496: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
497: if (eigi) *eigi = 0;
498: #else
499: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
500: if (eigi) *eigi = eps->eigi[k];
501: #endif
502: } else if (eps->problem_type==EPS_HAMILT) {
503: /* Hamiltonian eigenproblem */
504: k = eps->perm[i/2];
505: #if defined(PETSC_USE_COMPLEX)
506: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
507: if (eigi) *eigi = 0;
508: #else
509: if (eps->eigi[k]==0.0) { /* real eigenvalue */
510: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
511: if (eigi) *eigi = 0.0;
512: } else if (eps->eigr[k]==0.0) { /* purely imaginary eigenvalue */
513: if (eigr) *eigr = 0.0;
514: if (eigi) *eigi = (i%2)? -eps->eigi[k]: eps->eigi[k];
515: } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
516: iquad = i%2; /* index within the 4 values */
517: if (i>1) {
518: k2 = eps->perm[(i-2)/2];
519: if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
520: }
521: if (eigr) *eigr = (iquad<2)? -eps->eigr[k]: eps->eigr[k];
522: if (eigi) *eigi = (iquad%3)? -eps->eigi[k]: eps->eigi[k];
523: }
524: #endif
525: }
526: }
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*@
531: EPSGetEigenvector - Gets the `i`-th right eigenvector as computed by `EPSSolve()`.
533: Collective
535: Input Parameters:
536: + eps - the linear eigensolver context
537: - i - index of the solution
539: Output Parameters:
540: + Vr - real part of eigenvector
541: - Vi - imaginary part of eigenvector
543: Notes:
544: The caller must provide valid `Vec` objects, i.e., they must be created
545: by the calling program with e.g. `MatCreateVecs()`.
547: If the corresponding eigenvalue is real, then `Vi` is set to zero. If PETSc is
548: configured with complex scalars the eigenvector is stored
549: directly in `Vr` (`Vi` is set to zero). In any case, the user can pass `NULL` in `Vr`
550: or `Vi` if one of them is not required.
552: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
553: Eigenpairs are indexed according to the ordering criterion established
554: with `EPSSetWhichEigenpairs()`.
556: The 2-norm of the eigenvector is one unless the problem is generalized
557: Hermitian. In this case the eigenvector is normalized with respect to the
558: norm defined by the $B$ matrix.
560: In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
561: [](#sec:structured-vectors).
563: Level: beginner
565: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`, `EPSGetLeftEigenvector()`
566: @*/
567: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
568: {
569: PetscInt nconv;
571: PetscFunctionBegin;
576: EPSCheckSolved(eps,1);
577: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
578: PetscCall(EPS_GetActualConverged(eps,&nconv));
579: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
580: PetscCall(EPSComputeVectors(eps));
581: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: /*@
586: EPSGetLeftEigenvector - Gets the `i`-th left eigenvector as computed by `EPSSolve()`.
588: Collective
590: Input Parameters:
591: + eps - the linear eigensolver context
592: - i - index of the solution
594: Output Parameters:
595: + Wr - real part of left eigenvector
596: - Wi - imaginary part of left eigenvector
598: Notes:
599: The caller must provide valid `Vec` objects, i.e., they must be created
600: by the calling program with e.g. `MatCreateVecs()`.
602: If the corresponding eigenvalue is real, then `Wi` is set to zero. If PETSc is
603: configured with complex scalars the eigenvector is stored directly in `Wr`
604: (`Wi` is set to zero). In any case, the user can pass `NULL` in `Wr` or `Wi` if
605: one of them is not required.
607: The index `i` should be a value between 0 and `nconv`-1 (see `EPSGetConverged()`).
608: Eigensolutions are indexed according to the ordering criterion established
609: with `EPSSetWhichEigenpairs()`.
611: Left eigenvectors are available only if the `twosided` flag was set, see
612: `EPSSetTwoSided()`.
614: In case of structured eigenproblems such as `EPS_BSE`, see the discussion about
615: [](#sec:structured-vectors).
617: Level: intermediate
619: .seealso: [](ch:eps), `EPSGetEigenvector()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSSetTwoSided()`
620: @*/
621: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
622: {
623: PetscInt nconv;
624: PetscBool trivial,lrepreduced=PETSC_FALSE;
625: Mat H;
626: IS is[2];
627: Vec v0,v1,z,z0,z1;
629: PetscFunctionBegin;
634: EPSCheckSolved(eps,1);
635: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
636: PetscCall(EPS_GetActualConverged(eps,&nconv));
637: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
639: 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;
640: if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
642: PetscCall(EPSComputeVectors(eps));
643: if (trivial) {
644: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
645: PetscCall(STGetMatrix(eps->st,0,&H));
646: if (eps->problem_type==EPS_LREP) PetscCall(SlepcCheckMatLREPReduced(H,&lrepreduced));
647: if (eps->problem_type==EPS_BSE || !lrepreduced) { /* change sign of bottom part of the vector */
648: PetscCall(MatNestGetISs(H,is,NULL));
649: if (Wr) {
650: PetscCall(VecGetSubVector(Wr,is[1],&v1));
651: PetscCall(VecScale(v1,-1.0));
652: PetscCall(VecRestoreSubVector(Wr,is[1],&v1));
653: }
654: #if !defined(PETSC_USE_COMPLEX)
655: if (Wi) {
656: PetscCall(VecGetSubVector(Wi,is[1],&v1));
657: PetscCall(VecScale(v1,-1.0));
658: PetscCall(VecRestoreSubVector(Wi,is[1],&v1));
659: }
660: #endif
661: } else if (eps->problem_type==EPS_LREP) { /* swap the two parts of the vector */
662: if (Wr) {
663: PetscCall(VecDuplicate(Wr,&z));
664: PetscCall(VecCopy(Wr,z));
665: PetscCall(MatNestGetISs(H,is,NULL));
666: PetscCall(VecGetSubVector(Wr,is[0],&v0));
667: PetscCall(VecGetSubVector(Wr,is[1],&v1));
668: PetscCall(VecGetSubVector(z,is[0],&z0));
669: PetscCall(VecGetSubVector(z,is[1],&z1));
670: PetscCall(VecCopy(z0,v1));
671: PetscCall(VecCopy(z1,v0));
672: PetscCall(VecRestoreSubVector(Wr,is[0],&v0));
673: PetscCall(VecRestoreSubVector(Wr,is[1],&v1));
674: PetscCall(VecRestoreSubVector(z,is[0],&z0));
675: PetscCall(VecRestoreSubVector(z,is[1],&z1));
676: PetscCall(VecDestroy(&z));
677: }
678: }
679: } else {
680: PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
681: }
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: /*@
686: EPSGetErrorEstimate - Returns the error estimate associated to the `i`-th
687: computed eigenpair.
689: Not Collective
691: Input Parameters:
692: + eps - the linear eigensolver context
693: - i - index of eigenpair
695: Output Parameter:
696: . errest - the error estimate
698: Note:
699: This is the error estimate used internally by the eigensolver. The actual
700: error bound can be computed with `EPSComputeError()`. See discussion at
701: section [](#sec:errbnd).
703: Level: advanced
705: .seealso: [](ch:eps), [](#sec:errbnd), `EPSComputeError()`
706: @*/
707: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
708: {
709: PetscInt nconv;
711: PetscFunctionBegin;
713: PetscAssertPointer(errest,3);
714: EPSCheckSolved(eps,1);
715: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
716: PetscCall(EPS_GetActualConverged(eps,&nconv));
717: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
718: if (nconv==eps->nconv) {
719: *errest = eps->errest[eps->perm[i]];
720: } else {
721: PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
722: /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
723: *errest = eps->errest[eps->perm[i/2]];
724: }
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /*
729: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
730: associated with an eigenpair.
732: Input Parameters:
733: trans - whether A' must be used instead of A
734: kr,ki - eigenvalue
735: xr,xi - eigenvector
736: z - three work vectors (the second one not referenced in complex scalars)
737: */
738: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
739: {
740: PetscInt nmat;
741: Mat A,B;
742: Vec u,w;
743: PetscScalar alpha;
744: #if !defined(PETSC_USE_COMPLEX)
745: Vec v;
746: PetscReal ni,nr;
747: #endif
748: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
750: PetscFunctionBegin;
751: u = z[0]; w = z[2];
752: PetscCall(STGetNumMatrices(eps->st,&nmat));
753: PetscCall(STGetMatrix(eps->st,0,&A));
754: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
756: #if !defined(PETSC_USE_COMPLEX)
757: v = z[1];
758: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
759: #endif
760: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
761: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
762: if (nmat>1) PetscCall((*matmult)(B,xr,w));
763: else PetscCall(VecCopy(xr,w)); /* w=B*x */
764: alpha = trans? -PetscConj(kr): -kr;
765: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
766: }
767: PetscCall(VecNorm(u,NORM_2,norm));
768: #if !defined(PETSC_USE_COMPLEX)
769: } else {
770: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
771: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
772: if (nmat>1) PetscCall((*matmult)(B,xr,v));
773: else PetscCall(VecCopy(xr,v)); /* v=B*xr */
774: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
775: if (nmat>1) PetscCall((*matmult)(B,xi,w));
776: else PetscCall(VecCopy(xi,w)); /* w=B*xi */
777: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
778: }
779: PetscCall(VecNorm(u,NORM_2,&nr));
780: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
781: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
782: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
783: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
784: }
785: PetscCall(VecNorm(u,NORM_2,&ni));
786: *norm = SlepcAbsEigenvalue(nr,ni);
787: }
788: #endif
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: /*@
793: EPSComputeError - Computes the error (based on the residual norm) associated
794: with the `i`-th computed eigenpair.
796: Collective
798: Input Parameters:
799: + eps - the linear eigensolver context
800: . i - the solution index
801: - type - the type of error to compute, see `EPSErrorType`
803: Output Parameter:
804: . error - the error
806: Notes:
807: The error can be computed in various ways, all of them based on the residual
808: norm $\|Ax-\lambda Bx\|_2$ where $(\lambda,x)$ is the approximate eigenpair.
810: If the computation of left eigenvectors was enabled with `EPSSetTwoSided()`,
811: then the error will be computed using the maximum of the value above and
812: the left residual norm $\|y^*A-\lambda y^*B\|_2$, where $y$ is the approximate left
813: eigenvector.
815: Level: beginner
817: .seealso: [](ch:eps), `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`, `EPSSetTwoSided()`
818: @*/
819: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
820: {
821: Mat A,B;
822: Vec xr,xi,w[3];
823: PetscReal t,vecnorm=1.0,errorl;
824: PetscScalar kr,ki;
825: PetscBool flg;
827: PetscFunctionBegin;
831: PetscAssertPointer(error,4);
832: EPSCheckSolved(eps,1);
834: /* allocate work vectors */
835: #if defined(PETSC_USE_COMPLEX)
836: PetscCall(EPSSetWorkVecs(eps,3));
837: xi = NULL;
838: w[1] = NULL;
839: #else
840: PetscCall(EPSSetWorkVecs(eps,5));
841: xi = eps->work[3];
842: w[1] = eps->work[4];
843: #endif
844: xr = eps->work[0];
845: w[0] = eps->work[1];
846: w[2] = eps->work[2];
848: /* compute residual norm */
849: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
850: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
852: /* compute 2-norm of eigenvector */
853: if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
855: /* if two-sided, compute left residual norm and take the maximum */
856: if (eps->twosided) {
857: PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
858: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
859: *error = PetscMax(*error,errorl);
860: }
862: /* compute error */
863: switch (type) {
864: case EPS_ERROR_ABSOLUTE:
865: break;
866: case EPS_ERROR_RELATIVE:
867: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
868: break;
869: case EPS_ERROR_BACKWARD:
870: /* initialization of matrix norms */
871: if (!eps->nrma) {
872: PetscCall(STGetMatrix(eps->st,0,&A));
873: PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
874: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
875: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
876: }
877: if (eps->isgeneralized) {
878: if (!eps->nrmb) {
879: PetscCall(STGetMatrix(eps->st,1,&B));
880: PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
881: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
882: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
883: }
884: } else eps->nrmb = 1.0;
885: t = SlepcAbsEigenvalue(kr,ki);
886: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
887: break;
888: default:
889: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
890: }
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: /*
895: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
896: for the recurrence that builds the right subspace.
898: Collective
900: Input Parameters:
901: + eps - the linear eigensolver context
902: - i - iteration number
904: Output Parameter:
905: . breakdown - flag indicating that a breakdown has occurred
907: Notes:
908: The start vector is computed from another vector: for the first step (i=0),
909: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
910: vector is created. Then this vector is forced to be in the range of OP (only
911: for generalized definite problems) and orthonormalized with respect to all
912: V-vectors up to i-1. The resulting vector is placed in V[i].
914: The flag breakdown is set to true if either i=0 and the vector belongs to the
915: deflation space, or i>0 and the vector is linearly dependent with respect
916: to the V-vectors.
917: */
918: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
919: {
920: PetscReal norm;
921: PetscBool lindep;
922: Vec w,z;
924: PetscFunctionBegin;
928: /* For the first step, use the first initial vector, otherwise a random one */
929: if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
931: /* Force the vector to be in the range of OP for generalized problems with B-inner product */
932: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
933: PetscCall(BVCreateVec(eps->V,&w));
934: PetscCall(BVCopyVec(eps->V,i,w));
935: PetscCall(BVGetColumn(eps->V,i,&z));
936: PetscCall(STApply(eps->st,w,z));
937: PetscCall(BVRestoreColumn(eps->V,i,&z));
938: PetscCall(VecDestroy(&w));
939: }
941: /* Orthonormalize the vector with respect to previous vectors */
942: PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
943: if (breakdown) *breakdown = lindep;
944: else if (lindep || norm == 0.0) {
945: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
946: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
947: }
948: PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
949: PetscFunctionReturn(PETSC_SUCCESS);
950: }
952: /*
953: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
954: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
955: */
956: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
957: {
958: PetscReal norm;
959: PetscBool lindep;
960: Vec w,z;
962: PetscFunctionBegin;
966: /* For the first step, use the first initial vector, otherwise a random one */
967: if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
969: /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
970: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
971: PetscCall(BVCreateVec(eps->W,&w));
972: PetscCall(BVCopyVec(eps->W,i,w));
973: PetscCall(BVGetColumn(eps->W,i,&z));
974: PetscCall(STApplyHermitianTranspose(eps->st,w,z));
975: PetscCall(BVRestoreColumn(eps->W,i,&z));
976: PetscCall(VecDestroy(&w));
977: }
979: /* Orthonormalize the vector with respect to previous vectors */
980: PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
981: if (breakdown) *breakdown = lindep;
982: else if (lindep || norm == 0.0) {
983: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
984: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
985: }
986: PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }