Actual source code: epssolve.c
slepc-3.21.0 2024-03-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: #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
29: static PetscErrorCode EPSComputeValues(EPS eps)
30: {
31: PetscBool injective,iscomp,isfilter;
32: PetscInt i,n,aux,nconv0;
33: Mat A,B=NULL,G,Z;
35: PetscFunctionBegin;
36: switch (eps->categ) {
37: case EPS_CATEGORY_KRYLOV:
38: case EPS_CATEGORY_OTHER:
39: PetscCall(STIsInjective(eps->st,&injective));
40: if (injective) {
41: /* one-to-one mapping: backtransform eigenvalues */
42: PetscUseTypeMethod(eps,backtransform);
43: } else {
44: /* compute eigenvalues from Rayleigh quotient */
45: PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
46: if (!n) break;
47: PetscCall(EPSGetOperators(eps,&A,&B));
48: PetscCall(BVSetActiveColumns(eps->V,0,n));
49: PetscCall(DSGetCompact(eps->ds,&iscomp));
50: PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
51: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
52: PetscCall(BVMatProject(eps->V,A,eps->V,G));
53: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
54: if (B) {
55: PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
56: PetscCall(BVMatProject(eps->V,B,eps->V,G));
57: PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
58: }
59: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
60: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
61: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
62: PetscCall(DSSetCompact(eps->ds,iscomp));
63: if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
64: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
65: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
66: PetscCall(BVMultInPlace(eps->V,Z,0,n));
67: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
68: }
69: /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
70: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
71: if (isfilter) {
72: nconv0 = eps->nconv;
73: for (i=0;i<eps->nconv;i++) {
74: if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
75: eps->nconv--;
76: if (i<eps->nconv) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
77: }
78: }
79: if (nconv0>eps->nconv) PetscCall(PetscInfo(eps,"Discarded %" PetscInt_FMT " computed eigenvalues lying outside the interval\n",nconv0-eps->nconv));
80: }
81: }
82: break;
83: case EPS_CATEGORY_PRECOND:
84: case EPS_CATEGORY_CONTOUR:
85: /* eigenvalues already available as an output of the solver */
86: break;
87: }
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: /*@
92: EPSSolve - Solves the eigensystem.
94: Collective
96: Input Parameter:
97: . eps - eigensolver context obtained from EPSCreate()
99: Options Database Keys:
100: + -eps_view - print information about the solver used
101: . -eps_view_mat0 - view the first matrix (A)
102: . -eps_view_mat1 - view the second matrix (B)
103: . -eps_view_vectors - view the computed eigenvectors
104: . -eps_view_values - view the computed eigenvalues
105: . -eps_converged_reason - print reason for convergence, and number of iterations
106: . -eps_error_absolute - print absolute errors of each eigenpair
107: . -eps_error_relative - print relative errors of each eigenpair
108: - -eps_error_backward - print backward errors of each eigenpair
110: Notes:
111: All the command-line options listed above admit an optional argument specifying
112: the viewer type and options. For instance, use '-eps_view_mat0 binary:amatrix.bin'
113: to save the A matrix to a binary file, '-eps_view_values draw' to draw the computed
114: eigenvalues graphically, or '-eps_error_relative :myerr.m:ascii_matlab' to save
115: the errors in a file that can be executed in Matlab.
117: Level: beginner
119: .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
120: @*/
121: PetscErrorCode EPSSolve(EPS eps)
122: {
123: PetscInt i;
124: PetscBool hasname;
125: STMatMode matmode;
126: Mat A,B;
128: PetscFunctionBegin;
130: if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
131: PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
133: /* Call setup */
134: PetscCall(EPSSetUp(eps));
135: eps->nconv = 0;
136: eps->its = 0;
137: for (i=0;i<eps->ncv;i++) {
138: eps->eigr[i] = 0.0;
139: eps->eigi[i] = 0.0;
140: eps->errest[i] = 0.0;
141: eps->perm[i] = i;
142: }
143: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
144: PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
146: /* Call solver */
147: PetscUseTypeMethod(eps,solve);
148: PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
149: eps->state = EPS_STATE_SOLVED;
151: /* Only the first nconv columns contain useful information (except in CISS) */
152: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
153: if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
155: /* If inplace, purify eigenvectors before reverting operator */
156: PetscCall(STGetMatMode(eps->st,&matmode));
157: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
158: PetscCall(STPostSolve(eps->st));
160: /* Map eigenvalues back to the original problem if appropriate */
161: PetscCall(EPSComputeValues(eps));
163: #if !defined(PETSC_USE_COMPLEX)
164: /* Reorder conjugate eigenvalues (positive imaginary first) */
165: for (i=0;i<eps->nconv-1;i++) {
166: if (eps->eigi[i] != 0) {
167: if (eps->eigi[i] < 0) {
168: eps->eigi[i] = -eps->eigi[i];
169: eps->eigi[i+1] = -eps->eigi[i+1];
170: /* the next correction only works with eigenvectors */
171: PetscCall(EPSComputeVectors(eps));
172: PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
173: }
174: i++;
175: }
176: }
177: #endif
179: /* Sort eigenvalues according to eps->which parameter */
180: PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
181: PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
183: /* Various viewers */
184: PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
185: PetscCall(EPSConvergedReasonViewFromOptions(eps));
186: PetscCall(EPSErrorViewFromOptions(eps));
187: PetscCall(EPSValuesViewFromOptions(eps));
188: PetscCall(EPSVectorsViewFromOptions(eps));
190: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
191: if (hasname) {
192: PetscCall(EPSGetOperators(eps,&A,NULL));
193: PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
194: }
195: if (eps->isgeneralized) {
196: PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
197: if (hasname) {
198: PetscCall(EPSGetOperators(eps,NULL,&B));
199: PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
200: }
201: }
203: /* Remove deflation and initial subspaces */
204: if (eps->nds) {
205: PetscCall(BVSetNumConstraints(eps->V,0));
206: eps->nds = 0;
207: }
208: eps->nini = 0;
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /*@
213: EPSGetIterationNumber - Gets the current iteration number. If the
214: call to EPSSolve() is complete, then it returns the number of iterations
215: carried out by the solution method.
217: Not Collective
219: Input Parameter:
220: . eps - the eigensolver context
222: Output Parameter:
223: . its - number of iterations
225: Note:
226: During the i-th iteration this call returns i-1. If EPSSolve() is
227: complete, then parameter "its" contains either the iteration number at
228: which convergence was successfully reached, or failure was detected.
229: Call EPSGetConvergedReason() to determine if the solver converged or
230: failed and why.
232: Level: intermediate
234: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
235: @*/
236: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
237: {
238: PetscFunctionBegin;
240: PetscAssertPointer(its,2);
241: *its = eps->its;
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*@
246: EPSGetConverged - Gets the number of converged eigenpairs.
248: Not Collective
250: Input Parameter:
251: . eps - the eigensolver context
253: Output Parameter:
254: . nconv - number of converged eigenpairs
256: Note:
257: This function should be called after EPSSolve() has finished.
259: Level: beginner
261: .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
262: @*/
263: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
264: {
265: PetscFunctionBegin;
267: PetscAssertPointer(nconv,2);
268: EPSCheckSolved(eps,1);
269: *nconv = eps->nconv;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
275: stopped.
277: Not Collective
279: Input Parameter:
280: . eps - the eigensolver context
282: Output Parameter:
283: . reason - negative value indicates diverged, positive value converged
285: Options Database Key:
286: . -eps_converged_reason - print the reason to a viewer
288: Notes:
289: Possible values for reason are
290: + EPS_CONVERGED_TOL - converged up to tolerance
291: . EPS_CONVERGED_USER - converged due to a user-defined condition
292: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
293: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
294: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
296: Can only be called after the call to EPSSolve() is complete.
298: Level: intermediate
300: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
301: @*/
302: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
303: {
304: PetscFunctionBegin;
306: PetscAssertPointer(reason,2);
307: EPSCheckSolved(eps,1);
308: *reason = eps->reason;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@
313: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
314: subspace.
316: Collective
318: Input Parameter:
319: . eps - the eigensolver context
321: Output Parameter:
322: . v - an array of vectors
324: Notes:
325: This function should be called after EPSSolve() has finished.
327: The user should provide in v an array of nconv vectors, where nconv is
328: the value returned by EPSGetConverged().
330: The first k vectors returned in v span an invariant subspace associated
331: with the first k computed eigenvalues (note that this is not true if the
332: k-th eigenvalue is complex and matrix A is real; in this case the first
333: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
334: in X for all x in X (a similar definition applies for generalized
335: eigenproblems).
337: Level: intermediate
339: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
340: @*/
341: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
342: {
343: PetscInt i;
344: BV V=eps->V;
345: Vec w;
347: PetscFunctionBegin;
349: PetscAssertPointer(v,2);
351: EPSCheckSolved(eps,1);
352: PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
353: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
354: PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
355: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
356: PetscCall(BVCopy(eps->V,V));
357: for (i=0;i<eps->nconv;i++) {
358: PetscCall(BVGetColumn(V,i,&w));
359: PetscCall(VecPointwiseDivide(w,w,eps->D));
360: PetscCall(BVRestoreColumn(V,i,&w));
361: }
362: PetscCall(BVOrthogonalize(V,NULL));
363: }
364: for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
365: if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /*@C
370: EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
371: EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
373: Collective
375: Input Parameters:
376: + eps - eigensolver context
377: - i - index of the solution
379: Output Parameters:
380: + eigr - real part of eigenvalue
381: . eigi - imaginary part of eigenvalue
382: . Vr - real part of eigenvector
383: - Vi - imaginary part of eigenvector
385: Notes:
386: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
387: required. Otherwise, the caller must provide valid Vec objects, i.e.,
388: they must be created by the calling program with e.g. MatCreateVecs().
390: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
391: configured with complex scalars the eigenvalue is stored
392: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
393: set to zero). In both cases, the user can pass NULL in eigi and Vi.
395: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
396: Eigenpairs are indexed according to the ordering criterion established
397: with EPSSetWhichEigenpairs().
399: The 2-norm of the eigenvector is one unless the problem is generalized
400: Hermitian. In this case the eigenvector is normalized with respect to the
401: norm defined by the B matrix.
403: Level: beginner
405: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
406: EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
407: @*/
408: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
409: {
410: PetscFunctionBegin;
413: EPSCheckSolved(eps,1);
414: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
415: PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
416: PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
417: if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /*@C
422: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
424: Not Collective
426: Input Parameters:
427: + eps - eigensolver context
428: - i - index of the solution
430: Output Parameters:
431: + eigr - real part of eigenvalue
432: - eigi - imaginary part of eigenvalue
434: Notes:
435: If the eigenvalue is real, then eigi is set to zero. If PETSc is
436: configured with complex scalars the eigenvalue is stored
437: directly in eigr (eigi is set to zero).
439: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
440: Eigenpairs are indexed according to the ordering criterion established
441: with EPSSetWhichEigenpairs().
443: Level: beginner
445: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
446: @*/
447: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
448: {
449: PetscInt k;
451: PetscFunctionBegin;
453: EPSCheckSolved(eps,1);
454: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
455: PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
456: k = eps->perm[i];
457: #if defined(PETSC_USE_COMPLEX)
458: if (eigr) *eigr = eps->eigr[k];
459: if (eigi) *eigi = 0;
460: #else
461: if (eigr) *eigr = eps->eigr[k];
462: if (eigi) *eigi = eps->eigi[k];
463: #endif
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: /*@
468: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
470: Collective
472: Input Parameters:
473: + eps - eigensolver context
474: - i - index of the solution
476: Output Parameters:
477: + Vr - real part of eigenvector
478: - Vi - imaginary part of eigenvector
480: Notes:
481: The caller must provide valid Vec objects, i.e., they must be created
482: by the calling program with e.g. MatCreateVecs().
484: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
485: configured with complex scalars the eigenvector is stored
486: directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
487: or Vi if one of them is not required.
489: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
490: Eigenpairs are indexed according to the ordering criterion established
491: with EPSSetWhichEigenpairs().
493: The 2-norm of the eigenvector is one unless the problem is generalized
494: Hermitian. In this case the eigenvector is normalized with respect to the
495: norm defined by the B matrix.
497: Level: beginner
499: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
500: @*/
501: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
502: {
503: PetscInt k;
505: PetscFunctionBegin;
510: EPSCheckSolved(eps,1);
511: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
512: PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
513: PetscCall(EPSComputeVectors(eps));
514: k = eps->perm[i];
515: PetscCall(BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi));
516: PetscFunctionReturn(PETSC_SUCCESS);
517: }
519: /*@
520: EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
522: Collective
524: Input Parameters:
525: + eps - eigensolver context
526: - i - index of the solution
528: Output Parameters:
529: + Wr - real part of left eigenvector
530: - Wi - imaginary part of left eigenvector
532: Notes:
533: The caller must provide valid Vec objects, i.e., they must be created
534: by the calling program with e.g. MatCreateVecs().
536: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
537: configured with complex scalars the eigenvector is stored directly in Wr
538: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
539: one of them is not required.
541: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
542: Eigensolutions are indexed according to the ordering criterion established
543: with EPSSetWhichEigenpairs().
545: Left eigenvectors are available only if the twosided flag was set, see
546: EPSSetTwoSided().
548: Level: intermediate
550: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
551: @*/
552: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
553: {
554: PetscInt k;
556: PetscFunctionBegin;
561: EPSCheckSolved(eps,1);
562: PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
563: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
564: PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
565: PetscCall(EPSComputeVectors(eps));
566: k = eps->perm[i];
567: PetscCall(BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: /*@
572: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
573: computed eigenpair.
575: Not Collective
577: Input Parameters:
578: + eps - eigensolver context
579: - i - index of eigenpair
581: Output Parameter:
582: . errest - the error estimate
584: Notes:
585: This is the error estimate used internally by the eigensolver. The actual
586: error bound can be computed with EPSComputeError(). See also the users
587: manual for details.
589: Level: advanced
591: .seealso: EPSComputeError()
592: @*/
593: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
594: {
595: PetscFunctionBegin;
597: PetscAssertPointer(errest,3);
598: EPSCheckSolved(eps,1);
599: PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
600: PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
601: *errest = eps->errest[eps->perm[i]];
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*
606: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
607: associated with an eigenpair.
609: Input Parameters:
610: trans - whether A' must be used instead of A
611: kr,ki - eigenvalue
612: xr,xi - eigenvector
613: z - three work vectors (the second one not referenced in complex scalars)
614: */
615: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
616: {
617: PetscInt nmat;
618: Mat A,B;
619: Vec u,w;
620: PetscScalar alpha;
621: #if !defined(PETSC_USE_COMPLEX)
622: Vec v;
623: PetscReal ni,nr;
624: #endif
625: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
627: PetscFunctionBegin;
628: u = z[0]; w = z[2];
629: PetscCall(STGetNumMatrices(eps->st,&nmat));
630: PetscCall(STGetMatrix(eps->st,0,&A));
631: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
633: #if !defined(PETSC_USE_COMPLEX)
634: v = z[1];
635: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
636: #endif
637: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
638: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
639: if (nmat>1) PetscCall((*matmult)(B,xr,w));
640: else PetscCall(VecCopy(xr,w)); /* w=B*x */
641: alpha = trans? -PetscConj(kr): -kr;
642: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
643: }
644: PetscCall(VecNorm(u,NORM_2,norm));
645: #if !defined(PETSC_USE_COMPLEX)
646: } else {
647: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
648: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
649: if (nmat>1) PetscCall((*matmult)(B,xr,v));
650: else PetscCall(VecCopy(xr,v)); /* v=B*xr */
651: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
652: if (nmat>1) PetscCall((*matmult)(B,xi,w));
653: else PetscCall(VecCopy(xi,w)); /* w=B*xi */
654: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
655: }
656: PetscCall(VecNorm(u,NORM_2,&nr));
657: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
658: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
659: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
660: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
661: }
662: PetscCall(VecNorm(u,NORM_2,&ni));
663: *norm = SlepcAbsEigenvalue(nr,ni);
664: }
665: #endif
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: EPSComputeError - Computes the error (based on the residual norm) associated
671: with the i-th computed eigenpair.
673: Collective
675: Input Parameters:
676: + eps - the eigensolver context
677: . i - the solution index
678: - type - the type of error to compute
680: Output Parameter:
681: . error - the error
683: Notes:
684: The error can be computed in various ways, all of them based on the residual
685: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
687: Level: beginner
689: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
690: @*/
691: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
692: {
693: Mat A,B;
694: Vec xr,xi,w[3];
695: PetscReal t,vecnorm=1.0,errorl;
696: PetscScalar kr,ki;
697: PetscBool flg;
699: PetscFunctionBegin;
703: PetscAssertPointer(error,4);
704: EPSCheckSolved(eps,1);
706: /* allocate work vectors */
707: #if defined(PETSC_USE_COMPLEX)
708: PetscCall(EPSSetWorkVecs(eps,3));
709: xi = NULL;
710: w[1] = NULL;
711: #else
712: PetscCall(EPSSetWorkVecs(eps,5));
713: xi = eps->work[3];
714: w[1] = eps->work[4];
715: #endif
716: xr = eps->work[0];
717: w[0] = eps->work[1];
718: w[2] = eps->work[2];
720: /* compute residual norm */
721: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
722: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
724: /* compute 2-norm of eigenvector */
725: if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
727: /* if two-sided, compute left residual norm and take the maximum */
728: if (eps->twosided) {
729: PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
730: PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
731: *error = PetscMax(*error,errorl);
732: }
734: /* compute error */
735: switch (type) {
736: case EPS_ERROR_ABSOLUTE:
737: break;
738: case EPS_ERROR_RELATIVE:
739: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
740: break;
741: case EPS_ERROR_BACKWARD:
742: /* initialization of matrix norms */
743: if (!eps->nrma) {
744: PetscCall(STGetMatrix(eps->st,0,&A));
745: PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
746: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
747: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
748: }
749: if (eps->isgeneralized) {
750: if (!eps->nrmb) {
751: PetscCall(STGetMatrix(eps->st,1,&B));
752: PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
753: PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
754: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
755: }
756: } else eps->nrmb = 1.0;
757: t = SlepcAbsEigenvalue(kr,ki);
758: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
759: break;
760: default:
761: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
762: }
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: /*
767: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
768: for the recurrence that builds the right subspace.
770: Collective
772: Input Parameters:
773: + eps - the eigensolver context
774: - i - iteration number
776: Output Parameters:
777: . breakdown - flag indicating that a breakdown has occurred
779: Notes:
780: The start vector is computed from another vector: for the first step (i=0),
781: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
782: vector is created. Then this vector is forced to be in the range of OP (only
783: for generalized definite problems) and orthonormalized with respect to all
784: V-vectors up to i-1. The resulting vector is placed in V[i].
786: The flag breakdown is set to true if either i=0 and the vector belongs to the
787: deflation space, or i>0 and the vector is linearly dependent with respect
788: to the V-vectors.
789: */
790: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
791: {
792: PetscReal norm;
793: PetscBool lindep;
794: Vec w,z;
796: PetscFunctionBegin;
800: /* For the first step, use the first initial vector, otherwise a random one */
801: if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
803: /* Force the vector to be in the range of OP for definite generalized problems */
804: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
805: PetscCall(BVCreateVec(eps->V,&w));
806: PetscCall(BVCopyVec(eps->V,i,w));
807: PetscCall(BVGetColumn(eps->V,i,&z));
808: PetscCall(STApply(eps->st,w,z));
809: PetscCall(BVRestoreColumn(eps->V,i,&z));
810: PetscCall(VecDestroy(&w));
811: }
813: /* Orthonormalize the vector with respect to previous vectors */
814: PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
815: if (breakdown) *breakdown = lindep;
816: else if (lindep || norm == 0.0) {
817: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
818: PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
819: }
820: PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
821: PetscFunctionReturn(PETSC_SUCCESS);
822: }
824: /*
825: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
826: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
827: */
828: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
829: {
830: PetscReal norm;
831: PetscBool lindep;
833: PetscFunctionBegin;
837: /* For the first step, use the first initial vector, otherwise a random one */
838: if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
840: /* Orthonormalize the vector with respect to previous vectors */
841: PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
842: if (breakdown) *breakdown = lindep;
843: else if (lindep || norm == 0.0) {
844: PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
845: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
846: }
847: PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }