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 - 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));
134: /* Safeguard for matrices of size 0 */
135: if (eps->n == 0) {
136: eps->nconv = 0;
137: eps->reason = EPS_CONVERGED_TOL;
138: eps->state = EPS_STATE_SOLVED;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: eps->nconv = 0;
143: eps->its = 0;
144: for (i=0;i<eps->ncv;i++) {
145: eps->eigr[i] = 0.0;
146: eps->eigi[i] = 0.0;
147: eps->errest[i] = 0.0;
148: eps->perm[i] = i;
149: }
150: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
151: PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
153: /* Call solver */
154: PetscUseTypeMethod(eps,solve);
155: PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
156: eps->state = EPS_STATE_SOLVED;
158: /* Only the first nconv columns contain useful information (except in CISS) */
159: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
160: if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
162: /* If inplace, purify eigenvectors before reverting operator */
163: PetscCall(STGetMatMode(eps->st,&matmode));
164: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
165: PetscCall(STPostSolve(eps->st));
167: /* Map eigenvalues back to the original problem if appropriate */
168: PetscCall(EPSComputeValues(eps));
170: #if !defined(PETSC_USE_COMPLEX)
171: /* Reorder conjugate eigenvalues (positive imaginary first) */
172: for (i=0;i<eps->nconv-1;i++) {
173: if (eps->eigi[i] != 0 && (eps->problem_type!=EPS_HAMILT || eps->eigr[i]!=0)) {
174: /* conjugate eigenvalues */
175: if (eps->eigi[i] < 0) {
176: eps->eigi[i] = -eps->eigi[i];
177: eps->eigi[i+1] = -eps->eigi[i+1];
178: /* the next correction only works with eigenvectors */
179: PetscCall(EPSComputeVectors(eps));
180: PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
181: if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
182: }
183: i++;
184: }
185: }
186: #endif
188: /* Sort eigenvalues according to eps->which parameter */
189: if (eps->problem_type==EPS_HAMILT) PetscCall(SlepcSortEigenvaluesSpecial(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
190: else PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
191: PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
193: /* Various viewers */
194: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
195: PetscCall(EPSConvergedReasonViewFromOptions(eps));
196: PetscCall(EPSErrorViewFromOptions(eps));
197: PetscCall(EPSValuesViewFromOptions(eps));
198: PetscCall(EPSVectorsViewFromOptions(eps));
200: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
201: if (hasname) {
202: PetscCall(EPSGetOperators(eps,&A,NULL));
203: PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
204: }
205: if (eps->isgeneralized) {
206: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
207: if (hasname) {
208: PetscCall(EPSGetOperators(eps,NULL,&B));
209: PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
210: }
211: }
213: /* Remove deflation and initial subspaces */
214: if (eps->nds) {
215: PetscCall(BVSetNumConstraints(eps->V,0));
216: eps->nds = 0;
217: }
218: eps->nini = 0;
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@
223: EPSGetIterationNumber - Gets the current iteration number. If the
224: call to EPSSolve() is complete, then it returns the number of iterations
225: carried out by the solution method.
227: Not Collective
229: Input Parameter:
230: . eps - the eigensolver context
232: Output Parameter:
233: . its - number of iterations
235: Note:
236: During the i-th iteration this call returns i-1. If EPSSolve() is
237: complete, then parameter "its" contains either the iteration number at
238: which convergence was successfully reached, or failure was detected.
239: Call EPSGetConvergedReason() to determine if the solver converged or
240: failed and why.
242: Level: intermediate
244: .seealso: `EPSGetConvergedReason()`, `EPSSetTolerances()`
245: @*/
246: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
247: {
248: PetscFunctionBegin;
250: PetscAssertPointer(its,2);
251: *its = eps->its;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@
256: EPSGetConverged - Gets the number of converged eigenpairs.
258: Not Collective
260: Input Parameter:
261: . eps - the eigensolver context
263: Output Parameter:
264: . nconv - number of converged eigenpairs
266: Note:
267: This function should be called after EPSSolve() has finished.
269: Level: beginner
271: .seealso: `EPSSetDimensions()`, `EPSSolve()`, `EPSGetEigenpair()`
272: @*/
273: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
274: {
275: PetscFunctionBegin;
277: PetscAssertPointer(nconv,2);
278: EPSCheckSolved(eps,1);
279: PetscCall(EPS_GetActualConverged(eps,nconv));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*@
284: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
285: stopped.
287: Not Collective
289: Input Parameter:
290: . eps - the eigensolver context
292: Output Parameter:
293: . reason - negative value indicates diverged, positive value converged
295: Options Database Key:
296: . -eps_converged_reason - print the reason to a viewer
298: Notes:
299: Possible values for reason are
300: + EPS_CONVERGED_TOL - converged up to tolerance
301: . EPS_CONVERGED_USER - converged due to a user-defined condition
302: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
303: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
304: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
306: Can only be called after the call to EPSSolve() is complete.
308: Level: intermediate
310: .seealso: `EPSSetTolerances()`, `EPSSolve()`, `EPSConvergedReason`
311: @*/
312: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
313: {
314: PetscFunctionBegin;
316: PetscAssertPointer(reason,2);
317: EPSCheckSolved(eps,1);
318: *reason = eps->reason;
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@
323: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
324: subspace.
326: Collective
328: Input Parameter:
329: . eps - the eigensolver context
331: Output Parameter:
332: . v - an array of vectors
334: Notes:
335: This function should be called after EPSSolve() has finished.
337: The user should provide in v an array of nconv vectors, where nconv is
338: the value returned by EPSGetConverged().
340: The first k vectors returned in v span an invariant subspace associated
341: with the first k computed eigenvalues (note that this is not true if the
342: k-th eigenvalue is complex and matrix A is real; in this case the first
343: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
344: in X for all x in X (a similar definition applies for generalized
345: eigenproblems).
347: Level: intermediate
349: .seealso: `EPSGetEigenpair()`, `EPSGetConverged()`, `EPSSolve()`
350: @*/
351: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
352: {
353: PetscInt i;
354: BV V=eps->V;
355: Vec w;
357: PetscFunctionBegin;
359: PetscAssertPointer(v,2);
361: EPSCheckSolved(eps,1);
362: PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
363: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
364: PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
365: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
366: PetscCall(BVCopy(eps->V,V));
367: for (i=0;i<eps->nconv;i++) {
368: PetscCall(BVGetColumn(V,i,&w));
369: PetscCall(VecPointwiseDivide(w,w,eps->D));
370: PetscCall(BVRestoreColumn(V,i,&w));
371: }
372: PetscCall(BVOrthogonalize(V,NULL));
373: }
374: for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
375: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /*@
380: EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
381: EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
383: Collective
385: Input Parameters:
386: + eps - eigensolver context
387: - i - index of the solution
389: Output Parameters:
390: + eigr - real part of eigenvalue
391: . eigi - imaginary part of eigenvalue
392: . Vr - real part of eigenvector
393: - Vi - imaginary part of eigenvector
395: Notes:
396: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
397: required. Otherwise, the caller must provide valid Vec objects, i.e.,
398: they must be created by the calling program with e.g. MatCreateVecs().
400: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
401: configured with complex scalars the eigenvalue is stored
402: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
403: set to zero). In both cases, the user can pass NULL in eigi and Vi.
405: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
406: Eigenpairs are indexed according to the ordering criterion established
407: with EPSSetWhichEigenpairs().
409: The 2-norm of the eigenvector is one unless the problem is generalized
410: Hermitian. In this case the eigenvector is normalized with respect to the
411: norm defined by the B matrix.
413: Level: beginner
415: .seealso: `EPSGetEigenvalue()`, `EPSGetEigenvector()`, `EPSGetLeftEigenvector()`, `EPSSolve()`,
416: `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetInvariantSubspace()`
417: @*/
418: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
419: {
420: PetscInt nconv;
422: PetscFunctionBegin;
425: EPSCheckSolved(eps,1);
426: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
427: PetscCall(EPS_GetActualConverged(eps,&nconv));
428: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
429: PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
430: if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*@
435: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
437: Not Collective
439: Input Parameters:
440: + eps - eigensolver context
441: - i - index of the solution
443: Output Parameters:
444: + eigr - real part of eigenvalue
445: - eigi - imaginary part of eigenvalue
447: Notes:
448: If the eigenvalue is real, then eigi is set to zero. If PETSc is
449: configured with complex scalars the eigenvalue is stored
450: directly in eigr (eigi is set to zero).
452: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
453: Eigenpairs are indexed according to the ordering criterion established
454: with EPSSetWhichEigenpairs().
456: Level: beginner
458: .seealso: `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`
459: @*/
460: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
461: {
462: PetscInt k,nconv;
463: #if !defined(PETSC_USE_COMPLEX)
464: PetscInt k2, iquad;
465: #endif
467: PetscFunctionBegin;
469: EPSCheckSolved(eps,1);
470: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
471: PetscCall(EPS_GetActualConverged(eps,&nconv));
472: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
473: if (nconv==eps->nconv) {
474: k = eps->perm[i];
475: #if defined(PETSC_USE_COMPLEX)
476: if (eigr) *eigr = eps->eigr[k];
477: if (eigi) *eigi = 0;
478: #else
479: if (eigr) *eigr = eps->eigr[k];
480: if (eigi) *eigi = eps->eigi[k];
481: #endif
482: } else {
483: PetscCheck(eps->problem_type==EPS_BSE || eps->problem_type==EPS_HAMILT,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE or Hamiltonian");
484: if (eps->problem_type==EPS_BSE) {
485: /* BSE problem, even index is +lambda, odd index is -lambda */
486: k = eps->perm[i/2];
487: #if defined(PETSC_USE_COMPLEX)
488: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
489: if (eigi) *eigi = 0;
490: #else
491: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
492: if (eigi) *eigi = eps->eigi[k];
493: #endif
494: } else if (eps->problem_type==EPS_HAMILT) {
495: /* Hamiltonian eigenproblem */
496: k = eps->perm[i/2];
497: #if defined(PETSC_USE_COMPLEX)
498: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
499: if (eigi) *eigi = 0;
500: #else
501: if (eps->eigi[k]==0.0) { /* real eigenvalue */
502: if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
503: if (eigi) *eigi = 0.0;
504: } else if (eps->eigr[k]==0.0) { /* purely imaginary eigenvalue */
505: if (eigr) *eigr = 0.0;
506: if (eigi) *eigi = (i%2)? -eps->eigi[k]: eps->eigi[k];
507: } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
508: iquad = i%2; /* index within the 4 values */
509: if (i>1) {
510: k2 = eps->perm[(i-2)/2];
511: if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
512: }
513: if (eigr) *eigr = (iquad<2)? -eps->eigr[k]: eps->eigr[k];
514: if (eigi) *eigi = (iquad%3)? -eps->eigi[k]: eps->eigi[k];
515: }
516: #endif
517: }
518: }
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
525: Collective
527: Input Parameters:
528: + eps - eigensolver context
529: - i - index of the solution
531: Output Parameters:
532: + Vr - real part of eigenvector
533: - Vi - imaginary part of eigenvector
535: Notes:
536: The caller must provide valid Vec objects, i.e., they must be created
537: by the calling program with e.g. MatCreateVecs().
539: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
540: configured with complex scalars the eigenvector is stored
541: directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
542: or Vi if one of them is not required.
544: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
545: Eigenpairs are indexed according to the ordering criterion established
546: with EPSSetWhichEigenpairs().
548: The 2-norm of the eigenvector is one unless the problem is generalized
549: Hermitian. In this case the eigenvector is normalized with respect to the
550: norm defined by the B matrix.
552: Level: beginner
554: .seealso: `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`, `EPSGetLeftEigenvector()`
555: @*/
556: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
557: {
558: PetscInt nconv;
560: PetscFunctionBegin;
565: EPSCheckSolved(eps,1);
566: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
567: PetscCall(EPS_GetActualConverged(eps,&nconv));
568: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
569: PetscCall(EPSComputeVectors(eps));
570: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: /*@
575: EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
577: Collective
579: Input Parameters:
580: + eps - eigensolver context
581: - i - index of the solution
583: Output Parameters:
584: + Wr - real part of left eigenvector
585: - Wi - imaginary part of left eigenvector
587: Notes:
588: The caller must provide valid Vec objects, i.e., they must be created
589: by the calling program with e.g. MatCreateVecs().
591: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
592: configured with complex scalars the eigenvector is stored directly in Wr
593: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
594: one of them is not required.
596: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
597: Eigensolutions are indexed according to the ordering criterion established
598: with EPSSetWhichEigenpairs().
600: Left eigenvectors are available only if the twosided flag was set, see
601: EPSSetTwoSided().
603: Level: intermediate
605: .seealso: `EPSGetEigenvector()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSSetTwoSided()`
606: @*/
607: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
608: {
609: PetscInt nconv;
610: PetscBool trivial;
611: Mat H;
612: IS is[2];
613: Vec v;
615: PetscFunctionBegin;
620: EPSCheckSolved(eps,1);
621: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
622: PetscCall(EPS_GetActualConverged(eps,&nconv));
623: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
625: trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE)? PETSC_TRUE: PETSC_FALSE;
626: if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
628: PetscCall(EPSComputeVectors(eps));
629: if (trivial) {
630: PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
631: if (eps->problem_type==EPS_BSE) { /* change sign of bottom part of the vector */
632: PetscCall(STGetMatrix(eps->st,0,&H));
633: PetscCall(MatNestGetISs(H,is,NULL));
634: if (Wr) {
635: PetscCall(VecGetSubVector(Wr,is[1],&v));
636: PetscCall(VecScale(v,-1.0));
637: PetscCall(VecRestoreSubVector(Wr,is[1],&v));
638: }
639: #if !defined(PETSC_USE_COMPLEX)
640: if (Wi) {
641: PetscCall(VecGetSubVector(Wi,is[1],&v));
642: PetscCall(VecScale(v,-1.0));
643: PetscCall(VecRestoreSubVector(Wi,is[1],&v));
644: }
645: #endif
646: }
647: } else {
648: PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
649: }
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: /*@
654: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
655: computed eigenpair.
657: Not Collective
659: Input Parameters:
660: + eps - eigensolver context
661: - i - index of eigenpair
663: Output Parameter:
664: . errest - the error estimate
666: Notes:
667: This is the error estimate used internally by the eigensolver. The actual
668: error bound can be computed with EPSComputeError(). See also the users
669: manual for details.
671: Level: advanced
673: .seealso: `EPSComputeError()`
674: @*/
675: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
676: {
677: PetscInt nconv;
679: PetscFunctionBegin;
681: PetscAssertPointer(errest,3);
682: EPSCheckSolved(eps,1);
683: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
684: PetscCall(EPS_GetActualConverged(eps,&nconv));
685: PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
686: if (nconv==eps->nconv) {
687: *errest = eps->errest[eps->perm[i]];
688: } else {
689: PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
690: /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
691: *errest = eps->errest[eps->perm[i/2]];
692: }
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }
696: /*
697: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
698: associated with an eigenpair.
700: Input Parameters:
701: trans - whether A' must be used instead of A
702: kr,ki - eigenvalue
703: xr,xi - eigenvector
704: z - three work vectors (the second one not referenced in complex scalars)
705: */
706: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
707: {
708: PetscInt nmat;
709: Mat A,B;
710: Vec u,w;
711: PetscScalar alpha;
712: #if !defined(PETSC_USE_COMPLEX)
713: Vec v;
714: PetscReal ni,nr;
715: #endif
716: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
718: PetscFunctionBegin;
719: u = z[0]; w = z[2];
720: PetscCall(STGetNumMatrices(eps->st,&nmat));
721: PetscCall(STGetMatrix(eps->st,0,&A));
722: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
724: #if !defined(PETSC_USE_COMPLEX)
725: v = z[1];
726: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
727: #endif
728: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
729: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
730: if (nmat>1) PetscCall((*matmult)(B,xr,w));
731: else PetscCall(VecCopy(xr,w)); /* w=B*x */
732: alpha = trans? -PetscConj(kr): -kr;
733: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
734: }
735: PetscCall(VecNorm(u,NORM_2,norm));
736: #if !defined(PETSC_USE_COMPLEX)
737: } else {
738: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
739: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
740: if (nmat>1) PetscCall((*matmult)(B,xr,v));
741: else PetscCall(VecCopy(xr,v)); /* v=B*xr */
742: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
743: if (nmat>1) PetscCall((*matmult)(B,xi,w));
744: else PetscCall(VecCopy(xi,w)); /* w=B*xi */
745: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
746: }
747: PetscCall(VecNorm(u,NORM_2,&nr));
748: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
749: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
750: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
751: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
752: }
753: PetscCall(VecNorm(u,NORM_2,&ni));
754: *norm = SlepcAbsEigenvalue(nr,ni);
755: }
756: #endif
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*@
761: EPSComputeError - Computes the error (based on the residual norm) associated
762: with the i-th computed eigenpair.
764: Collective
766: Input Parameters:
767: + eps - the eigensolver context
768: . i - the solution index
769: - type - the type of error to compute
771: Output Parameter:
772: . error - the error
774: Notes:
775: The error can be computed in various ways, all of them based on the residual
776: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
778: Level: beginner
780: .seealso: `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`
781: @*/
782: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
783: {
784: Mat A,B;
785: Vec xr,xi,w[3];
786: PetscReal t,vecnorm=1.0,errorl;
787: PetscScalar kr,ki;
788: PetscBool flg;
790: PetscFunctionBegin;
794: PetscAssertPointer(error,4);
795: EPSCheckSolved(eps,1);
797: /* allocate work vectors */
798: #if defined(PETSC_USE_COMPLEX)
799: PetscCall(EPSSetWorkVecs(eps,3));
800: xi = NULL;
801: w[1] = NULL;
802: #else
803: PetscCall(EPSSetWorkVecs(eps,5));
804: xi = eps->work[3];
805: w[1] = eps->work[4];
806: #endif
807: xr = eps->work[0];
808: w[0] = eps->work[1];
809: w[2] = eps->work[2];
811: /* compute residual norm */
812: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
813: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
815: /* compute 2-norm of eigenvector */
816: if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
818: /* if two-sided, compute left residual norm and take the maximum */
819: if (eps->twosided) {
820: PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
821: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
822: *error = PetscMax(*error,errorl);
823: }
825: /* compute error */
826: switch (type) {
827: case EPS_ERROR_ABSOLUTE:
828: break;
829: case EPS_ERROR_RELATIVE:
830: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
831: break;
832: case EPS_ERROR_BACKWARD:
833: /* initialization of matrix norms */
834: if (!eps->nrma) {
835: PetscCall(STGetMatrix(eps->st,0,&A));
836: PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
837: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
838: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
839: }
840: if (eps->isgeneralized) {
841: if (!eps->nrmb) {
842: PetscCall(STGetMatrix(eps->st,1,&B));
843: PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
844: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
845: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
846: }
847: } else eps->nrmb = 1.0;
848: t = SlepcAbsEigenvalue(kr,ki);
849: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
850: break;
851: default:
852: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
853: }
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: /*
858: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
859: for the recurrence that builds the right subspace.
861: Collective
863: Input Parameters:
864: + eps - the eigensolver context
865: - i - iteration number
867: Output Parameters:
868: . breakdown - flag indicating that a breakdown has occurred
870: Notes:
871: The start vector is computed from another vector: for the first step (i=0),
872: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
873: vector is created. Then this vector is forced to be in the range of OP (only
874: for generalized definite problems) and orthonormalized with respect to all
875: V-vectors up to i-1. The resulting vector is placed in V[i].
877: The flag breakdown is set to true if either i=0 and the vector belongs to the
878: deflation space, or i>0 and the vector is linearly dependent with respect
879: to the V-vectors.
880: */
881: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
882: {
883: PetscReal norm;
884: PetscBool lindep;
885: Vec w,z;
887: PetscFunctionBegin;
891: /* For the first step, use the first initial vector, otherwise a random one */
892: if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
894: /* Force the vector to be in the range of OP for generalized problems with B-inner product */
895: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
896: PetscCall(BVCreateVec(eps->V,&w));
897: PetscCall(BVCopyVec(eps->V,i,w));
898: PetscCall(BVGetColumn(eps->V,i,&z));
899: PetscCall(STApply(eps->st,w,z));
900: PetscCall(BVRestoreColumn(eps->V,i,&z));
901: PetscCall(VecDestroy(&w));
902: }
904: /* Orthonormalize the vector with respect to previous vectors */
905: PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
906: if (breakdown) *breakdown = lindep;
907: else if (lindep || norm == 0.0) {
908: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
909: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
910: }
911: PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: /*
916: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
917: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
918: */
919: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
920: {
921: PetscReal norm;
922: PetscBool lindep;
923: Vec w,z;
925: PetscFunctionBegin;
929: /* For the first step, use the first initial vector, otherwise a random one */
930: if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
932: /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
933: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
934: PetscCall(BVCreateVec(eps->W,&w));
935: PetscCall(BVCopyVec(eps->W,i,w));
936: PetscCall(BVGetColumn(eps->W,i,&z));
937: PetscCall(STApplyHermitianTranspose(eps->st,w,z));
938: PetscCall(BVRestoreColumn(eps->W,i,&z));
939: PetscCall(VecDestroy(&w));
940: }
942: /* Orthonormalize the vector with respect to previous vectors */
943: PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
944: if (breakdown) *breakdown = lindep;
945: else if (lindep || norm == 0.0) {
946: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
947: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
948: }
949: PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }