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 6695 : PetscErrorCode EPSComputeVectors(EPS eps)
19 : {
20 6695 : PetscFunctionBegin;
21 6695 : EPSCheckSolved(eps,1);
22 6695 : if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
23 6695 : eps->state = EPS_STATE_EIGENVECTORS;
24 6695 : PetscFunctionReturn(PETSC_SUCCESS);
25 : }
26 :
27 : #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
28 :
29 855 : static PetscErrorCode EPSComputeValues(EPS eps)
30 : {
31 855 : PetscBool injective,iscomp,isfilter;
32 855 : PetscInt i,n,aux,nconv0;
33 855 : Mat A,B=NULL,G,Z;
34 :
35 855 : PetscFunctionBegin;
36 855 : switch (eps->categ) {
37 639 : case EPS_CATEGORY_KRYLOV:
38 : case EPS_CATEGORY_OTHER:
39 639 : PetscCall(STIsInjective(eps->st,&injective));
40 639 : if (injective) {
41 : /* one-to-one mapping: backtransform eigenvalues */
42 629 : 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 855 : 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 855 : PetscErrorCode EPSSolve(EPS eps)
122 : {
123 855 : PetscInt i;
124 855 : PetscBool hasname;
125 855 : STMatMode matmode;
126 855 : Mat A,B;
127 :
128 855 : PetscFunctionBegin;
129 855 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
130 855 : if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
131 855 : PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
132 :
133 : /* Call setup */
134 855 : PetscCall(EPSSetUp(eps));
135 855 : eps->nconv = 0;
136 855 : eps->its = 0;
137 27537 : for (i=0;i<eps->ncv;i++) {
138 26682 : eps->eigr[i] = 0.0;
139 26682 : eps->eigi[i] = 0.0;
140 26682 : eps->errest[i] = 0.0;
141 26682 : eps->perm[i] = i;
142 : }
143 855 : PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
144 855 : PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
145 :
146 : /* Call solver */
147 855 : PetscUseTypeMethod(eps,solve);
148 855 : PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
149 855 : eps->state = EPS_STATE_SOLVED;
150 :
151 : /* Only the first nconv columns contain useful information (except in CISS) */
152 855 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
153 855 : if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
154 :
155 : /* If inplace, purify eigenvectors before reverting operator */
156 855 : PetscCall(STGetMatMode(eps->st,&matmode));
157 855 : if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
158 855 : PetscCall(STPostSolve(eps->st));
159 :
160 : /* Map eigenvalues back to the original problem if appropriate */
161 855 : PetscCall(EPSComputeValues(eps));
162 :
163 : #if !defined(PETSC_USE_COMPLEX)
164 : /* Reorder conjugate eigenvalues (positive imaginary first) */
165 9087 : for (i=0;i<eps->nconv-1;i++) {
166 8232 : if (eps->eigi[i] != 0) {
167 121 : if (eps->eigi[i] < 0) {
168 21 : eps->eigi[i] = -eps->eigi[i];
169 21 : eps->eigi[i+1] = -eps->eigi[i+1];
170 : /* the next correction only works with eigenvectors */
171 21 : PetscCall(EPSComputeVectors(eps));
172 21 : PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
173 : }
174 121 : i++;
175 : }
176 : }
177 : #endif
178 :
179 : /* Sort eigenvalues according to eps->which parameter */
180 855 : PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
181 855 : PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
182 :
183 : /* Various viewers */
184 855 : PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
185 855 : PetscCall(EPSConvergedReasonViewFromOptions(eps));
186 855 : PetscCall(EPSErrorViewFromOptions(eps));
187 855 : PetscCall(EPSValuesViewFromOptions(eps));
188 855 : PetscCall(EPSVectorsViewFromOptions(eps));
189 :
190 855 : PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
191 855 : if (hasname) {
192 0 : PetscCall(EPSGetOperators(eps,&A,NULL));
193 0 : PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
194 : }
195 855 : if (eps->isgeneralized) {
196 210 : PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
197 210 : 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 855 : if (eps->nds) {
205 14 : PetscCall(BVSetNumConstraints(eps->V,0));
206 14 : eps->nds = 0;
207 : }
208 855 : eps->nini = 0;
209 855 : 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 146 : PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
237 : {
238 146 : PetscFunctionBegin;
239 146 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
240 146 : PetscAssertPointer(its,2);
241 146 : *its = eps->its;
242 146 : 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 501 : PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
264 : {
265 501 : PetscFunctionBegin;
266 501 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
267 501 : PetscAssertPointer(nconv,2);
268 501 : EPSCheckSolved(eps,1);
269 501 : *nconv = eps->nconv;
270 501 : 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 143 : PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
303 : {
304 143 : PetscFunctionBegin;
305 143 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
306 143 : PetscAssertPointer(reason,2);
307 143 : EPSCheckSolved(eps,1);
308 143 : *reason = eps->reason;
309 143 : 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 5 : for (i=0;i<eps->nconv;i++) {
358 4 : PetscCall(BVGetColumn(V,i,&w));
359 4 : PetscCall(VecPointwiseDivide(w,w,eps->D));
360 4 : PetscCall(BVRestoreColumn(V,i,&w));
361 : }
362 1 : PetscCall(BVOrthogonalize(V,NULL));
363 : }
364 60 : 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 5100 : PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
409 : {
410 5100 : PetscFunctionBegin;
411 5100 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
412 20400 : PetscValidLogicalCollectiveInt(eps,i,2);
413 5100 : EPSCheckSolved(eps,1);
414 5100 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
415 5100 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
416 5100 : PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
417 5100 : if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
418 5100 : 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 6705 : PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
448 : {
449 6705 : PetscInt k;
450 :
451 6705 : PetscFunctionBegin;
452 6705 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
453 6705 : EPSCheckSolved(eps,1);
454 6705 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
455 6705 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
456 6705 : k = eps->perm[i];
457 : #if defined(PETSC_USE_COMPLEX)
458 : if (eigr) *eigr = eps->eigr[k];
459 : if (eigi) *eigi = 0;
460 : #else
461 6705 : if (eigr) *eigr = eps->eigr[k];
462 6705 : if (eigi) *eigi = eps->eigi[k];
463 : #endif
464 6705 : 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 6570 : PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
502 : {
503 6570 : PetscInt k;
504 :
505 6570 : PetscFunctionBegin;
506 6570 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
507 26280 : PetscValidLogicalCollectiveInt(eps,i,2);
508 6570 : if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Vr,3); }
509 6570 : if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Vi,4); }
510 6570 : EPSCheckSolved(eps,1);
511 6570 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
512 6570 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
513 6570 : PetscCall(EPSComputeVectors(eps));
514 6570 : k = eps->perm[i];
515 6570 : PetscCall(BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi));
516 6570 : 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 72 : PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
553 : {
554 72 : PetscInt k;
555 :
556 72 : PetscFunctionBegin;
557 72 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
558 288 : PetscValidLogicalCollectiveInt(eps,i,2);
559 72 : if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Wr,3); }
560 72 : if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Wi,4); }
561 72 : EPSCheckSolved(eps,1);
562 72 : PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
563 72 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
564 72 : PetscCheck(i<eps->nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
565 72 : PetscCall(EPSComputeVectors(eps));
566 72 : k = eps->perm[i];
567 72 : PetscCall(BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi));
568 72 : 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 3960 : PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
616 : {
617 3960 : PetscInt nmat;
618 3960 : Mat A,B;
619 3960 : Vec u,w;
620 3960 : PetscScalar alpha;
621 : #if !defined(PETSC_USE_COMPLEX)
622 3960 : Vec v;
623 3960 : PetscReal ni,nr;
624 : #endif
625 3960 : PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
626 :
627 3960 : PetscFunctionBegin;
628 3960 : u = z[0]; w = z[2];
629 3960 : PetscCall(STGetNumMatrices(eps->st,&nmat));
630 3960 : PetscCall(STGetMatrix(eps->st,0,&A));
631 3960 : if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
632 :
633 : #if !defined(PETSC_USE_COMPLEX)
634 3960 : v = z[1];
635 3960 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
636 : #endif
637 3830 : PetscCall((*matmult)(A,xr,u)); /* u=A*x */
638 3830 : if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
639 3830 : if (nmat>1) PetscCall((*matmult)(B,xr,w));
640 1725 : else PetscCall(VecCopy(xr,w)); /* w=B*x */
641 3830 : alpha = trans? -PetscConj(kr): -kr;
642 3830 : PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
643 : }
644 3830 : PetscCall(VecNorm(u,NORM_2,norm));
645 : #if !defined(PETSC_USE_COMPLEX)
646 : } else {
647 130 : PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
648 130 : if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
649 130 : if (nmat>1) PetscCall((*matmult)(B,xr,v));
650 126 : else PetscCall(VecCopy(xr,v)); /* v=B*xr */
651 130 : PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
652 130 : if (nmat>1) PetscCall((*matmult)(B,xi,w));
653 126 : else PetscCall(VecCopy(xi,w)); /* w=B*xi */
654 130 : PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
655 : }
656 130 : PetscCall(VecNorm(u,NORM_2,&nr));
657 130 : PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
658 130 : if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
659 130 : PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
660 130 : PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
661 : }
662 130 : PetscCall(VecNorm(u,NORM_2,&ni));
663 130 : *norm = SlepcAbsEigenvalue(nr,ni);
664 : }
665 : #endif
666 3960 : 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 3663 : PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
692 : {
693 3663 : Mat A,B;
694 3663 : Vec xr,xi,w[3];
695 3663 : PetscReal t,vecnorm=1.0,errorl;
696 3663 : PetscScalar kr,ki;
697 3663 : PetscBool flg;
698 :
699 3663 : PetscFunctionBegin;
700 3663 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
701 14652 : PetscValidLogicalCollectiveInt(eps,i,2);
702 14652 : PetscValidLogicalCollectiveEnum(eps,type,3);
703 3663 : PetscAssertPointer(error,4);
704 3663 : EPSCheckSolved(eps,1);
705 :
706 : /* allocate work vectors */
707 : #if defined(PETSC_USE_COMPLEX)
708 : PetscCall(EPSSetWorkVecs(eps,3));
709 : xi = NULL;
710 : w[1] = NULL;
711 : #else
712 3663 : PetscCall(EPSSetWorkVecs(eps,5));
713 3663 : xi = eps->work[3];
714 3663 : w[1] = eps->work[4];
715 : #endif
716 3663 : xr = eps->work[0];
717 3663 : w[0] = eps->work[1];
718 3663 : w[2] = eps->work[2];
719 :
720 : /* compute residual norm */
721 3663 : PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
722 3663 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
723 :
724 : /* compute 2-norm of eigenvector */
725 3663 : 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 3663 : 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 51 : *error = PetscMax(*error,errorl);
732 : }
733 :
734 : /* compute error */
735 3663 : switch (type) {
736 : case EPS_ERROR_ABSOLUTE:
737 : break;
738 1991 : case EPS_ERROR_RELATIVE:
739 1991 : *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
740 1991 : break;
741 1638 : case EPS_ERROR_BACKWARD:
742 : /* initialization of matrix norms */
743 1638 : if (!eps->nrma) {
744 34 : PetscCall(STGetMatrix(eps->st,0,&A));
745 34 : PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
746 34 : PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
747 34 : PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
748 : }
749 1638 : if (eps->isgeneralized) {
750 1632 : if (!eps->nrmb) {
751 33 : PetscCall(STGetMatrix(eps->st,1,&B));
752 33 : PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
753 33 : PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
754 33 : PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
755 : }
756 6 : } else eps->nrmb = 1.0;
757 1638 : t = SlepcAbsEigenvalue(kr,ki);
758 1638 : *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
759 1638 : break;
760 0 : default:
761 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
762 : }
763 3663 : 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 636 : PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
791 : {
792 636 : PetscReal norm;
793 636 : PetscBool lindep;
794 636 : Vec w,z;
795 :
796 636 : PetscFunctionBegin;
797 636 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
798 2544 : PetscValidLogicalCollectiveInt(eps,i,2);
799 :
800 : /* For the first step, use the first initial vector, otherwise a random one */
801 636 : 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 636 : if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
805 112 : PetscCall(BVCreateVec(eps->V,&w));
806 112 : PetscCall(BVCopyVec(eps->V,i,w));
807 112 : PetscCall(BVGetColumn(eps->V,i,&z));
808 112 : PetscCall(STApply(eps->st,w,z));
809 112 : PetscCall(BVRestoreColumn(eps->V,i,&z));
810 112 : PetscCall(VecDestroy(&w));
811 : }
812 :
813 : /* Orthonormalize the vector with respect to previous vectors */
814 636 : PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
815 636 : if (breakdown) *breakdown = lindep;
816 599 : 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 636 : PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
821 636 : 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 24 : PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
829 : {
830 24 : PetscReal norm;
831 24 : PetscBool lindep;
832 :
833 24 : PetscFunctionBegin;
834 24 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
835 96 : PetscValidLogicalCollectiveInt(eps,i,2);
836 :
837 : /* For the first step, use the first initial vector, otherwise a random one */
838 24 : if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
839 :
840 : /* Orthonormalize the vector with respect to previous vectors */
841 24 : PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
842 24 : if (breakdown) *breakdown = lindep;
843 18 : 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 24 : PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
848 24 : PetscFunctionReturn(PETSC_SUCCESS);
849 : }
|