1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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 <petscdraw.h>
17: PetscErrorCode EPSComputeVectors(EPS eps) 18: {
22: EPSCheckSolved(eps,1);
23: if (eps->state==EPS_STATE_SOLVED && eps->ops->computevectors) {
24: (*eps->ops->computevectors)(eps);
25: }
26: eps->state = EPS_STATE_EIGENVECTORS;
27: return(0);
28: }
30: #define SWAP(a,b,t) {t=a;a=b;b=t;} 32: static PetscErrorCode EPSComputeValues(EPS eps) 33: {
35: PetscBool injective,iscomp,isfilter;
36: PetscInt i,n,aux,nconv0;
37: Mat A,B=NULL,G,Z;
40: switch (eps->categ) {
41: case EPS_CATEGORY_KRYLOV:
42: case EPS_CATEGORY_OTHER:
43: STIsInjective(eps->st,&injective);
44: if (injective) {
45: /* one-to-one mapping: backtransform eigenvalues */
46: if (eps->ops->backtransform) {
47: (*eps->ops->backtransform)(eps);
48: } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, spectral transform should have a backtransform operation");
49: } else {
50: /* compute eigenvalues from Rayleigh quotient */
51: DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);
52: if (!n) break;
53: EPSGetOperators(eps,&A,&B);
54: BVSetActiveColumns(eps->V,0,n);
55: DSGetCompact(eps->ds,&iscomp);
56: DSSetCompact(eps->ds,PETSC_FALSE);
57: DSGetMat(eps->ds,DS_MAT_A,&G);
58: BVMatProject(eps->V,A,eps->V,G);
59: DSRestoreMat(eps->ds,DS_MAT_A,&G);
60: if (B) {
61: DSGetMat(eps->ds,DS_MAT_B,&G);
62: BVMatProject(eps->V,B,eps->V,G);
63: DSRestoreMat(eps->ds,DS_MAT_A,&G);
64: }
65: DSSolve(eps->ds,eps->eigr,eps->eigi);
66: DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
67: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
68: DSSetCompact(eps->ds,iscomp);
69: if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
70: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
71: DSGetMat(eps->ds,DS_MAT_X,&Z);
72: BVMultInPlace(eps->V,Z,0,n);
73: MatDestroy(&Z);
74: }
75: /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
76: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter);
77: if (isfilter) {
78: nconv0 = eps->nconv;
79: for (i=0;i<eps->nconv;i++) {
80: if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
81: eps->nconv--;
82: if (i<eps->nconv) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
83: }
84: }
85: if (nconv0>eps->nconv) {
86: PetscInfo1(eps,"Discarded %D computed eigenvalues lying outside the interval\n",nconv0-eps->nconv);
87: }
88: }
89: }
90: break;
91: case EPS_CATEGORY_PRECOND:
92: case EPS_CATEGORY_CONTOUR:
93: /* eigenvalues already available as an output of the solver */
94: break;
95: }
96: return(0);
97: }
99: /*@
100: EPSSolve - Solves the eigensystem.
102: Collective on EPS104: Input Parameter:
105: . eps - eigensolver context obtained from EPSCreate()
107: Options Database Keys:
108: + -eps_view - print information about the solver used
109: . -eps_view_mat0 binary - save the first matrix (A) to the default binary viewer
110: . -eps_view_mat1 binary - save the second matrix (B) to the default binary viewer
111: . -eps_view_vectors binary - save the computed eigenvectors to the default binary viewer
112: . -eps_view_values - print computed eigenvalues
113: . -eps_converged_reason - print reason for convergence, and number of iterations
114: . -eps_error_absolute - print absolute errors of each eigenpair
115: . -eps_error_relative - print relative errors of each eigenpair
116: - -eps_error_backward - print backward errors of each eigenpair
118: Level: beginner
120: .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
121: @*/
122: PetscErrorCode EPSSolve(EPS eps)123: {
125: PetscInt i;
126: STMatMode matmode;
127: Mat A,B;
131: if (eps->state>=EPS_STATE_SOLVED) return(0);
132: PetscLogEventBegin(EPS_Solve,eps,0,0,0);
134: /* call setup */
135: EPSSetUp(eps);
136: eps->nconv = 0;
137: eps->its = 0;
138: for (i=0;i<eps->ncv;i++) {
139: eps->eigr[i] = 0.0;
140: eps->eigi[i] = 0.0;
141: eps->errest[i] = 0.0;
142: eps->perm[i] = i;
143: }
144: EPSViewFromOptions(eps,NULL,"-eps_view_pre");
145: RGViewFromOptions(eps->rg,NULL,"-rg_view");
147: /* call solver */
148: (*eps->ops->solve)(eps);
149: eps->state = EPS_STATE_SOLVED;
151: STGetMatMode(eps->st,&matmode);
152: if (matmode == ST_MATMODE_INPLACE && eps->ispositive) {
153: /* Purify eigenvectors before reverting operator */
154: EPSComputeVectors(eps);
155: }
156: STPostSolve(eps->st);
158: if (!eps->reason) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
160: /* Map eigenvalues back to the original problem if appropriate */
161: 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: EPSComputeVectors(eps);
172: BVScaleColumn(eps->V,i+1,-1.0);
173: }
174: i++;
175: }
176: }
177: #endif
179: /* sort eigenvalues according to eps->which parameter */
180: SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm);
181: PetscLogEventEnd(EPS_Solve,eps,0,0,0);
183: /* various viewers */
184: EPSViewFromOptions(eps,NULL,"-eps_view");
185: EPSReasonViewFromOptions(eps);
186: EPSErrorViewFromOptions(eps);
187: EPSValuesViewFromOptions(eps);
188: EPSVectorsViewFromOptions(eps);
189: EPSGetOperators(eps,&A,&B);
190: MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0");
191: if (eps->isgeneralized) {
192: MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1");
193: }
195: /* Remove deflation and initial subspaces */
196: if (eps->nds) {
197: BVSetNumConstraints(eps->V,0);
198: eps->nds = 0;
199: }
200: eps->nini = 0;
201: return(0);
202: }
204: /*@
205: EPSGetIterationNumber - Gets the current iteration number. If the
206: call to EPSSolve() is complete, then it returns the number of iterations
207: carried out by the solution method.
209: Not Collective
211: Input Parameter:
212: . eps - the eigensolver context
214: Output Parameter:
215: . its - number of iterations
217: Note:
218: During the i-th iteration this call returns i-1. If EPSSolve() is
219: complete, then parameter "its" contains either the iteration number at
220: which convergence was successfully reached, or failure was detected.
221: Call EPSGetConvergedReason() to determine if the solver converged or
222: failed and why.
224: Level: intermediate
226: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
227: @*/
228: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)229: {
233: *its = eps->its;
234: return(0);
235: }
237: /*@
238: EPSGetConverged - Gets the number of converged eigenpairs.
240: Not Collective
242: Input Parameter:
243: . eps - the eigensolver context
245: Output Parameter:
246: . nconv - number of converged eigenpairs
248: Note:
249: This function should be called after EPSSolve() has finished.
251: Level: beginner
253: .seealso: EPSSetDimensions(), EPSSolve()
254: @*/
255: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)256: {
260: EPSCheckSolved(eps,1);
261: *nconv = eps->nconv;
262: return(0);
263: }
265: /*@
266: EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
267: stopped.
269: Not Collective
271: Input Parameter:
272: . eps - the eigensolver context
274: Output Parameter:
275: . reason - negative value indicates diverged, positive value converged
277: Notes:
278: Possible values for reason are
279: + EPS_CONVERGED_TOL - converged up to tolerance
280: . EPS_CONVERGED_USER - converged due to a user-defined condition
281: . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
282: . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
283: - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
285: Can only be called after the call to EPSSolve() is complete.
287: Level: intermediate
289: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason290: @*/
291: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)292: {
296: EPSCheckSolved(eps,1);
297: *reason = eps->reason;
298: return(0);
299: }
301: /*@
302: EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
303: subspace.
305: Not Collective, but vectors are shared by all processors that share the EPS307: Input Parameter:
308: . eps - the eigensolver context
310: Output Parameter:
311: . v - an array of vectors
313: Notes:
314: This function should be called after EPSSolve() has finished.
316: The user should provide in v an array of nconv vectors, where nconv is
317: the value returned by EPSGetConverged().
319: The first k vectors returned in v span an invariant subspace associated
320: with the first k computed eigenvalues (note that this is not true if the
321: k-th eigenvalue is complex and matrix A is real; in this case the first
322: k+1 vectors should be used). An invariant subspace X of A satisfies Ax
323: in X for all x in X (a similar definition applies for generalized
324: eigenproblems).
326: Level: intermediate
328: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
329: @*/
330: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec *v)331: {
333: PetscInt i;
339: EPSCheckSolved(eps,1);
340: 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");
341: for (i=0;i<eps->nconv;i++) {
342: BVCopyVec(eps->V,i,v[i]);
343: if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
344: VecPointwiseDivide(v[i],v[i],eps->D);
345: VecNormalize(v[i],NULL);
346: }
347: }
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 EPS357: 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(), EPSSolve(),
388: EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
389: @*/
390: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)391: {
397: EPSCheckSolved(eps,1);
398: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
399: EPSGetEigenvalue(eps,i,eigr,eigi);
400: if (Vr || Vi) { EPSGetEigenvector(eps,i,Vr,Vi); }
401: return(0);
402: }
404: /*@C
405: EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
407: Not Collective
409: Input Parameters:
410: + eps - eigensolver context
411: - i - index of the solution
413: Output Parameters:
414: + eigr - real part of eigenvalue
415: - eigi - imaginary part of eigenvalue
417: Notes:
418: If the eigenvalue is real, then eigi is set to zero. If PETSc is
419: configured with complex scalars the eigenvalue is stored
420: directly in eigr (eigi is set to zero).
422: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
423: Eigenpairs are indexed according to the ordering criterion established
424: with EPSSetWhichEigenpairs().
426: Level: beginner
428: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
429: @*/
430: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)431: {
432: PetscInt k;
436: EPSCheckSolved(eps,1);
437: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
438: k = eps->perm[i];
439: #if defined(PETSC_USE_COMPLEX)
440: if (eigr) *eigr = eps->eigr[k];
441: if (eigi) *eigi = 0;
442: #else
443: if (eigr) *eigr = eps->eigr[k];
444: if (eigi) *eigi = eps->eigi[k];
445: #endif
446: return(0);
447: }
449: /*@
450: EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
452: Logically Collective on EPS454: Input Parameters:
455: + eps - eigensolver context
456: - i - index of the solution
458: Output Parameters:
459: + Vr - real part of eigenvector
460: - Vi - imaginary part of eigenvector
462: Notes:
463: The caller must provide valid Vec objects, i.e., they must be created
464: by the calling program with e.g. MatCreateVecs().
466: If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
467: configured with complex scalars the eigenvector is stored
468: directly in Vr (Vi is set to zero). In both cases, the user can pass NULL in Vi.
470: The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
471: Eigenpairs are indexed according to the ordering criterion established
472: with EPSSetWhichEigenpairs().
474: The 2-norm of the eigenvector is one unless the problem is generalized
475: Hermitian. In this case the eigenvector is normalized with respect to the
476: norm defined by the B matrix.
478: Level: beginner
480: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
481: @*/
482: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)483: {
485: PetscInt k;
493: EPSCheckSolved(eps,1);
494: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
495: EPSComputeVectors(eps);
496: k = eps->perm[i];
497: #if defined(PETSC_USE_COMPLEX)
498: BVCopyVec(eps->V,k,Vr);
499: if (Vi) { VecSet(Vi,0.0); }
500: #else
501: if (eps->eigi[k] > 0) { /* first value of conjugate pair */
502: BVCopyVec(eps->V,k,Vr);
503: if (Vi) {
504: BVCopyVec(eps->V,k+1,Vi);
505: }
506: } else if (eps->eigi[k] < 0) { /* second value of conjugate pair */
507: BVCopyVec(eps->V,k-1,Vr);
508: if (Vi) {
509: BVCopyVec(eps->V,k,Vi);
510: VecScale(Vi,-1.0);
511: }
512: } else { /* real eigenvalue */
513: BVCopyVec(eps->V,k,Vr);
514: if (Vi) { VecSet(Vi,0.0); }
515: }
516: #endif
517: return(0);
518: }
520: /*@
521: EPSGetErrorEstimate - Returns the error estimate associated to the i-th
522: computed eigenpair.
524: Not Collective
526: Input Parameter:
527: + eps - eigensolver context
528: - i - index of eigenpair
530: Output Parameter:
531: . errest - the error estimate
533: Notes:
534: This is the error estimate used internally by the eigensolver. The actual
535: error bound can be computed with EPSComputeError(). See also the users
536: manual for details.
538: Level: advanced
540: .seealso: EPSComputeError()
541: @*/
542: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)543: {
547: EPSCheckSolved(eps,1);
548: if (i<0 || i>=eps->nconv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
549: *errest = eps->errest[eps->perm[i]];
550: return(0);
551: }
553: /*
554: EPSComputeResidualNorm_Private - Computes the norm of the residual vector
555: associated with an eigenpair.
557: Input Parameters:
558: kr,ki - eigenvalue
559: xr,xi - eigenvector
560: z - three work vectors (the second one not referenced in complex scalars)
561: */
562: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)563: {
565: PetscInt nmat;
566: Mat A,B;
567: Vec u,w;
568: #if !defined(PETSC_USE_COMPLEX)
569: Vec v;
570: PetscReal ni,nr;
571: #endif
574: u = z[0]; w = z[2];
575: STGetNumMatrices(eps->st,&nmat);
576: STGetMatrix(eps->st,0,&A);
577: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
579: #if !defined(PETSC_USE_COMPLEX)
580: v = z[1];
581: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
582: #endif
583: MatMult(A,xr,u); /* u=A*x */
584: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
585: if (nmat>1) { MatMult(B,xr,w); }
586: else { VecCopy(xr,w); } /* w=B*x */
587: VecAXPY(u,-kr,w); /* u=A*x-k*B*x */
588: }
589: VecNorm(u,NORM_2,norm);
590: #if !defined(PETSC_USE_COMPLEX)
591: } else {
592: MatMult(A,xr,u); /* u=A*xr */
593: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
594: if (nmat>1) { MatMult(B,xr,v); }
595: else { VecCopy(xr,v); } /* v=B*xr */
596: VecAXPY(u,-kr,v); /* u=A*xr-kr*B*xr */
597: if (nmat>1) { MatMult(B,xi,w); }
598: else { VecCopy(xi,w); } /* w=B*xi */
599: VecAXPY(u,ki,w); /* u=A*xr-kr*B*xr+ki*B*xi */
600: }
601: VecNorm(u,NORM_2,&nr);
602: MatMult(A,xi,u); /* u=A*xi */
603: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
604: VecAXPY(u,-kr,w); /* u=A*xi-kr*B*xi */
605: VecAXPY(u,-ki,v); /* u=A*xi-kr*B*xi-ki*B*xr */
606: }
607: VecNorm(u,NORM_2,&ni);
608: *norm = SlepcAbsEigenvalue(nr,ni);
609: }
610: #endif
611: return(0);
612: }
614: /*@
615: EPSComputeError - Computes the error (based on the residual norm) associated
616: with the i-th computed eigenpair.
618: Collective on EPS620: Input Parameter:
621: + eps - the eigensolver context
622: . i - the solution index
623: - type - the type of error to compute
625: Output Parameter:
626: . error - the error
628: Notes:
629: The error can be computed in various ways, all of them based on the residual
630: norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
632: Level: beginner
634: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
635: @*/
636: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)637: {
639: Mat A,B;
640: Vec xr,xi,w[3];
641: PetscReal t,vecnorm=1.0;
642: PetscScalar kr,ki;
643: PetscBool flg;
650: EPSCheckSolved(eps,1);
652: /* allocate work vectors */
653: #if defined(PETSC_USE_COMPLEX)
654: EPSSetWorkVecs(eps,3);
655: xi = NULL;
656: w[1] = NULL;
657: #else
658: EPSSetWorkVecs(eps,5);
659: xi = eps->work[3];
660: w[1] = eps->work[4];
661: #endif
662: xr = eps->work[0];
663: w[0] = eps->work[1];
664: w[2] = eps->work[2];
666: /* compute residual norms */
667: EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
668: EPSComputeResidualNorm_Private(eps,kr,ki,xr,xi,w,error);
670: /* compute 2-norm of eigenvector */
671: if (eps->problem_type==EPS_GHEP) {
672: VecNorm(xr,NORM_2,&vecnorm);
673: }
675: /* compute error */
676: switch (type) {
677: case EPS_ERROR_ABSOLUTE:
678: break;
679: case EPS_ERROR_RELATIVE:
680: *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
681: break;
682: case EPS_ERROR_BACKWARD:
683: /* initialization of matrix norms */
684: if (!eps->nrma) {
685: STGetMatrix(eps->st,0,&A);
686: MatHasOperation(A,MATOP_NORM,&flg);
687: if (!flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
688: MatNorm(A,NORM_INFINITY,&eps->nrma);
689: }
690: if (eps->isgeneralized) {
691: if (!eps->nrmb) {
692: STGetMatrix(eps->st,1,&B);
693: MatHasOperation(B,MATOP_NORM,&flg);
694: if (!flg) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
695: MatNorm(B,NORM_INFINITY,&eps->nrmb);
696: }
697: } else eps->nrmb = 1.0;
698: t = SlepcAbsEigenvalue(kr,ki);
699: *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
700: break;
701: default:702: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
703: }
704: return(0);
705: }
707: /*
708: EPSGetStartVector - Generate a suitable vector to be used as the starting vector
709: for the recurrence that builds the right subspace.
711: Collective on EPS and Vec
713: Input Parameters:
714: + eps - the eigensolver context
715: - i - iteration number
717: Output Parameters:
718: . breakdown - flag indicating that a breakdown has occurred
720: Notes:
721: The start vector is computed from another vector: for the first step (i=0),
722: the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
723: vector is created. Then this vector is forced to be in the range of OP (only
724: for generalized definite problems) and orthonormalized with respect to all
725: V-vectors up to i-1. The resulting vector is placed in V[i].
727: The flag breakdown is set to true if either i=0 and the vector belongs to the
728: deflation space, or i>0 and the vector is linearly dependent with respect
729: to the V-vectors.
730: */
731: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)732: {
734: PetscReal norm;
735: PetscBool lindep;
736: Vec w,z;
742: /* For the first step, use the first initial vector, otherwise a random one */
743: if (i>0 || eps->nini==0) {
744: BVSetRandomColumn(eps->V,i);
745: }
747: /* Force the vector to be in the range of OP for definite generalized problems */
748: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
749: BVCreateVec(eps->V,&w);
750: BVCopyVec(eps->V,i,w);
751: BVGetColumn(eps->V,i,&z);
752: STApply(eps->st,w,z);
753: BVRestoreColumn(eps->V,i,&z);
754: VecDestroy(&w);
755: }
757: /* Orthonormalize the vector with respect to previous vectors */
758: BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep);
759: if (breakdown) *breakdown = lindep;
760: else if (lindep || norm == 0.0) {
761: if (i==0) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Initial vector is zero or belongs to the deflation space");
762: else SETERRQ(PetscObjectComm((PetscObject)eps),1,"Unable to generate more start vectors");
763: }
764: BVScaleColumn(eps->V,i,1.0/norm);
765: return(0);
766: }