Actual source code: epssolve.c
slepc-3.23.2 2025-06-30
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 - eigensolver context obtained from EPSCreate()
97: Options Database Keys:
98: + -eps_view - print information about the solver used
99: . -eps_view_mat0 - view the first matrix (A)
100: . -eps_view_mat1 - view the second matrix (B)
101: . -eps_view_vectors - view the computed eigenvectors
102: . -eps_view_values - view the computed eigenvalues
103: . -eps_converged_reason - print reason for convergence, and number of iterations
104: . -eps_error_absolute - print absolute errors of each eigenpair
105: . -eps_error_relative - print relative errors of each eigenpair
106: - -eps_error_backward - print backward errors of each eigenpair
108: Notes:
109: All the command-line options listed above admit an optional argument specifying
110: the viewer type and options. For instance, use '-eps_view_mat0 binary:amatrix.bin'
111: to save the A matrix to a binary file, '-eps_view_values draw' to draw the computed
112: eigenvalues graphically, or '-eps_error_relative :myerr.m:ascii_matlab' to save
113: the errors in a file that can be executed in Matlab.
115: Level: beginner
117: .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
118: @*/
119: PetscErrorCode EPSSolve(EPS eps)
120: {
121: PetscInt i;
122: PetscBool hasname;
123: STMatMode matmode;
124: Mat A,B;
126: PetscFunctionBegin;
128: if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
129: PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
131: /* Call setup */
132: PetscCall(EPSSetUp(eps));
133: eps->nconv = 0;
134: eps->its = 0;
135: for (i=0;i<eps->ncv;i++) {
136: eps->eigr[i] = 0.0;
137: eps->eigi[i] = 0.0;
138: eps->errest[i] = 0.0;
139: eps->perm[i] = i;
140: }
141: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
142: PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
144: /* Call solver */
145: PetscUseTypeMethod(eps,solve);
146: PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
147: eps->state = EPS_STATE_SOLVED;
149: /* Only the first nconv columns contain useful information (except in CISS) */
150: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
151: if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
153: /* If inplace, purify eigenvectors before reverting operator */
154: PetscCall(STGetMatMode(eps->st,&matmode));
155: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
156: PetscCall(STPostSolve(eps->st));
158: /* Map eigenvalues back to the original problem if appropriate */
159: PetscCall(EPSComputeValues(eps));
161: #if !defined(PETSC_USE_COMPLEX)
162: /* Reorder conjugate eigenvalues (positive imaginary first) */
163: for (i=0;i<eps->nconv-1;i++) {
164: if (eps->eigi[i] != 0) {
165: if (eps->eigi[i] < 0) {
166: eps->eigi[i] = -eps->eigi[i];
167: eps->eigi[i+1] = -eps->eigi[i+1];
168: /* the next correction only works with eigenvectors */
169: PetscCall(EPSComputeVectors(eps));
170: PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
171: if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
172: }
173: i++;
174: }
175: }
176: #endif
178: /* Sort eigenvalues according to eps->which parameter */
179: PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
180: PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
182: /* Various viewers */
183: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
184: PetscCall(EPSConvergedReasonViewFromOptions(eps));
185: PetscCall(EPSErrorViewFromOptions(eps));
186: PetscCall(EPSValuesViewFromOptions(eps));
187: PetscCall(EPSVectorsViewFromOptions(eps));
189: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
190: if (hasname) {
191: PetscCall(EPSGetOperators(eps,&A,NULL));
192: PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
193: }
194: if (eps->isgeneralized) {
195: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
196: if (hasname) {
197: PetscCall(EPSGetOperators(eps,NULL,&B));
198: PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
199: }
200: }
202: /* Remove deflation and initial subspaces */
203: if (eps->nds) {
204: PetscCall(BVSetNumConstraints(eps->V,0));
205: eps->nds = 0;
206: }
207: eps->nini = 0;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /*@
212: EPSGetIterationNumber - Gets the current iteration number. If the
213: call to EPSSolve() is complete, then it returns the number of iterations
214: carried out by the solution method.
216: Not Collective
218: Input Parameter:
219: . eps - the eigensolver context
221: Output Parameter:
222: . its - number of iterations
224: Note:
225: During the i-th iteration this call returns i-1. If EPSSolve() is
226: complete, then parameter "its" contains either the iteration number at
227: which convergence was successfully reached, or failure was detected.
228: Call EPSGetConvergedReason() to determine if the solver converged or
229: failed and why.
231: Level: intermediate
233: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
234: @*/
235: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
236: {
237: PetscFunctionBegin;
239: PetscAssertPointer(its,2);
240: *its = eps->its;
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*@
245: EPSGetConverged - Gets the number of converged eigenpairs.
247: Not Collective
249: Input Parameter:
250: . eps - the eigensolver context
252: Output Parameter:
253: . nconv - number of converged eigenpairs
255: Note:
256: This function should be called after EPSSolve() has finished.
258: Level: beginner
260: .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
261: @*/
262: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
263: {
264: PetscFunctionBegin;
266: PetscAssertPointer(nconv,2);
267: EPSCheckSolved(eps,1);
268: PetscCall(EPS_GetActualConverged(eps,nconv));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: /*@
273: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
274: stopped.
276: Not Collective
278: Input Parameter:
279: . eps - the eigensolver context
281: Output Parameter:
282: . reason - negative value indicates diverged, positive value converged
284: Options Database Key:
285: . -eps_converged_reason - print the reason to a viewer
287: Notes:
288: Possible values for reason are
289: + EPS_CONVERGED_TOL - converged up to tolerance
290: . EPS_CONVERGED_USER - converged due to a user-defined condition
291: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
292: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
293: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
295: Can only be called after the call to EPSSolve() is complete.
297: Level: intermediate
299: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
300: @*/
301: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
302: {
303: PetscFunctionBegin;
305: PetscAssertPointer(reason,2);
306: EPSCheckSolved(eps,1);
307: *reason = eps->reason;
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
313: subspace.
315: Collective
317: Input Parameter:
318: . eps - the eigensolver context
320: Output Parameter:
321: . v - an array of vectors
323: Notes:
324: This function should be called after EPSSolve() has finished.
326: The user should provide in v an array of nconv vectors, where nconv is
327: the value returned by EPSGetConverged().
329: The first k vectors returned in v span an invariant subspace associated
330: with the first k computed eigenvalues (note that this is not true if the
331: k-th eigenvalue is complex and matrix A is real; in this case the first
332: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
333: in X for all x in X (a similar definition applies for generalized
334: eigenproblems).
336: Level: intermediate
338: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
339: @*/
340: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
341: {
342: PetscInt i;
343: BV V=eps->V;
344: Vec w;
346: PetscFunctionBegin;
348: PetscAssertPointer(v,2);
350: EPSCheckSolved(eps,1);
351: PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
352: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
353: PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
354: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
355: PetscCall(BVCopy(eps->V,V));
356: for (i=0;i<eps->nconv;i++) {
357: PetscCall(BVGetColumn(V,i,&w));
358: PetscCall(VecPointwiseDivide(w,w,eps->D));
359: PetscCall(BVRestoreColumn(V,i,&w));
360: }
361: PetscCall(BVOrthogonalize(V,NULL));
362: }
363: for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
364: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*@
369: EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
370: EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
372: Collective
374: Input Parameters:
375: + eps - eigensolver context
376: - i - index of the solution
378: Output Parameters:
379: + eigr - real part of eigenvalue
380: . eigi - imaginary part of eigenvalue
381: . Vr - real part of eigenvector
382: - Vi - imaginary part of eigenvector
384: Notes:
385: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
386: required. Otherwise, the caller must provide valid Vec objects, i.e.,
387: they must be created by the calling program with e.g. MatCreateVecs().
389: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
390: configured with complex scalars the eigenvalue is stored
391: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
392: set to zero). In both cases, the user can pass NULL in eigi and Vi.
394: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
395: Eigenpairs are indexed according to the ordering criterion established
396: with EPSSetWhichEigenpairs().
398: The 2-norm of the eigenvector is one unless the problem is generalized
399: Hermitian. In this case the eigenvector is normalized with respect to the
400: norm defined by the B matrix.
402: Level: beginner
404: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
405: EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
406: @*/
407: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
408: {
409: PetscInt nconv;
411: PetscFunctionBegin;
414: EPSCheckSolved(eps,1);
415: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
416: PetscCall(EPS_GetActualConverged(eps,&nconv));
417: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
418: PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
419: if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /*@
424: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
426: Not Collective
428: Input Parameters:
429: + eps - eigensolver context
430: - i - index of the solution
432: Output Parameters:
433: + eigr - real part of eigenvalue
434: - eigi - imaginary part of eigenvalue
436: Notes:
437: If the eigenvalue is real, then eigi is set to zero. If PETSc is
438: configured with complex scalars the eigenvalue is stored
439: directly in eigr (eigi is set to zero).
441: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
442: Eigenpairs are indexed according to the ordering criterion established
443: with EPSSetWhichEigenpairs().
445: Level: beginner
447: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
448: @*/
449: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
450: {
451: PetscInt k,nconv;
453: PetscFunctionBegin;
455: EPSCheckSolved(eps,1);
456: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
457: PetscCall(EPS_GetActualConverged(eps,&nconv));
458: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
459: if (nconv==eps->nconv) {
460: k = eps->perm[i];
461: #if defined(PETSC_USE_COMPLEX)
462: if (eigr) *eigr = eps->eigr[k];
463: if (eigi) *eigi = 0;
464: #else
465: if (eigr) *eigr = eps->eigr[k];
466: if (eigi) *eigi = eps->eigi[k];
467: #endif
468: } else {
469: PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
470: /* BSE problem, even index is +lambda, odd index is -lambda */
471: k = eps->perm[i/2];
472: #if defined(PETSC_USE_COMPLEX)
473: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
474: if (eigi) *eigi = 0;
475: #else
476: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
477: if (eigi) *eigi = eps->eigi[k];
478: #endif
479: }
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: /*@
484: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
486: Collective
488: Input Parameters:
489: + eps - eigensolver context
490: - i - index of the solution
492: Output Parameters:
493: + Vr - real part of eigenvector
494: - Vi - imaginary part of eigenvector
496: Notes:
497: The caller must provide valid Vec objects, i.e., they must be created
498: by the calling program with e.g. MatCreateVecs().
500: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
501: configured with complex scalars the eigenvector is stored
502: directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
503: or Vi if one of them is not required.
505: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
506: Eigenpairs are indexed according to the ordering criterion established
507: with EPSSetWhichEigenpairs().
509: The 2-norm of the eigenvector is one unless the problem is generalized
510: Hermitian. In this case the eigenvector is normalized with respect to the
511: norm defined by the B matrix.
513: Level: beginner
515: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
516: @*/
517: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
518: {
519: PetscInt nconv;
521: PetscFunctionBegin;
526: EPSCheckSolved(eps,1);
527: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
528: PetscCall(EPS_GetActualConverged(eps,&nconv));
529: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
530: PetscCall(EPSComputeVectors(eps));
531: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*@
536: EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
538: Collective
540: Input Parameters:
541: + eps - eigensolver context
542: - i - index of the solution
544: Output Parameters:
545: + Wr - real part of left eigenvector
546: - Wi - imaginary part of left eigenvector
548: Notes:
549: The caller must provide valid Vec objects, i.e., they must be created
550: by the calling program with e.g. MatCreateVecs().
552: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
553: configured with complex scalars the eigenvector is stored directly in Wr
554: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
555: one of them is not required.
557: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
558: Eigensolutions are indexed according to the ordering criterion established
559: with EPSSetWhichEigenpairs().
561: Left eigenvectors are available only if the twosided flag was set, see
562: EPSSetTwoSided().
564: Level: intermediate
566: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
567: @*/
568: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
569: {
570: PetscInt nconv;
571: PetscBool trivial;
572: Mat H;
573: IS is[2];
574: Vec v;
576: PetscFunctionBegin;
581: EPSCheckSolved(eps,1);
582: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
583: PetscCall(EPS_GetActualConverged(eps,&nconv));
584: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
586: trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE)? PETSC_TRUE: PETSC_FALSE;
587: if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
589: PetscCall(EPSComputeVectors(eps));
590: if (trivial) {
591: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
592: if (eps->problem_type==EPS_BSE) { /* change sign of bottom part of the vector */
593: PetscCall(STGetMatrix(eps->st,0,&H));
594: PetscCall(MatNestGetISs(H,is,NULL));
595: if (Wr) {
596: PetscCall(VecGetSubVector(Wr,is[1],&v));
597: PetscCall(VecScale(v,-1.0));
598: PetscCall(VecRestoreSubVector(Wr,is[1],&v));
599: }
600: #if !defined(PETSC_USE_COMPLEX)
601: if (Wi) {
602: PetscCall(VecGetSubVector(Wi,is[1],&v));
603: PetscCall(VecScale(v,-1.0));
604: PetscCall(VecRestoreSubVector(Wi,is[1],&v));
605: }
606: #endif
607: }
608: } else {
609: PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
610: }
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /*@
615: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
616: computed eigenpair.
618: Not Collective
620: Input Parameters:
621: + eps - eigensolver context
622: - i - index of eigenpair
624: Output Parameter:
625: . errest - the error estimate
627: Notes:
628: This is the error estimate used internally by the eigensolver. The actual
629: error bound can be computed with EPSComputeError(). See also the users
630: manual for details.
632: Level: advanced
634: .seealso: EPSComputeError()
635: @*/
636: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
637: {
638: PetscInt nconv;
640: PetscFunctionBegin;
642: PetscAssertPointer(errest,3);
643: EPSCheckSolved(eps,1);
644: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
645: PetscCall(EPS_GetActualConverged(eps,&nconv));
646: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
647: if (nconv==eps->nconv) {
648: *errest = eps->errest[eps->perm[i]];
649: } else {
650: PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
651: /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
652: *errest = eps->errest[eps->perm[i/2]];
653: }
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: /*
658: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
659: associated with an eigenpair.
661: Input Parameters:
662: trans - whether A' must be used instead of A
663: kr,ki - eigenvalue
664: xr,xi - eigenvector
665: z - three work vectors (the second one not referenced in complex scalars)
666: */
667: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
668: {
669: PetscInt nmat;
670: Mat A,B;
671: Vec u,w;
672: PetscScalar alpha;
673: #if !defined(PETSC_USE_COMPLEX)
674: Vec v;
675: PetscReal ni,nr;
676: #endif
677: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
679: PetscFunctionBegin;
680: u = z[0]; w = z[2];
681: PetscCall(STGetNumMatrices(eps->st,&nmat));
682: PetscCall(STGetMatrix(eps->st,0,&A));
683: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
685: #if !defined(PETSC_USE_COMPLEX)
686: v = z[1];
687: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
688: #endif
689: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
690: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
691: if (nmat>1) PetscCall((*matmult)(B,xr,w));
692: else PetscCall(VecCopy(xr,w)); /* w=B*x */
693: alpha = trans? -PetscConj(kr): -kr;
694: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
695: }
696: PetscCall(VecNorm(u,NORM_2,norm));
697: #if !defined(PETSC_USE_COMPLEX)
698: } else {
699: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
700: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
701: if (nmat>1) PetscCall((*matmult)(B,xr,v));
702: else PetscCall(VecCopy(xr,v)); /* v=B*xr */
703: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
704: if (nmat>1) PetscCall((*matmult)(B,xi,w));
705: else PetscCall(VecCopy(xi,w)); /* w=B*xi */
706: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
707: }
708: PetscCall(VecNorm(u,NORM_2,&nr));
709: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
710: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
711: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
712: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
713: }
714: PetscCall(VecNorm(u,NORM_2,&ni));
715: *norm = SlepcAbsEigenvalue(nr,ni);
716: }
717: #endif
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: /*@
722: EPSComputeError - Computes the error (based on the residual norm) associated
723: with the i-th computed eigenpair.
725: Collective
727: Input Parameters:
728: + eps - the eigensolver context
729: . i - the solution index
730: - type - the type of error to compute
732: Output Parameter:
733: . error - the error
735: Notes:
736: The error can be computed in various ways, all of them based on the residual
737: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
739: Level: beginner
741: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
742: @*/
743: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
744: {
745: Mat A,B;
746: Vec xr,xi,w[3];
747: PetscReal t,vecnorm=1.0,errorl;
748: PetscScalar kr,ki;
749: PetscBool flg;
751: PetscFunctionBegin;
755: PetscAssertPointer(error,4);
756: EPSCheckSolved(eps,1);
758: /* allocate work vectors */
759: #if defined(PETSC_USE_COMPLEX)
760: PetscCall(EPSSetWorkVecs(eps,3));
761: xi = NULL;
762: w[1] = NULL;
763: #else
764: PetscCall(EPSSetWorkVecs(eps,5));
765: xi = eps->work[3];
766: w[1] = eps->work[4];
767: #endif
768: xr = eps->work[0];
769: w[0] = eps->work[1];
770: w[2] = eps->work[2];
772: /* compute residual norm */
773: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
774: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
776: /* compute 2-norm of eigenvector */
777: if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
779: /* if two-sided, compute left residual norm and take the maximum */
780: if (eps->twosided) {
781: PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
782: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
783: *error = PetscMax(*error,errorl);
784: }
786: /* compute error */
787: switch (type) {
788: case EPS_ERROR_ABSOLUTE:
789: break;
790: case EPS_ERROR_RELATIVE:
791: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
792: break;
793: case EPS_ERROR_BACKWARD:
794: /* initialization of matrix norms */
795: if (!eps->nrma) {
796: PetscCall(STGetMatrix(eps->st,0,&A));
797: PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
798: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
799: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
800: }
801: if (eps->isgeneralized) {
802: if (!eps->nrmb) {
803: PetscCall(STGetMatrix(eps->st,1,&B));
804: PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
805: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
806: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
807: }
808: } else eps->nrmb = 1.0;
809: t = SlepcAbsEigenvalue(kr,ki);
810: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
811: break;
812: default:
813: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
814: }
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: /*
819: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
820: for the recurrence that builds the right subspace.
822: Collective
824: Input Parameters:
825: + eps - the eigensolver context
826: - i - iteration number
828: Output Parameters:
829: . breakdown - flag indicating that a breakdown has occurred
831: Notes:
832: The start vector is computed from another vector: for the first step (i=0),
833: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
834: vector is created. Then this vector is forced to be in the range of OP (only
835: for generalized definite problems) and orthonormalized with respect to all
836: V-vectors up to i-1. The resulting vector is placed in V[i].
838: The flag breakdown is set to true if either i=0 and the vector belongs to the
839: deflation space, or i>0 and the vector is linearly dependent with respect
840: to the V-vectors.
841: */
842: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
843: {
844: PetscReal norm;
845: PetscBool lindep;
846: Vec w,z;
848: PetscFunctionBegin;
852: /* For the first step, use the first initial vector, otherwise a random one */
853: if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
855: /* Force the vector to be in the range of OP for generalized problems with B-inner product */
856: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
857: PetscCall(BVCreateVec(eps->V,&w));
858: PetscCall(BVCopyVec(eps->V,i,w));
859: PetscCall(BVGetColumn(eps->V,i,&z));
860: PetscCall(STApply(eps->st,w,z));
861: PetscCall(BVRestoreColumn(eps->V,i,&z));
862: PetscCall(VecDestroy(&w));
863: }
865: /* Orthonormalize the vector with respect to previous vectors */
866: PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
867: if (breakdown) *breakdown = lindep;
868: else if (lindep || norm == 0.0) {
869: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
870: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
871: }
872: PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*
877: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
878: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
879: */
880: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
881: {
882: PetscReal norm;
883: PetscBool lindep;
884: Vec w,z;
886: PetscFunctionBegin;
890: /* For the first step, use the first initial vector, otherwise a random one */
891: if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
893: /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
894: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
895: PetscCall(BVCreateVec(eps->W,&w));
896: PetscCall(BVCopyVec(eps->W,i,w));
897: PetscCall(BVGetColumn(eps->W,i,&z));
898: PetscCall(STApplyHermitianTranspose(eps->st,w,z));
899: PetscCall(BVRestoreColumn(eps->W,i,&z));
900: PetscCall(VecDestroy(&w));
901: }
903: /* Orthonormalize the vector with respect to previous vectors */
904: PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
905: if (breakdown) *breakdown = lindep;
906: else if (lindep || norm == 0.0) {
907: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
908: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
909: }
910: PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }