Actual source code: epssolve.c
slepc-3.18.0 2022-10-01
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: EPSCheckSolved(eps,1);
21: if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
22: eps->state = EPS_STATE_EIGENVECTORS;
23: return 0;
24: }
26: #define SWAP(a,b,t) {t=a;a=b;b=t;}
28: static PetscErrorCode EPSComputeValues(EPS eps)
29: {
30: PetscBool injective,iscomp,isfilter;
31: PetscInt i,n,aux,nconv0;
32: Mat A,B=NULL,G,Z;
34: switch (eps->categ) {
35: case EPS_CATEGORY_KRYLOV:
36: case EPS_CATEGORY_OTHER:
37: 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: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL);
44: if (!n) break;
45: EPSGetOperators(eps,&A,&B);
46: BVSetActiveColumns(eps->V,0,n);
47: DSGetCompact(eps->ds,&iscomp);
48: DSSetCompact(eps->ds,PETSC_FALSE);
49: DSGetMat(eps->ds,DS_MAT_A,&G);
50: BVMatProject(eps->V,A,eps->V,G);
51: DSRestoreMat(eps->ds,DS_MAT_A,&G);
52: if (B) {
53: DSGetMat(eps->ds,DS_MAT_B,&G);
54: BVMatProject(eps->V,B,eps->V,G);
55: DSRestoreMat(eps->ds,DS_MAT_A,&G);
56: }
57: DSSolve(eps->ds,eps->eigr,eps->eigi);
58: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
59: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
60: DSSetCompact(eps->ds,iscomp);
61: if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
62: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
63: DSGetMat(eps->ds,DS_MAT_X,&Z);
64: BVMultInPlace(eps->V,Z,0,n);
65: DSRestoreMat(eps->ds,DS_MAT_X,&Z);
66: }
67: /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
68: 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) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
75: }
76: }
77: if (nconv0>eps->nconv) 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: return 0;
87: }
89: /*@
90: EPSSolve - Solves the eigensystem.
92: Collective on eps
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: STMatMode matmode;
123: Mat A,B;
126: if (eps->state>=EPS_STATE_SOLVED) return 0;
127: PetscLogEventBegin(EPS_Solve,eps,0,0,0);
129: /* Call setup */
130: EPSSetUp(eps);
131: eps->nconv = 0;
132: eps->its = 0;
133: for (i=0;i<eps->ncv;i++) {
134: eps->eigr[i] = 0.0;
135: eps->eigi[i] = 0.0;
136: eps->errest[i] = 0.0;
137: eps->perm[i] = i;
138: }
139: EPSViewFromOptions(eps,NULL,"-eps_view_pre");
140: RGViewFromOptions(eps->rg,NULL,"-rg_view");
142: /* Call solver */
143: PetscUseTypeMethod(eps,solve);
145: eps->state = EPS_STATE_SOLVED;
147: /* Only the first nconv columns contain useful information (except in CISS) */
148: BVSetActiveColumns(eps->V,0,eps->nconv);
149: if (eps->twosided) BVSetActiveColumns(eps->W,0,eps->nconv);
151: /* If inplace, purify eigenvectors before reverting operator */
152: STGetMatMode(eps->st,&matmode);
153: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) EPSComputeVectors(eps);
154: STPostSolve(eps->st);
156: /* Map eigenvalues back to the original problem if appropriate */
157: EPSComputeValues(eps);
159: #if !defined(PETSC_USE_COMPLEX)
160: /* Reorder conjugate eigenvalues (positive imaginary first) */
161: for (i=0;i<eps->nconv-1;i++) {
162: if (eps->eigi[i] != 0) {
163: if (eps->eigi[i] < 0) {
164: eps->eigi[i] = -eps->eigi[i];
165: eps->eigi[i+1] = -eps->eigi[i+1];
166: /* the next correction only works with eigenvectors */
167: EPSComputeVectors(eps);
168: BVScaleColumn(eps->V,i+1,-1.0);
169: }
170: i++;
171: }
172: }
173: #endif
175: /* Sort eigenvalues according to eps->which parameter */
176: SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm);
177: PetscLogEventEnd(EPS_Solve,eps,0,0,0);
179: /* Various viewers */
180: EPSViewFromOptions(eps,NULL,"-eps_view");
181: EPSConvergedReasonViewFromOptions(eps);
182: EPSErrorViewFromOptions(eps);
183: EPSValuesViewFromOptions(eps);
184: EPSVectorsViewFromOptions(eps);
185: EPSGetOperators(eps,&A,&B);
186: MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0");
187: if (eps->isgeneralized) MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1");
189: /* Remove deflation and initial subspaces */
190: if (eps->nds) {
191: BVSetNumConstraints(eps->V,0);
192: eps->nds = 0;
193: }
194: eps->nini = 0;
195: return 0;
196: }
198: /*@
199: EPSGetIterationNumber - Gets the current iteration number. If the
200: call to EPSSolve() is complete, then it returns the number of iterations
201: carried out by the solution method.
203: Not Collective
205: Input Parameter:
206: . eps - the eigensolver context
208: Output Parameter:
209: . its - number of iterations
211: Note:
212: During the i-th iteration this call returns i-1. If EPSSolve() is
213: complete, then parameter "its" contains either the iteration number at
214: which convergence was successfully reached, or failure was detected.
215: Call EPSGetConvergedReason() to determine if the solver converged or
216: failed and why.
218: Level: intermediate
220: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
221: @*/
222: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
223: {
226: *its = eps->its;
227: return 0;
228: }
230: /*@
231: EPSGetConverged - Gets the number of converged eigenpairs.
233: Not Collective
235: Input Parameter:
236: . eps - the eigensolver context
238: Output Parameter:
239: . nconv - number of converged eigenpairs
241: Note:
242: This function should be called after EPSSolve() has finished.
244: Level: beginner
246: .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
247: @*/
248: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
249: {
252: EPSCheckSolved(eps,1);
253: *nconv = eps->nconv;
254: return 0;
255: }
257: /*@
258: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
259: stopped.
261: Not Collective
263: Input Parameter:
264: . eps - the eigensolver context
266: Output Parameter:
267: . reason - negative value indicates diverged, positive value converged
269: Options Database Key:
270: . -eps_converged_reason - print the reason to a viewer
272: Notes:
273: Possible values for reason are
274: + EPS_CONVERGED_TOL - converged up to tolerance
275: . EPS_CONVERGED_USER - converged due to a user-defined condition
276: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
277: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
278: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
280: Can only be called after the call to EPSSolve() is complete.
282: Level: intermediate
284: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
285: @*/
286: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
287: {
290: EPSCheckSolved(eps,1);
291: *reason = eps->reason;
292: return 0;
293: }
295: /*@
296: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
297: subspace.
299: Not Collective, but vectors are shared by all processors that share the EPS
301: Input Parameter:
302: . eps - the eigensolver context
304: Output Parameter:
305: . v - an array of vectors
307: Notes:
308: This function should be called after EPSSolve() has finished.
310: The user should provide in v an array of nconv vectors, where nconv is
311: the value returned by EPSGetConverged().
313: The first k vectors returned in v span an invariant subspace associated
314: with the first k computed eigenvalues (note that this is not true if the
315: k-th eigenvalue is complex and matrix A is real; in this case the first
316: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
317: in X for all x in X (a similar definition applies for generalized
318: eigenproblems).
320: Level: intermediate
322: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
323: @*/
324: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
325: {
326: PetscInt i;
327: BV V=eps->V;
328: Vec w;
333: EPSCheckSolved(eps,1);
335: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
336: BVDuplicateResize(eps->V,eps->nconv,&V);
337: BVSetActiveColumns(eps->V,0,eps->nconv);
338: BVCopy(eps->V,V);
339: for (i=0;i<eps->nconv;i++) {
340: BVGetColumn(V,i,&w);
341: VecPointwiseDivide(w,w,eps->D);
342: BVRestoreColumn(V,i,&w);
343: }
344: BVOrthogonalize(V,NULL);
345: }
346: for (i=0;i<eps->nconv;i++) BVCopyVec(V,i,v[i]);
347: if (eps->balance!=EPS_BALANCE_NONE && eps->D) BVDestroy(&V);
348: return 0;
349: }
351: /*@C
352: EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
353: EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
355: Logically Collective on eps
357: Input Parameters:
358: + eps - eigensolver context
359: - i - index of the solution
361: Output Parameters:
362: + eigr - real part of eigenvalue
363: . eigi - imaginary part of eigenvalue
364: . Vr - real part of eigenvector
365: - Vi - imaginary part of eigenvector
367: Notes:
368: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
369: required. Otherwise, the caller must provide valid Vec objects, i.e.,
370: they must be created by the calling program with e.g. MatCreateVecs().
372: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
373: configured with complex scalars the eigenvalue is stored
374: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
375: set to zero). In both cases, the user can pass NULL in eigi and Vi.
377: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
378: Eigenpairs are indexed according to the ordering criterion established
379: with EPSSetWhichEigenpairs().
381: The 2-norm of the eigenvector is one unless the problem is generalized
382: Hermitian. In this case the eigenvector is normalized with respect to the
383: norm defined by the B matrix.
385: Level: beginner
387: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
388: EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
389: @*/
390: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
391: {
394: EPSCheckSolved(eps,1);
397: EPSGetEigenvalue(eps,i,eigr,eigi);
398: if (Vr || Vi) EPSGetEigenvector(eps,i,Vr,Vi);
399: return 0;
400: }
402: /*@C
403: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
405: Not Collective
407: Input Parameters:
408: + eps - eigensolver context
409: - i - index of the solution
411: Output Parameters:
412: + eigr - real part of eigenvalue
413: - eigi - imaginary part of eigenvalue
415: Notes:
416: If the eigenvalue is real, then eigi is set to zero. If PETSc is
417: configured with complex scalars the eigenvalue is stored
418: directly in eigr (eigi is set to zero).
420: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
421: Eigenpairs are indexed according to the ordering criterion established
422: with EPSSetWhichEigenpairs().
424: Level: beginner
426: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
427: @*/
428: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
429: {
430: PetscInt k;
433: EPSCheckSolved(eps,1);
436: k = eps->perm[i];
437: #if defined(PETSC_USE_COMPLEX)
438: if (eigr) *eigr = eps->eigr[k];
439: if (eigi) *eigi = 0;
440: #else
441: if (eigr) *eigr = eps->eigr[k];
442: if (eigi) *eigi = eps->eigi[k];
443: #endif
444: return 0;
445: }
447: /*@
448: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
450: Logically Collective on eps
452: Input Parameters:
453: + eps - eigensolver context
454: - i - index of the solution
456: Output Parameters:
457: + Vr - real part of eigenvector
458: - Vi - imaginary part of eigenvector
460: Notes:
461: The caller must provide valid Vec objects, i.e., they must be created
462: by the calling program with e.g. MatCreateVecs().
464: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
465: configured with complex scalars the eigenvector is stored
466: directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
467: or Vi if one of them is not required.
469: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
470: Eigenpairs are indexed according to the ordering criterion established
471: with EPSSetWhichEigenpairs().
473: The 2-norm of the eigenvector is one unless the problem is generalized
474: Hermitian. In this case the eigenvector is normalized with respect to the
475: norm defined by the B matrix.
477: Level: beginner
479: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
480: @*/
481: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
482: {
483: PetscInt k;
489: EPSCheckSolved(eps,1);
492: EPSComputeVectors(eps);
493: k = eps->perm[i];
494: BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi);
495: return 0;
496: }
498: /*@
499: EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
501: Logically Collective on eps
503: Input Parameters:
504: + eps - eigensolver context
505: - i - index of the solution
507: Output Parameters:
508: + Wr - real part of left eigenvector
509: - Wi - imaginary part of left eigenvector
511: Notes:
512: The caller must provide valid Vec objects, i.e., they must be created
513: by the calling program with e.g. MatCreateVecs().
515: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
516: configured with complex scalars the eigenvector is stored directly in Wr
517: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
518: one of them is not required.
520: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
521: Eigensolutions are indexed according to the ordering criterion established
522: with EPSSetWhichEigenpairs().
524: Left eigenvectors are available only if the twosided flag was set, see
525: EPSSetTwoSided().
527: Level: intermediate
529: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
530: @*/
531: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
532: {
533: PetscInt k;
539: EPSCheckSolved(eps,1);
543: EPSComputeVectors(eps);
544: k = eps->perm[i];
545: BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi);
546: return 0;
547: }
549: /*@
550: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
551: computed eigenpair.
553: Not Collective
555: Input Parameters:
556: + eps - eigensolver context
557: - i - index of eigenpair
559: Output Parameter:
560: . errest - the error estimate
562: Notes:
563: This is the error estimate used internally by the eigensolver. The actual
564: error bound can be computed with EPSComputeError(). See also the users
565: manual for details.
567: Level: advanced
569: .seealso: EPSComputeError()
570: @*/
571: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
572: {
575: EPSCheckSolved(eps,1);
578: *errest = eps->errest[eps->perm[i]];
579: return 0;
580: }
582: /*
583: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
584: associated with an eigenpair.
586: Input Parameters:
587: trans - whether A' must be used instead of A
588: kr,ki - eigenvalue
589: xr,xi - eigenvector
590: z - three work vectors (the second one not referenced in complex scalars)
591: */
592: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
593: {
594: PetscInt nmat;
595: Mat A,B;
596: Vec u,w;
597: PetscScalar alpha;
598: #if !defined(PETSC_USE_COMPLEX)
599: Vec v;
600: PetscReal ni,nr;
601: #endif
602: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
604: u = z[0]; w = z[2];
605: STGetNumMatrices(eps->st,&nmat);
606: STGetMatrix(eps->st,0,&A);
607: if (nmat>1) STGetMatrix(eps->st,1,&B);
609: #if !defined(PETSC_USE_COMPLEX)
610: v = z[1];
611: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
612: #endif
613: (*matmult)(A,xr,u); /* u=A*x */
614: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
615: if (nmat>1) (*matmult)(B,xr,w);
616: else VecCopy(xr,w); /* w=B*x */
617: alpha = trans? -PetscConj(kr): -kr;
618: VecAXPY(u,alpha,w); /* u=A*x-k*B*x */
619: }
620: VecNorm(u,NORM_2,norm);
621: #if !defined(PETSC_USE_COMPLEX)
622: } else {
623: (*matmult)(A,xr,u); /* u=A*xr */
624: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
625: if (nmat>1) (*matmult)(B,xr,v);
626: else VecCopy(xr,v); /* v=B*xr */
627: VecAXPY(u,-kr,v); /* u=A*xr-kr*B*xr */
628: if (nmat>1) (*matmult)(B,xi,w);
629: else VecCopy(xi,w); /* w=B*xi */
630: VecAXPY(u,trans?-ki:ki,w); /* u=A*xr-kr*B*xr+ki*B*xi */
631: }
632: VecNorm(u,NORM_2,&nr);
633: (*matmult)(A,xi,u); /* u=A*xi */
634: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
635: VecAXPY(u,-kr,w); /* u=A*xi-kr*B*xi */
636: VecAXPY(u,trans?ki:-ki,v); /* u=A*xi-kr*B*xi-ki*B*xr */
637: }
638: VecNorm(u,NORM_2,&ni);
639: *norm = SlepcAbsEigenvalue(nr,ni);
640: }
641: #endif
642: return 0;
643: }
645: /*@
646: EPSComputeError - Computes the error (based on the residual norm) associated
647: with the i-th computed eigenpair.
649: Collective on eps
651: Input Parameters:
652: + eps - the eigensolver context
653: . i - the solution index
654: - type - the type of error to compute
656: Output Parameter:
657: . error - the error
659: Notes:
660: The error can be computed in various ways, all of them based on the residual
661: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
663: Level: beginner
665: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
666: @*/
667: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
668: {
669: Mat A,B;
670: Vec xr,xi,w[3];
671: PetscReal t,vecnorm=1.0,errorl;
672: PetscScalar kr,ki;
673: PetscBool flg;
679: EPSCheckSolved(eps,1);
681: /* allocate work vectors */
682: #if defined(PETSC_USE_COMPLEX)
683: EPSSetWorkVecs(eps,3);
684: xi = NULL;
685: w[1] = NULL;
686: #else
687: EPSSetWorkVecs(eps,5);
688: xi = eps->work[3];
689: w[1] = eps->work[4];
690: #endif
691: xr = eps->work[0];
692: w[0] = eps->work[1];
693: w[2] = eps->work[2];
695: /* compute residual norm */
696: EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
697: EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error);
699: /* compute 2-norm of eigenvector */
700: if (eps->problem_type==EPS_GHEP) VecNorm(xr,NORM_2,&vecnorm);
702: /* if two-sided, compute left residual norm and take the maximum */
703: if (eps->twosided) {
704: EPSGetLeftEigenvector(eps,i,xr,xi);
705: EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl);
706: *error = PetscMax(*error,errorl);
707: }
709: /* compute error */
710: switch (type) {
711: case EPS_ERROR_ABSOLUTE:
712: break;
713: case EPS_ERROR_RELATIVE:
714: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
715: break;
716: case EPS_ERROR_BACKWARD:
717: /* initialization of matrix norms */
718: if (!eps->nrma) {
719: STGetMatrix(eps->st,0,&A);
720: MatHasOperation(A,MATOP_NORM,&flg);
722: MatNorm(A,NORM_INFINITY,&eps->nrma);
723: }
724: if (eps->isgeneralized) {
725: if (!eps->nrmb) {
726: STGetMatrix(eps->st,1,&B);
727: MatHasOperation(B,MATOP_NORM,&flg);
729: MatNorm(B,NORM_INFINITY,&eps->nrmb);
730: }
731: } else eps->nrmb = 1.0;
732: t = SlepcAbsEigenvalue(kr,ki);
733: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
734: break;
735: default:
736: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
737: }
738: return 0;
739: }
741: /*
742: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
743: for the recurrence that builds the right subspace.
745: Collective on eps
747: Input Parameters:
748: + eps - the eigensolver context
749: - i - iteration number
751: Output Parameters:
752: . breakdown - flag indicating that a breakdown has occurred
754: Notes:
755: The start vector is computed from another vector: for the first step (i=0),
756: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
757: vector is created. Then this vector is forced to be in the range of OP (only
758: for generalized definite problems) and orthonormalized with respect to all
759: V-vectors up to i-1. The resulting vector is placed in V[i].
761: The flag breakdown is set to true if either i=0 and the vector belongs to the
762: deflation space, or i>0 and the vector is linearly dependent with respect
763: to the V-vectors.
764: */
765: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
766: {
767: PetscReal norm;
768: PetscBool lindep;
769: Vec w,z;
774: /* For the first step, use the first initial vector, otherwise a random one */
775: if (i>0 || eps->nini==0) BVSetRandomColumn(eps->V,i);
777: /* Force the vector to be in the range of OP for definite generalized problems */
778: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
779: BVCreateVec(eps->V,&w);
780: BVCopyVec(eps->V,i,w);
781: BVGetColumn(eps->V,i,&z);
782: STApply(eps->st,w,z);
783: BVRestoreColumn(eps->V,i,&z);
784: VecDestroy(&w);
785: }
787: /* Orthonormalize the vector with respect to previous vectors */
788: BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep);
789: if (breakdown) *breakdown = lindep;
790: else if (lindep || norm == 0.0) {
793: }
794: BVScaleColumn(eps->V,i,1.0/norm);
795: return 0;
796: }
798: /*
799: EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
800: vector for the recurrence that builds the left subspace. See EPSGetStartVector().
801: */
802: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
803: {
804: PetscReal norm;
805: PetscBool lindep;
810: /* For the first step, use the first initial vector, otherwise a random one */
811: if (i>0 || eps->ninil==0) BVSetRandomColumn(eps->W,i);
813: /* Orthonormalize the vector with respect to previous vectors */
814: BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep);
815: if (breakdown) *breakdown = lindep;
816: else if (lindep || norm == 0.0) {
818: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
819: }
820: BVScaleColumn(eps->W,i,1.0/norm);
821: return 0;
822: }