Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 : #include <slepc/private/bvimpl.h>
16 : #include <petscdraw.h>
17 :
18 6195 : PetscErrorCode EPSComputeVectors(EPS eps)
19 : {
20 6195 : PetscFunctionBegin;
21 6195 : EPSCheckSolved(eps,1);
22 6195 : if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
23 6195 : eps->state = EPS_STATE_EIGENVECTORS;
24 6195 : PetscFunctionReturn(PETSC_SUCCESS);
25 : }
26 :
27 : #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
28 :
29 830 : static PetscErrorCode EPSComputeValues(EPS eps)
30 : {
31 830 : PetscBool injective,iscomp,isfilter;
32 830 : PetscInt i,n,aux,nconv0;
33 830 : Mat A,B=NULL,G,Z;
34 :
35 830 : PetscFunctionBegin;
36 830 : switch (eps->categ) {
37 608 : case EPS_CATEGORY_KRYLOV:
38 : case EPS_CATEGORY_OTHER:
39 608 : PetscCall(STIsInjective(eps->st,&injective));
40 608 : if (injective) {
41 : /* one-to-one mapping: backtransform eigenvalues */
42 598 : PetscUseTypeMethod(eps,backtransform);
43 : } else {
44 : /* compute eigenvalues from Rayleigh quotient */
45 10 : PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
46 10 : if (!n) break;
47 10 : PetscCall(EPSGetOperators(eps,&A,&B));
48 10 : PetscCall(BVSetActiveColumns(eps->V,0,n));
49 10 : PetscCall(DSGetCompact(eps->ds,&iscomp));
50 10 : PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
51 10 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
52 10 : PetscCall(BVMatProject(eps->V,A,eps->V,G));
53 10 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
54 10 : if (B) {
55 0 : PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
56 0 : PetscCall(BVMatProject(eps->V,B,eps->V,G));
57 0 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
58 : }
59 10 : PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
60 10 : PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
61 10 : PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
62 10 : PetscCall(DSSetCompact(eps->ds,iscomp));
63 10 : if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
64 10 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
65 10 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
66 10 : PetscCall(BVMultInPlace(eps->V,Z,0,n));
67 10 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
68 : }
69 : /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
70 10 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
71 10 : if (isfilter) {
72 9 : nconv0 = eps->nconv;
73 82 : for (i=0;i<eps->nconv;i++) {
74 73 : if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
75 0 : eps->nconv--;
76 0 : if (i<eps->nconv) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
77 : }
78 : }
79 9 : 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 830 : PetscFunctionReturn(PETSC_SUCCESS);
89 : }
90 :
91 : /*@
92 : EPSSolve - Solves the eigensystem.
93 :
94 : Collective
95 :
96 : Input Parameter:
97 : . eps - eigensolver context obtained from EPSCreate()
98 :
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
109 :
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.
116 :
117 : Level: beginner
118 :
119 : .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
120 : @*/
121 830 : PetscErrorCode EPSSolve(EPS eps)
122 : {
123 830 : PetscInt i;
124 830 : PetscBool hasname;
125 830 : STMatMode matmode;
126 830 : Mat A,B;
127 :
128 830 : PetscFunctionBegin;
129 830 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
130 830 : if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
131 830 : PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
132 :
133 : /* Call setup */
134 830 : PetscCall(EPSSetUp(eps));
135 830 : eps->nconv = 0;
136 830 : eps->its = 0;
137 29289 : for (i=0;i<eps->ncv;i++) {
138 28459 : eps->eigr[i] = 0.0;
139 28459 : eps->eigi[i] = 0.0;
140 28459 : eps->errest[i] = 0.0;
141 28459 : eps->perm[i] = i;
142 : }
143 830 : PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
144 830 : PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
145 :
146 : /* Call solver */
147 830 : PetscUseTypeMethod(eps,solve);
148 830 : PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
149 830 : eps->state = EPS_STATE_SOLVED;
150 :
151 : /* Only the first nconv columns contain useful information (except in CISS) */
152 830 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
153 830 : if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
154 :
155 : /* If inplace, purify eigenvectors before reverting operator */
156 830 : PetscCall(STGetMatMode(eps->st,&matmode));
157 830 : if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
158 830 : PetscCall(STPostSolve(eps->st));
159 :
160 : /* Map eigenvalues back to the original problem if appropriate */
161 830 : PetscCall(EPSComputeValues(eps));
162 :
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
178 :
179 : /* Sort eigenvalues according to eps->which parameter */
180 830 : PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
181 830 : PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
182 :
183 : /* Various viewers */
184 830 : PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
185 830 : PetscCall(EPSConvergedReasonViewFromOptions(eps));
186 830 : PetscCall(EPSErrorViewFromOptions(eps));
187 830 : PetscCall(EPSValuesViewFromOptions(eps));
188 830 : PetscCall(EPSVectorsViewFromOptions(eps));
189 :
190 830 : PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
191 830 : if (hasname) {
192 0 : PetscCall(EPSGetOperators(eps,&A,NULL));
193 0 : PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
194 : }
195 830 : if (eps->isgeneralized) {
196 181 : PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
197 181 : if (hasname) {
198 0 : PetscCall(EPSGetOperators(eps,NULL,&B));
199 0 : PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
200 : }
201 : }
202 :
203 : /* Remove deflation and initial subspaces */
204 830 : if (eps->nds) {
205 14 : PetscCall(BVSetNumConstraints(eps->V,0));
206 14 : eps->nds = 0;
207 : }
208 830 : eps->nini = 0;
209 830 : PetscFunctionReturn(PETSC_SUCCESS);
210 : }
211 :
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.
216 :
217 : Not Collective
218 :
219 : Input Parameter:
220 : . eps - the eigensolver context
221 :
222 : Output Parameter:
223 : . its - number of iterations
224 :
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.
231 :
232 : Level: intermediate
233 :
234 : .seealso: EPSGetConvergedReason(), EPSSetTolerances()
235 : @*/
236 134 : PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
237 : {
238 134 : PetscFunctionBegin;
239 134 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
240 134 : PetscAssertPointer(its,2);
241 134 : *its = eps->its;
242 134 : PetscFunctionReturn(PETSC_SUCCESS);
243 : }
244 :
245 : /*@
246 : EPSGetConverged - Gets the number of converged eigenpairs.
247 :
248 : Not Collective
249 :
250 : Input Parameter:
251 : . eps - the eigensolver context
252 :
253 : Output Parameter:
254 : . nconv - number of converged eigenpairs
255 :
256 : Note:
257 : This function should be called after EPSSolve() has finished.
258 :
259 : Level: beginner
260 :
261 : .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
262 : @*/
263 479 : PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
264 : {
265 479 : PetscFunctionBegin;
266 479 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
267 479 : PetscAssertPointer(nconv,2);
268 479 : EPSCheckSolved(eps,1);
269 479 : *nconv = eps->nconv;
270 479 : PetscFunctionReturn(PETSC_SUCCESS);
271 : }
272 :
273 : /*@
274 : EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
275 : stopped.
276 :
277 : Not Collective
278 :
279 : Input Parameter:
280 : . eps - the eigensolver context
281 :
282 : Output Parameter:
283 : . reason - negative value indicates diverged, positive value converged
284 :
285 : Options Database Key:
286 : . -eps_converged_reason - print the reason to a viewer
287 :
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
295 :
296 : Can only be called after the call to EPSSolve() is complete.
297 :
298 : Level: intermediate
299 :
300 : .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
301 : @*/
302 132 : PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
303 : {
304 132 : PetscFunctionBegin;
305 132 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
306 132 : PetscAssertPointer(reason,2);
307 132 : EPSCheckSolved(eps,1);
308 132 : *reason = eps->reason;
309 132 : PetscFunctionReturn(PETSC_SUCCESS);
310 : }
311 :
312 : /*@
313 : EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
314 : subspace.
315 :
316 : Collective
317 :
318 : Input Parameter:
319 : . eps - the eigensolver context
320 :
321 : Output Parameter:
322 : . v - an array of vectors
323 :
324 : Notes:
325 : This function should be called after EPSSolve() has finished.
326 :
327 : The user should provide in v an array of nconv vectors, where nconv is
328 : the value returned by EPSGetConverged().
329 :
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).
336 :
337 : Level: intermediate
338 :
339 : .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
340 : @*/
341 8 : PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
342 : {
343 8 : PetscInt i;
344 8 : BV V=eps->V;
345 8 : Vec w;
346 :
347 8 : PetscFunctionBegin;
348 8 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
349 8 : PetscAssertPointer(v,2);
350 8 : PetscValidHeaderSpecific(*v,VEC_CLASSID,2);
351 8 : EPSCheckSolved(eps,1);
352 8 : 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 8 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
354 1 : PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
355 1 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
356 1 : PetscCall(BVCopy(eps->V,V));
357 6 : for (i=0;i<eps->nconv;i++) {
358 5 : PetscCall(BVGetColumn(V,i,&w));
359 5 : PetscCall(VecPointwiseDivide(w,w,eps->D));
360 5 : PetscCall(BVRestoreColumn(V,i,&w));
361 : }
362 1 : PetscCall(BVOrthogonalize(V,NULL));
363 : }
364 65 : for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
365 8 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
366 8 : PetscFunctionReturn(PETSC_SUCCESS);
367 : }
368 :
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.
372 :
373 : Collective
374 :
375 : Input Parameters:
376 : + eps - eigensolver context
377 : - i - index of the solution
378 :
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
384 :
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().
389 :
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.
394 :
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().
398 :
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.
402 :
403 : Level: beginner
404 :
405 : .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
406 : EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
407 : @*/
408 4961 : PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
409 : {
410 4961 : PetscFunctionBegin;
411 4961 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
412 19844 : PetscValidLogicalCollectiveInt(eps,i,2);
413 4961 : EPSCheckSolved(eps,1);
414 4961 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
415 4961 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
416 4961 : PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
417 4961 : if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
418 4961 : PetscFunctionReturn(PETSC_SUCCESS);
419 : }
420 :
421 : /*@C
422 : EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
423 :
424 : Not Collective
425 :
426 : Input Parameters:
427 : + eps - eigensolver context
428 : - i - index of the solution
429 :
430 : Output Parameters:
431 : + eigr - real part of eigenvalue
432 : - eigi - imaginary part of eigenvalue
433 :
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).
438 :
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().
442 :
443 : Level: beginner
444 :
445 : .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
446 : @*/
447 6173 : PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
448 : {
449 6173 : PetscInt k;
450 :
451 6173 : PetscFunctionBegin;
452 6173 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
453 6173 : EPSCheckSolved(eps,1);
454 6173 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
455 6173 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
456 6173 : k = eps->perm[i];
457 : #if defined(PETSC_USE_COMPLEX)
458 6173 : if (eigr) *eigr = eps->eigr[k];
459 6173 : if (eigi) *eigi = 0;
460 : #else
461 : if (eigr) *eigr = eps->eigr[k];
462 : if (eigi) *eigi = eps->eigi[k];
463 : #endif
464 6173 : PetscFunctionReturn(PETSC_SUCCESS);
465 : }
466 :
467 : /*@
468 : EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
469 :
470 : Collective
471 :
472 : Input Parameters:
473 : + eps - eigensolver context
474 : - i - index of the solution
475 :
476 : Output Parameters:
477 : + Vr - real part of eigenvector
478 : - Vi - imaginary part of eigenvector
479 :
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().
483 :
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.
488 :
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().
492 :
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.
496 :
497 : Level: beginner
498 :
499 : .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
500 : @*/
501 6090 : PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
502 : {
503 6090 : PetscInt k;
504 :
505 6090 : PetscFunctionBegin;
506 6090 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
507 24360 : PetscValidLogicalCollectiveInt(eps,i,2);
508 6090 : if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Vr,3); }
509 6090 : if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Vi,4); }
510 6090 : EPSCheckSolved(eps,1);
511 6090 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
512 6090 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
513 6090 : PetscCall(EPSComputeVectors(eps));
514 6090 : k = eps->perm[i];
515 6090 : PetscCall(BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi));
516 6090 : PetscFunctionReturn(PETSC_SUCCESS);
517 : }
518 :
519 : /*@
520 : EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
521 :
522 : Collective
523 :
524 : Input Parameters:
525 : + eps - eigensolver context
526 : - i - index of the solution
527 :
528 : Output Parameters:
529 : + Wr - real part of left eigenvector
530 : - Wi - imaginary part of left eigenvector
531 :
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().
535 :
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.
540 :
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().
544 :
545 : Left eigenvectors are available only if the twosided flag was set, see
546 : EPSSetTwoSided().
547 :
548 : Level: intermediate
549 :
550 : .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
551 : @*/
552 71 : PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
553 : {
554 71 : PetscInt k;
555 :
556 71 : PetscFunctionBegin;
557 71 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
558 284 : PetscValidLogicalCollectiveInt(eps,i,2);
559 71 : if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Wr,3); }
560 71 : if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Wi,4); }
561 71 : EPSCheckSolved(eps,1);
562 71 : PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
563 71 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
564 71 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
565 71 : PetscCall(EPSComputeVectors(eps));
566 71 : k = eps->perm[i];
567 71 : PetscCall(BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi));
568 71 : PetscFunctionReturn(PETSC_SUCCESS);
569 : }
570 :
571 : /*@
572 : EPSGetErrorEstimate - Returns the error estimate associated to the i-th
573 : computed eigenpair.
574 :
575 : Not Collective
576 :
577 : Input Parameters:
578 : + eps - eigensolver context
579 : - i - index of eigenpair
580 :
581 : Output Parameter:
582 : . errest - the error estimate
583 :
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.
588 :
589 : Level: advanced
590 :
591 : .seealso: EPSComputeError()
592 : @*/
593 46 : PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
594 : {
595 46 : PetscFunctionBegin;
596 46 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
597 46 : PetscAssertPointer(errest,3);
598 46 : EPSCheckSolved(eps,1);
599 46 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
600 46 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
601 46 : *errest = eps->errest[eps->perm[i]];
602 46 : PetscFunctionReturn(PETSC_SUCCESS);
603 : }
604 :
605 : /*
606 : EPSComputeResidualNorm_Private - Computes the norm of the residual vector
607 : associated with an eigenpair.
608 :
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 3964 : PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
616 : {
617 3964 : PetscInt nmat;
618 3964 : Mat A,B;
619 3964 : Vec u,w;
620 3964 : PetscScalar alpha;
621 : #if !defined(PETSC_USE_COMPLEX)
622 : Vec v;
623 : PetscReal ni,nr;
624 : #endif
625 3964 : PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
626 :
627 3964 : PetscFunctionBegin;
628 3964 : u = z[0]; w = z[2];
629 3964 : PetscCall(STGetNumMatrices(eps->st,&nmat));
630 3964 : PetscCall(STGetMatrix(eps->st,0,&A));
631 3964 : if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
632 :
633 : #if !defined(PETSC_USE_COMPLEX)
634 : v = z[1];
635 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
636 : #endif
637 3964 : PetscCall((*matmult)(A,xr,u)); /* u=A*x */
638 3964 : if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
639 3964 : if (nmat>1) PetscCall((*matmult)(B,xr,w));
640 1863 : else PetscCall(VecCopy(xr,w)); /* w=B*x */
641 3964 : alpha = trans? -PetscConj(kr): -kr;
642 3964 : PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
643 : }
644 3964 : 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 3964 : PetscFunctionReturn(PETSC_SUCCESS);
667 : }
668 :
669 : /*@
670 : EPSComputeError - Computes the error (based on the residual norm) associated
671 : with the i-th computed eigenpair.
672 :
673 : Collective
674 :
675 : Input Parameters:
676 : + eps - the eigensolver context
677 : . i - the solution index
678 : - type - the type of error to compute
679 :
680 : Output Parameter:
681 : . error - the error
682 :
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.
686 :
687 : Level: beginner
688 :
689 : .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
690 : @*/
691 3576 : PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
692 : {
693 3576 : Mat A,B;
694 3576 : Vec xr,xi,w[3];
695 3576 : PetscReal t,vecnorm=1.0,errorl;
696 3576 : PetscScalar kr,ki;
697 3576 : PetscBool flg;
698 :
699 3576 : PetscFunctionBegin;
700 3576 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
701 14304 : PetscValidLogicalCollectiveInt(eps,i,2);
702 14304 : PetscValidLogicalCollectiveEnum(eps,type,3);
703 3576 : PetscAssertPointer(error,4);
704 3576 : EPSCheckSolved(eps,1);
705 :
706 : /* allocate work vectors */
707 : #if defined(PETSC_USE_COMPLEX)
708 3576 : PetscCall(EPSSetWorkVecs(eps,3));
709 3576 : xi = NULL;
710 3576 : w[1] = NULL;
711 : #else
712 : PetscCall(EPSSetWorkVecs(eps,5));
713 : xi = eps->work[3];
714 : w[1] = eps->work[4];
715 : #endif
716 3576 : xr = eps->work[0];
717 3576 : w[0] = eps->work[1];
718 3576 : w[2] = eps->work[2];
719 :
720 : /* compute residual norm */
721 3576 : PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
722 3576 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
723 :
724 : /* compute 2-norm of eigenvector */
725 3576 : if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
726 :
727 : /* if two-sided, compute left residual norm and take the maximum */
728 3576 : if (eps->twosided) {
729 36 : PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
730 36 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
731 50 : *error = PetscMax(*error,errorl);
732 : }
733 :
734 : /* compute error */
735 3576 : switch (type) {
736 : case EPS_ERROR_ABSOLUTE:
737 : break;
738 1872 : case EPS_ERROR_RELATIVE:
739 1872 : *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
740 1872 : break;
741 1674 : case EPS_ERROR_BACKWARD:
742 : /* initialization of matrix norms */
743 1674 : if (!eps->nrma) {
744 43 : PetscCall(STGetMatrix(eps->st,0,&A));
745 43 : PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
746 43 : PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
747 43 : PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
748 : }
749 1674 : if (eps->isgeneralized) {
750 1629 : if (!eps->nrmb) {
751 32 : PetscCall(STGetMatrix(eps->st,1,&B));
752 32 : PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
753 32 : PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
754 32 : PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
755 : }
756 45 : } else eps->nrmb = 1.0;
757 1674 : t = SlepcAbsEigenvalue(kr,ki);
758 1674 : *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
759 1674 : break;
760 0 : default:
761 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
762 : }
763 3576 : PetscFunctionReturn(PETSC_SUCCESS);
764 : }
765 :
766 : /*
767 : EPSGetStartVector - Generate a suitable vector to be used as the starting vector
768 : for the recurrence that builds the right subspace.
769 :
770 : Collective
771 :
772 : Input Parameters:
773 : + eps - the eigensolver context
774 : - i - iteration number
775 :
776 : Output Parameters:
777 : . breakdown - flag indicating that a breakdown has occurred
778 :
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].
785 :
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 614 : PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
791 : {
792 614 : PetscReal norm;
793 614 : PetscBool lindep;
794 614 : Vec w,z;
795 :
796 614 : PetscFunctionBegin;
797 614 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
798 2456 : PetscValidLogicalCollectiveInt(eps,i,2);
799 :
800 : /* For the first step, use the first initial vector, otherwise a random one */
801 614 : if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
802 :
803 : /* Force the vector to be in the range of OP for definite generalized problems */
804 614 : if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
805 92 : PetscCall(BVCreateVec(eps->V,&w));
806 92 : PetscCall(BVCopyVec(eps->V,i,w));
807 92 : PetscCall(BVGetColumn(eps->V,i,&z));
808 92 : PetscCall(STApply(eps->st,w,z));
809 92 : PetscCall(BVRestoreColumn(eps->V,i,&z));
810 92 : PetscCall(VecDestroy(&w));
811 : }
812 :
813 : /* Orthonormalize the vector with respect to previous vectors */
814 614 : PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
815 614 : if (breakdown) *breakdown = lindep;
816 578 : else if (lindep || norm == 0.0) {
817 0 : PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
818 0 : PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
819 : }
820 614 : PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
821 614 : PetscFunctionReturn(PETSC_SUCCESS);
822 : }
823 :
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 21 : PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
829 : {
830 21 : PetscReal norm;
831 21 : PetscBool lindep;
832 :
833 21 : PetscFunctionBegin;
834 21 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
835 84 : PetscValidLogicalCollectiveInt(eps,i,2);
836 :
837 : /* For the first step, use the first initial vector, otherwise a random one */
838 21 : if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
839 :
840 : /* Orthonormalize the vector with respect to previous vectors */
841 21 : PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
842 21 : if (breakdown) *breakdown = lindep;
843 16 : else if (lindep || norm == 0.0) {
844 0 : PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
845 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
846 : }
847 21 : PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
848 21 : PetscFunctionReturn(PETSC_SUCCESS);
849 : }
|