1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepceps.h" I*/
15: #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
16: #include <petscdraw.h>
18: PetscErrorCode EPSComputeVectors(EPS eps) 19: {
23: EPSCheckSolved(eps,1);
24: if (eps->state==EPS_STATE_SOLVED && eps->ops->computevectors) {
25: (*eps->ops->computevectors)(eps);
26: }
27: eps->state = EPS_STATE_EIGENVECTORS;
28: return(0);
29: }
31: #define SWAP(a,b,t) {t=a;a=b;b=t;} 33: static PetscErrorCode EPSComputeValues(EPS eps) 34: {
36: PetscBool injective,iscomp,isfilter;
37: PetscInt i,n,aux,nconv0;
38: Mat A,B=NULL,G,Z;
41: switch (eps->categ) {
42: case EPS_CATEGORY_KRYLOV:
43: case EPS_CATEGORY_OTHER:
44: STIsInjective(eps->st,&injective);
45: if (injective) {
46: /* one-to-one mapping: backtransform eigenvalues */
47: if (eps->ops->backtransform) {
48: (*eps->ops->backtransform)(eps);
49: } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, spectral transform should have a backtransform operation");
50: } else {
51: /* compute eigenvalues from Rayleigh quotient */
52: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);
53: if (!n) break;
54: EPSGetOperators(eps,&A,&B);
55: BVSetActiveColumns(eps->V,0,n);
56: DSGetCompact(eps->ds,&iscomp);
57: DSSetCompact(eps->ds,PETSC_FALSE);
58: DSGetMat(eps->ds,DS_MAT_A,&G);
59: BVMatProject(eps->V,A,eps->V,G);
60: DSRestoreMat(eps->ds,DS_MAT_A,&G);
61: if (B) {
62: DSGetMat(eps->ds,DS_MAT_B,&G);
63: BVMatProject(eps->V,B,eps->V,G);
64: DSRestoreMat(eps->ds,DS_MAT_A,&G);
65: }
66: DSSolve(eps->ds,eps->eigr,eps->eigi);
67: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
68: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
69: DSSetCompact(eps->ds,iscomp);
70: if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
71: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
72: DSGetMat(eps->ds,DS_MAT_X,&Z);
73: BVMultInPlace(eps->V,Z,0,n);
74: MatDestroy(&Z);
75: }
76: /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
77: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter);
78: if (isfilter) {
79: nconv0 = eps->nconv;
80: for (i=0;i<eps->nconv;i++) {
81: if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
82: eps->nconv--;
83: if (i<eps->nconv) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
84: }
85: }
86: if (nconv0>eps->nconv) {
87: PetscInfo1(eps,"Discarded %D computed eigenvalues lying outside the interval\n",nconv0-eps->nconv);
88: }
89: }
90: }
91: break;
92: case EPS_CATEGORY_PRECOND:
93: case EPS_CATEGORY_CONTOUR:
94: /* eigenvalues already available as an output of the solver */
95: break;
96: }
97: return(0);
98: }
100: /*@
101: EPSSolve - Solves the eigensystem.
103: Collective on EPS105: Input Parameter:
106: . eps - eigensolver context obtained from EPSCreate()
108: Options Database Keys:
109: + -eps_view - print information about the solver used
110: . -eps_view_mat0 binary - save the first matrix (A) to the default binary viewer
111: . -eps_view_mat1 binary - save the second matrix (B) to the default binary viewer
112: . -eps_view_vectors binary - save the computed eigenvectors to the default binary viewer
113: . -eps_view_values - print computed eigenvalues
114: . -eps_converged_reason - print reason for convergence, and number of iterations
115: . -eps_error_absolute - print absolute errors of each eigenpair
116: . -eps_error_relative - print relative errors of each eigenpair
117: - -eps_error_backward - print backward errors of each eigenpair
119: Level: beginner
121: .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
122: @*/
123: PetscErrorCode EPSSolve(EPS eps)124: {
126: PetscInt i;
127: STMatMode matmode;
128: Mat A,B;
132: if (eps->state>=EPS_STATE_SOLVED) return(0);
133: PetscLogEventBegin(EPS_Solve,eps,0,0,0);
135: /* call setup */
136: EPSSetUp(eps);
137: eps->nconv = 0;
138: eps->its = 0;
139: for (i=0;i<eps->ncv;i++) {
140: eps->eigr[i] = 0.0;
141: eps->eigi[i] = 0.0;
142: eps->errest[i] = 0.0;
143: eps->perm[i] = i;
144: }
145: EPSViewFromOptions(eps,NULL,"-eps_view_pre");
146: RGViewFromOptions(eps->rg,NULL,"-rg_view");
148: /* call solver */
149: (*eps->ops->solve)(eps);
150: eps->state = EPS_STATE_SOLVED;
152: STGetMatMode(eps->st,&matmode);
153: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) {
154: /* Purify eigenvectors before reverting operator */
155: EPSComputeVectors(eps);
156: }
157: STPostSolve(eps->st);
159: if (!eps->reason) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
161: /* Map eigenvalues back to the original problem if appropriate */
162: EPSComputeValues(eps);
164: #if !defined(PETSC_USE_COMPLEX)
165: /* reorder conjugate eigenvalues (positive imaginary first) */
166: for (i=0;i<eps->nconv-1;i++) {
167: if (eps->eigi[i] != 0) {
168: if (eps->eigi[i] < 0) {
169: eps->eigi[i] = -eps->eigi[i];
170: eps->eigi[i+1] = -eps->eigi[i+1];
171: /* the next correction only works with eigenvectors */
172: EPSComputeVectors(eps);
173: BVScaleColumn(eps->V,i+1,-1.0);
174: }
175: i++;
176: }
177: }
178: #endif
180: /* sort eigenvalues according to eps->which parameter */
181: SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm);
182: PetscLogEventEnd(EPS_Solve,eps,0,0,0);
184: /* various viewers */
185: EPSViewFromOptions(eps,NULL,"-eps_view");
186: EPSReasonViewFromOptions(eps);
187: EPSErrorViewFromOptions(eps);
188: EPSValuesViewFromOptions(eps);
189: EPSVectorsViewFromOptions(eps);
190: EPSGetOperators(eps,&A,&B);
191: MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0");
192: if (eps->isgeneralized) {
193: MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1");
194: }
196: /* Remove deflation and initial subspaces */
197: if (eps->nds) {
198: BVSetNumConstraints(eps->V,0);
199: eps->nds = 0;
200: }
201: eps->nini = 0;
202: return(0);
203: }
205: /*@
206: EPSGetIterationNumber - Gets the current iteration number. If the
207: call to EPSSolve() is complete, then it returns the number of iterations
208: carried out by the solution method.
210: Not Collective
212: Input Parameter:
213: . eps - the eigensolver context
215: Output Parameter:
216: . its - number of iterations
218: Note:
219: During the i-th iteration this call returns i-1. If EPSSolve() is
220: complete, then parameter "its" contains either the iteration number at
221: which convergence was successfully reached, or failure was detected.
222: Call EPSGetConvergedReason() to determine if the solver converged or
223: failed and why.
225: Level: intermediate
227: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
228: @*/
229: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)230: {
234: *its = eps->its;
235: return(0);
236: }
238: /*@
239: EPSGetConverged - Gets the number of converged eigenpairs.
241: Not Collective
243: Input Parameter:
244: . eps - the eigensolver context
246: Output Parameter:
247: . nconv - number of converged eigenpairs
249: Note:
250: This function should be called after EPSSolve() has finished.
252: Level: beginner
254: .seealso: EPSSetDimensions(), EPSSolve()
255: @*/
256: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)257: {
261: EPSCheckSolved(eps,1);
262: *nconv = eps->nconv;
263: return(0);
264: }
266: /*@
267: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
268: stopped.
270: Not Collective
272: Input Parameter:
273: . eps - the eigensolver context
275: Output Parameter:
276: . reason - negative value indicates diverged, positive value converged
278: Notes:
279: Possible values for reason are
280: + EPS_CONVERGED_TOL - converged up to tolerance
281: . EPS_CONVERGED_USER - converged due to a user-defined condition
282: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
283: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
284: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
286: Can only be called after the call to EPSSolve() is complete.
288: Level: intermediate
290: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason291: @*/
292: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)293: {
297: EPSCheckSolved(eps,1);
298: *reason = eps->reason;
299: return(0);
300: }
302: /*@
303: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
304: subspace.
306: Not Collective, but vectors are shared by all processors that share the EPS308: Input Parameter:
309: . eps - the eigensolver context
311: Output Parameter:
312: . v - an array of vectors
314: Notes:
315: This function should be called after EPSSolve() has finished.
317: The user should provide in v an array of nconv vectors, where nconv is
318: the value returned by EPSGetConverged().
320: The first k vectors returned in v span an invariant subspace associated
321: with the first k computed eigenvalues (note that this is not true if the
322: k-th eigenvalue is complex and matrix A is real; in this case the first
323: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
324: in X for all x in X (a similar definition applies for generalized
325: eigenproblems).
327: Level: intermediate
329: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
330: @*/
331: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec *v)332: {
334: PetscInt i;
340: EPSCheckSolved(eps,1);
341: if (!eps->ishermitian && eps->state==EPS_STATE_EIGENVECTORS) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
342: for (i=0;i<eps->nconv;i++) {
343: BVCopyVec(eps->V,i,v[i]);
344: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
345: VecPointwiseDivide(v[i],v[i],eps->D);
346: VecNormalize(v[i],NULL);
347: }
348: }
349: return(0);
350: }
352: /*@C
353: EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
354: EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
356: Logically Collective on EPS358: Input Parameters:
359: + eps - eigensolver context
360: - i - index of the solution
362: Output Parameters:
363: + eigr - real part of eigenvalue
364: . eigi - imaginary part of eigenvalue
365: . Vr - real part of eigenvector
366: - Vi - imaginary part of eigenvector
368: Notes:
369: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
370: required. Otherwise, the caller must provide valid Vec objects, i.e.,
371: they must be created by the calling program with e.g. MatCreateVecs().
373: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
374: configured with complex scalars the eigenvalue is stored
375: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
376: set to zero). In both cases, the user can pass NULL in eigi and Vi.
378: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
379: Eigenpairs are indexed according to the ordering criterion established
380: with EPSSetWhichEigenpairs().
382: The 2-norm of the eigenvector is one unless the problem is generalized
383: Hermitian. In this case the eigenvector is normalized with respect to the
384: norm defined by the B matrix.
386: Level: beginner
388: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
389: EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
390: @*/
391: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)392: {
398: EPSCheckSolved(eps,1);
399: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
400: EPSGetEigenvalue(eps,i,eigr,eigi);
401: if (Vr || Vi) { EPSGetEigenvector(eps,i,Vr,Vi); }
402: return(0);
403: }
405: /*@C
406: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
408: Not Collective
410: Input Parameters:
411: + eps - eigensolver context
412: - i - index of the solution
414: Output Parameters:
415: + eigr - real part of eigenvalue
416: - eigi - imaginary part of eigenvalue
418: Notes:
419: If the eigenvalue is real, then eigi is set to zero. If PETSc is
420: configured with complex scalars the eigenvalue is stored
421: directly in eigr (eigi is set to zero).
423: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
424: Eigenpairs are indexed according to the ordering criterion established
425: with EPSSetWhichEigenpairs().
427: Level: beginner
429: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
430: @*/
431: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)432: {
433: PetscInt k;
437: EPSCheckSolved(eps,1);
438: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
439: k = eps->perm[i];
440: #if defined(PETSC_USE_COMPLEX)
441: if (eigr) *eigr = eps->eigr[k];
442: if (eigi) *eigi = 0;
443: #else
444: if (eigr) *eigr = eps->eigr[k];
445: if (eigi) *eigi = eps->eigi[k];
446: #endif
447: return(0);
448: }
450: /*@
451: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
453: Logically Collective on EPS455: Input Parameters:
456: + eps - eigensolver context
457: - i - index of the solution
459: Output Parameters:
460: + Vr - real part of eigenvector
461: - Vi - imaginary part of eigenvector
463: Notes:
464: The caller must provide valid Vec objects, i.e., they must be created
465: by the calling program with e.g. MatCreateVecs().
467: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
468: configured with complex scalars the eigenvector is stored
469: directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
470: or Vi if one of them is not required.
472: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
473: Eigenpairs are indexed according to the ordering criterion established
474: with EPSSetWhichEigenpairs().
476: The 2-norm of the eigenvector is one unless the problem is generalized
477: Hermitian. In this case the eigenvector is normalized with respect to the
478: norm defined by the B matrix.
480: Level: beginner
482: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
483: @*/
484: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)485: {
487: PetscInt k;
494: EPSCheckSolved(eps,1);
495: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
496: EPSComputeVectors(eps);
497: k = eps->perm[i];
498: BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi);
499: return(0);
500: }
502: /*@
503: EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
505: Logically Collective on EPS507: Input Parameters:
508: + eps - eigensolver context
509: - i - index of the solution
511: Output Parameters:
512: + Wr - real part of left eigenvector
513: - Wi - imaginary part of left eigenvector
515: Notes:
516: The caller must provide valid Vec objects, i.e., they must be created
517: by the calling program with e.g. MatCreateVecs().
519: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
520: configured with complex scalars the eigenvector is stored directly in Wr
521: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
522: one of them is not required.
524: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
525: Eigensolutions are indexed according to the ordering criterion established
526: with EPSSetWhichEigenpairs().
528: Left eigenvectors are available only if the twosided flag was set, see
529: EPSSetTwoSided().
531: Level: intermediate
533: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
534: @*/
535: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)536: {
538: PetscInt k;
545: EPSCheckSolved(eps,1);
546: if (!eps->twosided) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
547: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
548: EPSComputeVectors(eps);
549: k = eps->perm[i];
550: BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi);
551: return(0);
552: }
554: /*@
555: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
556: computed eigenpair.
558: Not Collective
560: Input Parameter:
561: + eps - eigensolver context
562: - i - index of eigenpair
564: Output Parameter:
565: . errest - the error estimate
567: Notes:
568: This is the error estimate used internally by the eigensolver. The actual
569: error bound can be computed with EPSComputeError(). See also the users
570: manual for details.
572: Level: advanced
574: .seealso: EPSComputeError()
575: @*/
576: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)577: {
581: EPSCheckSolved(eps,1);
582: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
583: *errest = eps->errest[eps->perm[i]];
584: return(0);
585: }
587: /*
588: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
589: associated with an eigenpair.
591: Input Parameters:
592: trans - whether A' must be used instead of A
593: kr,ki - eigenvalue
594: xr,xi - eigenvector
595: z - three work vectors (the second one not referenced in complex scalars)
596: */
597: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)598: {
600: PetscInt nmat;
601: Mat A,B;
602: Vec u,w;
603: PetscScalar alpha;
604: #if !defined(PETSC_USE_COMPLEX)
605: Vec v;
606: PetscReal ni,nr;
607: #endif
608: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
611: u = z[0]; w = z[2];
612: STGetNumMatrices(eps->st,&nmat);
613: STGetMatrix(eps->st,0,&A);
614: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
616: #if !defined(PETSC_USE_COMPLEX)
617: v = z[1];
618: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
619: #endif
620: (*matmult)(A,xr,u); /* u=A*x */
621: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
622: if (nmat>1) { (*matmult)(B,xr,w); }
623: else { VecCopy(xr,w); } /* w=B*x */
624: alpha = trans? -PetscConj(kr): -kr;
625: VecAXPY(u,alpha,w); /* u=A*x-k*B*x */
626: }
627: VecNorm(u,NORM_2,norm);
628: #if !defined(PETSC_USE_COMPLEX)
629: } else {
630: (*matmult)(A,xr,u); /* u=A*xr */
631: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
632: if (nmat>1) { (*matmult)(B,xr,v); }
633: else { VecCopy(xr,v); } /* v=B*xr */
634: VecAXPY(u,-kr,v); /* u=A*xr-kr*B*xr */
635: if (nmat>1) { (*matmult)(B,xi,w); }
636: else { VecCopy(xi,w); } /* w=B*xi */
637: VecAXPY(u,trans?-ki:ki,w); /* u=A*xr-kr*B*xr+ki*B*xi */
638: }
639: VecNorm(u,NORM_2,&nr);
640: (*matmult)(A,xi,u); /* u=A*xi */
641: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
642: VecAXPY(u,-kr,w); /* u=A*xi-kr*B*xi */
643: VecAXPY(u,trans?ki:-ki,v); /* u=A*xi-kr*B*xi-ki*B*xr */
644: }
645: VecNorm(u,NORM_2,&ni);
646: *norm = SlepcAbsEigenvalue(nr,ni);
647: }
648: #endif
649: return(0);
650: }
652: /*@
653: EPSComputeError - Computes the error (based on the residual norm) associated
654: with the i-th computed eigenpair.
656: Collective on EPS658: Input Parameter:
659: + eps - the eigensolver context
660: . i - the solution index
661: - type - the type of error to compute
663: Output Parameter:
664: . error - the error
666: Notes:
667: The error can be computed in various ways, all of them based on the residual
668: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
670: Level: beginner
672: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
673: @*/
674: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)675: {
677: Mat A,B;
678: Vec xr,xi,w[3];
679: PetscReal t,vecnorm=1.0,errorl;
680: PetscScalar kr,ki;
681: PetscBool flg;
688: EPSCheckSolved(eps,1);
690: /* allocate work vectors */
691: #if defined(PETSC_USE_COMPLEX)
692: EPSSetWorkVecs(eps,3);
693: xi = NULL;
694: w[1] = NULL;
695: #else
696: EPSSetWorkVecs(eps,5);
697: xi = eps->work[3];
698: w[1] = eps->work[4];
699: #endif
700: xr = eps->work[0];
701: w[0] = eps->work[1];
702: w[2] = eps->work[2];
704: /* compute residual norm */
705: EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
706: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error);
708: /* compute 2-norm of eigenvector */
709: if (eps->problem_type==EPS_GHEP) {
710: VecNorm(xr,NORM_2,&vecnorm);
711: }
713: /* if two-sided, compute left residual norm and take the maximum */
714: if (eps->twosided) {
715: EPSGetLeftEigenvector(eps,i,xr,xi);
716: EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl);
717: *error = PetscMax(*error,errorl);
718: }
720: /* compute error */
721: switch (type) {
722: case EPS_ERROR_ABSOLUTE:
723: break;
724: case EPS_ERROR_RELATIVE:
725: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
726: break;
727: case EPS_ERROR_BACKWARD:
728: /* initialization of matrix norms */
729: if (!eps->nrma) {
730: STGetMatrix(eps->st,0,&A);
731: MatHasOperation(A,MATOP_NORM,&flg);
732: if (!flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
733: MatNorm(A,NORM_INFINITY,&eps->nrma);
734: }
735: if (eps->isgeneralized) {
736: if (!eps->nrmb) {
737: STGetMatrix(eps->st,1,&B);
738: MatHasOperation(B,MATOP_NORM,&flg);
739: if (!flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
740: MatNorm(B,NORM_INFINITY,&eps->nrmb);
741: }
742: } else eps->nrmb = 1.0;
743: t = SlepcAbsEigenvalue(kr,ki);
744: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
745: break;
746: default:747: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
748: }
749: return(0);
750: }
752: /*
753: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
754: for the recurrence that builds the right subspace.
756: Collective on EPS and Vec
758: Input Parameters:
759: + eps - the eigensolver context
760: - i - iteration number
762: Output Parameters:
763: . breakdown - flag indicating that a breakdown has occurred
765: Notes:
766: The start vector is computed from another vector: for the first step (i=0),
767: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
768: vector is created. Then this vector is forced to be in the range of OP (only
769: for generalized definite problems) and orthonormalized with respect to all
770: V-vectors up to i-1. The resulting vector is placed in V[i].
772: The flag breakdown is set to true if either i=0 and the vector belongs to the
773: deflation space, or i>0 and the vector is linearly dependent with respect
774: to the V-vectors.
775: */
776: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)777: {
779: PetscReal norm;
780: PetscBool lindep;
781: Vec w,z;
787: /* For the first step, use the first initial vector, otherwise a random one */
788: if (i>0 || eps->nini==0) {
789: BVSetRandomColumn(eps->V,i);
790: }
792: /* Force the vector to be in the range of OP for definite generalized problems */
793: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
794: BVCreateVec(eps->V,&w);
795: BVCopyVec(eps->V,i,w);
796: BVGetColumn(eps->V,i,&z);
797: STApply(eps->st,w,z);
798: BVRestoreColumn(eps->V,i,&z);
799: VecDestroy(&w);
800: }
802: /* Orthonormalize the vector with respect to previous vectors */
803: BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep);
804: if (breakdown) *breakdown = lindep;
805: else if (lindep || norm == 0.0) {
806: if (i==0) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Initial vector is zero or belongs to the deflation space");
807: else SETERRQ(PetscObjectComm((PetscObject)eps),1,"Unable to generate more start vectors");
808: }
809: BVScaleColumn(eps->V,i,1.0/norm);
810: return(0);
811: }