Actual source code: epssolve.c
slepc-3.22.1 2024-10-28
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: }
172: i++;
173: }
174: }
175: #endif
177: /* Sort eigenvalues according to eps->which parameter */
178: PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
179: PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
181: /* Various viewers */
182: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
183: PetscCall(EPSConvergedReasonViewFromOptions(eps));
184: PetscCall(EPSErrorViewFromOptions(eps));
185: PetscCall(EPSValuesViewFromOptions(eps));
186: PetscCall(EPSVectorsViewFromOptions(eps));
188: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
189: if (hasname) {
190: PetscCall(EPSGetOperators(eps,&A,NULL));
191: PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
192: }
193: if (eps->isgeneralized) {
194: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
195: if (hasname) {
196: PetscCall(EPSGetOperators(eps,NULL,&B));
197: PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
198: }
199: }
201: /* Remove deflation and initial subspaces */
202: if (eps->nds) {
203: PetscCall(BVSetNumConstraints(eps->V,0));
204: eps->nds = 0;
205: }
206: eps->nini = 0;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: /*@
211: EPSGetIterationNumber - Gets the current iteration number. If the
212: call to EPSSolve() is complete, then it returns the number of iterations
213: carried out by the solution method.
215: Not Collective
217: Input Parameter:
218: . eps - the eigensolver context
220: Output Parameter:
221: . its - number of iterations
223: Note:
224: During the i-th iteration this call returns i-1. If EPSSolve() is
225: complete, then parameter "its" contains either the iteration number at
226: which convergence was successfully reached, or failure was detected.
227: Call EPSGetConvergedReason() to determine if the solver converged or
228: failed and why.
230: Level: intermediate
232: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
233: @*/
234: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
235: {
236: PetscFunctionBegin;
238: PetscAssertPointer(its,2);
239: *its = eps->its;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@
244: EPSGetConverged - Gets the number of converged eigenpairs.
246: Not Collective
248: Input Parameter:
249: . eps - the eigensolver context
251: Output Parameter:
252: . nconv - number of converged eigenpairs
254: Note:
255: This function should be called after EPSSolve() has finished.
257: Level: beginner
259: .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
260: @*/
261: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
262: {
263: PetscFunctionBegin;
265: PetscAssertPointer(nconv,2);
266: EPSCheckSolved(eps,1);
267: PetscCall(EPS_GetActualConverged(eps,nconv));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*@
272: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
273: stopped.
275: Not Collective
277: Input Parameter:
278: . eps - the eigensolver context
280: Output Parameter:
281: . reason - negative value indicates diverged, positive value converged
283: Options Database Key:
284: . -eps_converged_reason - print the reason to a viewer
286: Notes:
287: Possible values for reason are
288: + EPS_CONVERGED_TOL - converged up to tolerance
289: . EPS_CONVERGED_USER - converged due to a user-defined condition
290: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
291: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
292: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
294: Can only be called after the call to EPSSolve() is complete.
296: Level: intermediate
298: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
299: @*/
300: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
301: {
302: PetscFunctionBegin;
304: PetscAssertPointer(reason,2);
305: EPSCheckSolved(eps,1);
306: *reason = eps->reason;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
312: subspace.
314: Collective
316: Input Parameter:
317: . eps - the eigensolver context
319: Output Parameter:
320: . v - an array of vectors
322: Notes:
323: This function should be called after EPSSolve() has finished.
325: The user should provide in v an array of nconv vectors, where nconv is
326: the value returned by EPSGetConverged().
328: The first k vectors returned in v span an invariant subspace associated
329: with the first k computed eigenvalues (note that this is not true if the
330: k-th eigenvalue is complex and matrix A is real; in this case the first
331: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
332: in X for all x in X (a similar definition applies for generalized
333: eigenproblems).
335: Level: intermediate
337: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
338: @*/
339: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
340: {
341: PetscInt i;
342: BV V=eps->V;
343: Vec w;
345: PetscFunctionBegin;
347: PetscAssertPointer(v,2);
349: EPSCheckSolved(eps,1);
350: PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
351: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
352: PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
353: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
354: PetscCall(BVCopy(eps->V,V));
355: for (i=0;i<eps->nconv;i++) {
356: PetscCall(BVGetColumn(V,i,&w));
357: PetscCall(VecPointwiseDivide(w,w,eps->D));
358: PetscCall(BVRestoreColumn(V,i,&w));
359: }
360: PetscCall(BVOrthogonalize(V,NULL));
361: }
362: for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
363: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@
368: EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
369: EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
371: Collective
373: Input Parameters:
374: + eps - eigensolver context
375: - i - index of the solution
377: Output Parameters:
378: + eigr - real part of eigenvalue
379: . eigi - imaginary part of eigenvalue
380: . Vr - real part of eigenvector
381: - Vi - imaginary part of eigenvector
383: Notes:
384: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
385: required. Otherwise, the caller must provide valid Vec objects, i.e.,
386: they must be created by the calling program with e.g. MatCreateVecs().
388: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
389: configured with complex scalars the eigenvalue is stored
390: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
391: set to zero). In both cases, the user can pass NULL in eigi and Vi.
393: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
394: Eigenpairs are indexed according to the ordering criterion established
395: with EPSSetWhichEigenpairs().
397: The 2-norm of the eigenvector is one unless the problem is generalized
398: Hermitian. In this case the eigenvector is normalized with respect to the
399: norm defined by the B matrix.
401: Level: beginner
403: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
404: EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
405: @*/
406: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
407: {
408: PetscInt nconv;
410: PetscFunctionBegin;
413: EPSCheckSolved(eps,1);
414: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
415: PetscCall(EPS_GetActualConverged(eps,&nconv));
416: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
417: PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
418: if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
425: Not Collective
427: Input Parameters:
428: + eps - eigensolver context
429: - i - index of the solution
431: Output Parameters:
432: + eigr - real part of eigenvalue
433: - eigi - imaginary part of eigenvalue
435: Notes:
436: If the eigenvalue is real, then eigi is set to zero. If PETSc is
437: configured with complex scalars the eigenvalue is stored
438: directly in eigr (eigi is set to zero).
440: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
441: Eigenpairs are indexed according to the ordering criterion established
442: with EPSSetWhichEigenpairs().
444: Level: beginner
446: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
447: @*/
448: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
449: {
450: PetscInt k,nconv;
452: PetscFunctionBegin;
454: EPSCheckSolved(eps,1);
455: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
456: PetscCall(EPS_GetActualConverged(eps,&nconv));
457: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
458: if (nconv==eps->nconv) {
459: k = eps->perm[i];
460: #if defined(PETSC_USE_COMPLEX)
461: if (eigr) *eigr = eps->eigr[k];
462: if (eigi) *eigi = 0;
463: #else
464: if (eigr) *eigr = eps->eigr[k];
465: if (eigi) *eigi = eps->eigi[k];
466: #endif
467: } else {
468: PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
469: /* BSE problem, even index is +lambda, odd index is -lambda */
470: k = eps->perm[i/2];
471: #if defined(PETSC_USE_COMPLEX)
472: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
473: if (eigi) *eigi = 0;
474: #else
475: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
476: if (eigi) *eigi = eps->eigi[k];
477: #endif
478: }
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*@
483: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
485: Collective
487: Input Parameters:
488: + eps - eigensolver context
489: - i - index of the solution
491: Output Parameters:
492: + Vr - real part of eigenvector
493: - Vi - imaginary part of eigenvector
495: Notes:
496: The caller must provide valid Vec objects, i.e., they must be created
497: by the calling program with e.g. MatCreateVecs().
499: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
500: configured with complex scalars the eigenvector is stored
501: directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
502: or Vi if one of them is not required.
504: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
505: Eigenpairs are indexed according to the ordering criterion established
506: with EPSSetWhichEigenpairs().
508: The 2-norm of the eigenvector is one unless the problem is generalized
509: Hermitian. In this case the eigenvector is normalized with respect to the
510: norm defined by the B matrix.
512: Level: beginner
514: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
515: @*/
516: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
517: {
518: PetscInt nconv;
520: PetscFunctionBegin;
525: EPSCheckSolved(eps,1);
526: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
527: PetscCall(EPS_GetActualConverged(eps,&nconv));
528: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
529: PetscCall(EPSComputeVectors(eps));
530: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: /*@
535: EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
537: Collective
539: Input Parameters:
540: + eps - eigensolver context
541: - i - index of the solution
543: Output Parameters:
544: + Wr - real part of left eigenvector
545: - Wi - imaginary part of left eigenvector
547: Notes:
548: The caller must provide valid Vec objects, i.e., they must be created
549: by the calling program with e.g. MatCreateVecs().
551: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
552: configured with complex scalars the eigenvector is stored directly in Wr
553: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
554: one of them is not required.
556: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
557: Eigensolutions are indexed according to the ordering criterion established
558: with EPSSetWhichEigenpairs().
560: Left eigenvectors are available only if the twosided flag was set, see
561: EPSSetTwoSided().
563: Level: intermediate
565: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
566: @*/
567: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
568: {
569: PetscInt nconv;
570: PetscBool trivial;
571: Mat H;
572: IS is[2];
573: Vec v;
575: PetscFunctionBegin;
580: EPSCheckSolved(eps,1);
581: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
582: PetscCall(EPS_GetActualConverged(eps,&nconv));
583: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
585: trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE)? PETSC_TRUE: PETSC_FALSE;
586: if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
588: PetscCall(EPSComputeVectors(eps));
589: if (trivial) {
590: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
591: if (eps->problem_type==EPS_BSE) { /* change sign of bottom part of the vector */
592: PetscCall(STGetMatrix(eps->st,0,&H));
593: PetscCall(MatNestGetISs(H,is,NULL));
594: if (Wr) {
595: PetscCall(VecGetSubVector(Wr,is[1],&v));
596: PetscCall(VecScale(v,-1.0));
597: PetscCall(VecRestoreSubVector(Wr,is[1],&v));
598: }
599: #if !defined(PETSC_USE_COMPLEX)
600: if (Wi) {
601: PetscCall(VecGetSubVector(Wi,is[1],&v));
602: PetscCall(VecScale(v,-1.0));
603: PetscCall(VecRestoreSubVector(Wi,is[1],&v));
604: }
605: #endif
606: }
607: } else {
608: PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
609: }
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: /*@
614: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
615: computed eigenpair.
617: Not Collective
619: Input Parameters:
620: + eps - eigensolver context
621: - i - index of eigenpair
623: Output Parameter:
624: . errest - the error estimate
626: Notes:
627: This is the error estimate used internally by the eigensolver. The actual
628: error bound can be computed with EPSComputeError(). See also the users
629: manual for details.
631: Level: advanced
633: .seealso: EPSComputeError()
634: @*/
635: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
636: {
637: PetscInt nconv;
639: PetscFunctionBegin;
641: PetscAssertPointer(errest,3);
642: EPSCheckSolved(eps,1);
643: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
644: PetscCall(EPS_GetActualConverged(eps,&nconv));
645: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
646: if (nconv==eps->nconv) {
647: *errest = eps->errest[eps->perm[i]];
648: } else {
649: PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
650: /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
651: *errest = eps->errest[eps->perm[i/2]];
652: }
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: /*
657: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
658: associated with an eigenpair.
660: Input Parameters:
661: trans - whether A' must be used instead of A
662: kr,ki - eigenvalue
663: xr,xi - eigenvector
664: z - three work vectors (the second one not referenced in complex scalars)
665: */
666: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
667: {
668: PetscInt nmat;
669: Mat A,B;
670: Vec u,w;
671: PetscScalar alpha;
672: #if !defined(PETSC_USE_COMPLEX)
673: Vec v;
674: PetscReal ni,nr;
675: #endif
676: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
678: PetscFunctionBegin;
679: u = z[0]; w = z[2];
680: PetscCall(STGetNumMatrices(eps->st,&nmat));
681: PetscCall(STGetMatrix(eps->st,0,&A));
682: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
684: #if !defined(PETSC_USE_COMPLEX)
685: v = z[1];
686: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
687: #endif
688: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
689: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
690: if (nmat>1) PetscCall((*matmult)(B,xr,w));
691: else PetscCall(VecCopy(xr,w)); /* w=B*x */
692: alpha = trans? -PetscConj(kr): -kr;
693: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
694: }
695: PetscCall(VecNorm(u,NORM_2,norm));
696: #if !defined(PETSC_USE_COMPLEX)
697: } else {
698: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
699: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
700: if (nmat>1) PetscCall((*matmult)(B,xr,v));
701: else PetscCall(VecCopy(xr,v)); /* v=B*xr */
702: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
703: if (nmat>1) PetscCall((*matmult)(B,xi,w));
704: else PetscCall(VecCopy(xi,w)); /* w=B*xi */
705: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
706: }
707: PetscCall(VecNorm(u,NORM_2,&nr));
708: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
709: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
710: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
711: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
712: }
713: PetscCall(VecNorm(u,NORM_2,&ni));
714: *norm = SlepcAbsEigenvalue(nr,ni);
715: }
716: #endif
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: EPSComputeError - Computes the error (based on the residual norm) associated
722: with the i-th computed eigenpair.
724: Collective
726: Input Parameters:
727: + eps - the eigensolver context
728: . i - the solution index
729: - type - the type of error to compute
731: Output Parameter:
732: . error - the error
734: Notes:
735: The error can be computed in various ways, all of them based on the residual
736: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
738: Level: beginner
740: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
741: @*/
742: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
743: {
744: Mat A,B;
745: Vec xr,xi,w[3];
746: PetscReal t,vecnorm=1.0,errorl;
747: PetscScalar kr,ki;
748: PetscBool flg;
750: PetscFunctionBegin;
754: PetscAssertPointer(error,4);
755: EPSCheckSolved(eps,1);
757: /* allocate work vectors */
758: #if defined(PETSC_USE_COMPLEX)
759: PetscCall(EPSSetWorkVecs(eps,3));
760: xi = NULL;
761: w[1] = NULL;
762: #else
763: PetscCall(EPSSetWorkVecs(eps,5));
764: xi = eps->work[3];
765: w[1] = eps->work[4];
766: #endif
767: xr = eps->work[0];
768: w[0] = eps->work[1];
769: w[2] = eps->work[2];
771: /* compute residual norm */
772: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
773: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
775: /* compute 2-norm of eigenvector */
776: if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
778: /* if two-sided, compute left residual norm and take the maximum */
779: if (eps->twosided) {
780: PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
781: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
782: *error = PetscMax(*error,errorl);
783: }
785: /* compute error */
786: switch (type) {
787: case EPS_ERROR_ABSOLUTE:
788: break;
789: case EPS_ERROR_RELATIVE:
790: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
791: break;
792: case EPS_ERROR_BACKWARD:
793: /* initialization of matrix norms */
794: if (!eps->nrma) {
795: PetscCall(STGetMatrix(eps->st,0,&A));
796: PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
797: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
798: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
799: }
800: if (eps->isgeneralized) {
801: if (!eps->nrmb) {
802: PetscCall(STGetMatrix(eps->st,1,&B));
803: PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
804: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
805: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
806: }
807: } else eps->nrmb = 1.0;
808: t = SlepcAbsEigenvalue(kr,ki);
809: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
810: break;
811: default:
812: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
813: }
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: /*
818: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
819: for the recurrence that builds the right subspace.
821: Collective
823: Input Parameters:
824: + eps - the eigensolver context
825: - i - iteration number
827: Output Parameters:
828: . breakdown - flag indicating that a breakdown has occurred
830: Notes:
831: The start vector is computed from another vector: for the first step (i=0),
832: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
833: vector is created. Then this vector is forced to be in the range of OP (only
834: for generalized definite problems) and orthonormalized with respect to all
835: V-vectors up to i-1. The resulting vector is placed in V[i].
837: The flag breakdown is set to true if either i=0 and the vector belongs to the
838: deflation space, or i>0 and the vector is linearly dependent with respect
839: to the V-vectors.
840: */
841: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
842: {
843: PetscReal norm;
844: PetscBool lindep;
845: Vec w,z;
847: PetscFunctionBegin;
851: /* For the first step, use the first initial vector, otherwise a random one */
852: if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
854: /* Force the vector to be in the range of OP for definite generalized problems */
855: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
856: PetscCall(BVCreateVec(eps->V,&w));
857: PetscCall(BVCopyVec(eps->V,i,w));
858: PetscCall(BVGetColumn(eps->V,i,&z));
859: PetscCall(STApply(eps->st,w,z));
860: PetscCall(BVRestoreColumn(eps->V,i,&z));
861: PetscCall(VecDestroy(&w));
862: }
864: /* Orthonormalize the vector with respect to previous vectors */
865: PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
866: if (breakdown) *breakdown = lindep;
867: else if (lindep || norm == 0.0) {
868: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
869: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
870: }
871: PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: /*
876: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
877: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
878: */
879: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
880: {
881: PetscReal norm;
882: PetscBool lindep;
884: PetscFunctionBegin;
888: /* For the first step, use the first initial vector, otherwise a random one */
889: if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
891: /* Orthonormalize the vector with respect to previous vectors */
892: PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
893: if (breakdown) *breakdown = lindep;
894: else if (lindep || norm == 0.0) {
895: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
896: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
897: }
898: PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }