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 7073 : PetscErrorCode EPSComputeVectors(EPS eps)
19 : {
20 7073 : PetscFunctionBegin;
21 7073 : EPSCheckSolved(eps,1);
22 7073 : if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
23 7073 : eps->state = EPS_STATE_EIGENVECTORS;
24 7073 : PetscFunctionReturn(PETSC_SUCCESS);
25 : }
26 :
27 869 : static PetscErrorCode EPSComputeValues(EPS eps)
28 : {
29 869 : PetscBool injective,iscomp,isfilter;
30 869 : PetscInt i,n,aux,nconv0;
31 869 : Mat A,B=NULL,G,Z;
32 :
33 869 : PetscFunctionBegin;
34 869 : switch (eps->categ) {
35 653 : case EPS_CATEGORY_KRYLOV:
36 : case EPS_CATEGORY_OTHER:
37 653 : PetscCall(STIsInjective(eps->st,&injective));
38 653 : if (injective) {
39 : /* one-to-one mapping: backtransform eigenvalues */
40 643 : PetscUseTypeMethod(eps,backtransform);
41 : } else {
42 : /* compute eigenvalues from Rayleigh quotient */
43 10 : PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
44 10 : if (!n) break;
45 10 : PetscCall(EPSGetOperators(eps,&A,&B));
46 10 : PetscCall(BVSetActiveColumns(eps->V,0,n));
47 10 : PetscCall(DSGetCompact(eps->ds,&iscomp));
48 10 : PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
49 10 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
50 10 : PetscCall(BVMatProject(eps->V,A,eps->V,G));
51 10 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
52 10 : if (B) {
53 0 : PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
54 0 : PetscCall(BVMatProject(eps->V,B,eps->V,G));
55 0 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
56 : }
57 10 : PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
58 10 : PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
59 10 : PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
60 10 : PetscCall(DSSetCompact(eps->ds,iscomp));
61 10 : if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
62 10 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
63 10 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
64 10 : PetscCall(BVMultInPlace(eps->V,Z,0,n));
65 10 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
66 : }
67 : /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
68 10 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
69 10 : if (isfilter) {
70 9 : nconv0 = eps->nconv;
71 82 : for (i=0;i<eps->nconv;i++) {
72 73 : if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
73 0 : eps->nconv--;
74 0 : if (i<eps->nconv) { SlepcSwap(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
75 : }
76 : }
77 9 : if (nconv0>eps->nconv) PetscCall(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 869 : PetscFunctionReturn(PETSC_SUCCESS);
87 : }
88 :
89 : /*@
90 : EPSSolve - Solves the eigensystem.
91 :
92 : Collective
93 :
94 : Input Parameter:
95 : . eps - eigensolver context obtained from EPSCreate()
96 :
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
107 :
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.
114 :
115 : Level: beginner
116 :
117 : .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
118 : @*/
119 869 : PetscErrorCode EPSSolve(EPS eps)
120 : {
121 869 : PetscInt i;
122 869 : PetscBool hasname;
123 869 : STMatMode matmode;
124 869 : Mat A,B;
125 :
126 869 : PetscFunctionBegin;
127 869 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
128 869 : if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
129 869 : PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));
130 :
131 : /* Call setup */
132 869 : PetscCall(EPSSetUp(eps));
133 869 : eps->nconv = 0;
134 869 : eps->its = 0;
135 27775 : for (i=0;i<eps->ncv;i++) {
136 26906 : eps->eigr[i] = 0.0;
137 26906 : eps->eigi[i] = 0.0;
138 26906 : eps->errest[i] = 0.0;
139 26906 : eps->perm[i] = i;
140 : }
141 869 : PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
142 869 : PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));
143 :
144 : /* Call solver */
145 869 : PetscUseTypeMethod(eps,solve);
146 869 : PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
147 869 : eps->state = EPS_STATE_SOLVED;
148 :
149 : /* Only the first nconv columns contain useful information (except in CISS) */
150 869 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
151 869 : if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));
152 :
153 : /* If inplace, purify eigenvectors before reverting operator */
154 869 : PetscCall(STGetMatMode(eps->st,&matmode));
155 869 : if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
156 869 : PetscCall(STPostSolve(eps->st));
157 :
158 : /* Map eigenvalues back to the original problem if appropriate */
159 869 : PetscCall(EPSComputeValues(eps));
160 :
161 : #if !defined(PETSC_USE_COMPLEX)
162 : /* Reorder conjugate eigenvalues (positive imaginary first) */
163 9186 : for (i=0;i<eps->nconv-1;i++) {
164 8317 : if (eps->eigi[i] != 0) {
165 121 : if (eps->eigi[i] < 0) {
166 21 : eps->eigi[i] = -eps->eigi[i];
167 21 : eps->eigi[i+1] = -eps->eigi[i+1];
168 : /* the next correction only works with eigenvectors */
169 21 : PetscCall(EPSComputeVectors(eps));
170 21 : PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
171 : }
172 121 : i++;
173 : }
174 : }
175 : #endif
176 :
177 : /* Sort eigenvalues according to eps->which parameter */
178 869 : PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
179 869 : PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));
180 :
181 : /* Various viewers */
182 869 : PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
183 869 : PetscCall(EPSConvergedReasonViewFromOptions(eps));
184 869 : PetscCall(EPSErrorViewFromOptions(eps));
185 869 : PetscCall(EPSValuesViewFromOptions(eps));
186 869 : PetscCall(EPSVectorsViewFromOptions(eps));
187 :
188 869 : PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
189 869 : if (hasname) {
190 0 : PetscCall(EPSGetOperators(eps,&A,NULL));
191 0 : PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
192 : }
193 869 : if (eps->isgeneralized) {
194 210 : PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
195 210 : if (hasname) {
196 0 : PetscCall(EPSGetOperators(eps,NULL,&B));
197 0 : PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
198 : }
199 : }
200 :
201 : /* Remove deflation and initial subspaces */
202 869 : if (eps->nds) {
203 14 : PetscCall(BVSetNumConstraints(eps->V,0));
204 14 : eps->nds = 0;
205 : }
206 869 : eps->nini = 0;
207 869 : PetscFunctionReturn(PETSC_SUCCESS);
208 : }
209 :
210 : /*@
211 : EPSGetIterationNumber - Gets the current iteration number. If the
212 : call to EPSSolve() is complete, then it returns the number of iterations
213 : carried out by the solution method.
214 :
215 : Not Collective
216 :
217 : Input Parameter:
218 : . eps - the eigensolver context
219 :
220 : Output Parameter:
221 : . its - number of iterations
222 :
223 : Note:
224 : During the i-th iteration this call returns i-1. If EPSSolve() is
225 : complete, then parameter "its" contains either the iteration number at
226 : which convergence was successfully reached, or failure was detected.
227 : Call EPSGetConvergedReason() to determine if the solver converged or
228 : failed and why.
229 :
230 : Level: intermediate
231 :
232 : .seealso: EPSGetConvergedReason(), EPSSetTolerances()
233 : @*/
234 146 : PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
235 : {
236 146 : PetscFunctionBegin;
237 146 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
238 146 : PetscAssertPointer(its,2);
239 146 : *its = eps->its;
240 146 : PetscFunctionReturn(PETSC_SUCCESS);
241 : }
242 :
243 : /*@
244 : EPSGetConverged - Gets the number of converged eigenpairs.
245 :
246 : Not Collective
247 :
248 : Input Parameter:
249 : . eps - the eigensolver context
250 :
251 : Output Parameter:
252 : . nconv - number of converged eigenpairs
253 :
254 : Note:
255 : This function should be called after EPSSolve() has finished.
256 :
257 : Level: beginner
258 :
259 : .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
260 : @*/
261 515 : PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
262 : {
263 515 : PetscFunctionBegin;
264 515 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
265 515 : PetscAssertPointer(nconv,2);
266 515 : EPSCheckSolved(eps,1);
267 515 : PetscCall(EPS_GetActualConverged(eps,nconv));
268 515 : PetscFunctionReturn(PETSC_SUCCESS);
269 : }
270 :
271 : /*@
272 : EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
273 : stopped.
274 :
275 : Not Collective
276 :
277 : Input Parameter:
278 : . eps - the eigensolver context
279 :
280 : Output Parameter:
281 : . reason - negative value indicates diverged, positive value converged
282 :
283 : Options Database Key:
284 : . -eps_converged_reason - print the reason to a viewer
285 :
286 : Notes:
287 : Possible values for reason are
288 : + EPS_CONVERGED_TOL - converged up to tolerance
289 : . EPS_CONVERGED_USER - converged due to a user-defined condition
290 : . EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
291 : . EPS_DIVERGED_BREAKDOWN - generic breakdown in method
292 : - EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
293 :
294 : Can only be called after the call to EPSSolve() is complete.
295 :
296 : Level: intermediate
297 :
298 : .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
299 : @*/
300 143 : PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
301 : {
302 143 : PetscFunctionBegin;
303 143 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
304 143 : PetscAssertPointer(reason,2);
305 143 : EPSCheckSolved(eps,1);
306 143 : *reason = eps->reason;
307 143 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 :
310 : /*@
311 : EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
312 : subspace.
313 :
314 : Collective
315 :
316 : Input Parameter:
317 : . eps - the eigensolver context
318 :
319 : Output Parameter:
320 : . v - an array of vectors
321 :
322 : Notes:
323 : This function should be called after EPSSolve() has finished.
324 :
325 : The user should provide in v an array of nconv vectors, where nconv is
326 : the value returned by EPSGetConverged().
327 :
328 : The first k vectors returned in v span an invariant subspace associated
329 : with the first k computed eigenvalues (note that this is not true if the
330 : k-th eigenvalue is complex and matrix A is real; in this case the first
331 : k+1 vectors should be used). An invariant subspace X of A satisfies Ax
332 : in X for all x in X (a similar definition applies for generalized
333 : eigenproblems).
334 :
335 : Level: intermediate
336 :
337 : .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
338 : @*/
339 8 : PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
340 : {
341 8 : PetscInt i;
342 8 : BV V=eps->V;
343 8 : Vec w;
344 :
345 8 : PetscFunctionBegin;
346 8 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
347 8 : PetscAssertPointer(v,2);
348 8 : PetscValidHeaderSpecific(*v,VEC_CLASSID,2);
349 8 : EPSCheckSolved(eps,1);
350 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");
351 8 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
352 1 : PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
353 1 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
354 1 : PetscCall(BVCopy(eps->V,V));
355 5 : for (i=0;i<eps->nconv;i++) {
356 4 : PetscCall(BVGetColumn(V,i,&w));
357 4 : PetscCall(VecPointwiseDivide(w,w,eps->D));
358 4 : PetscCall(BVRestoreColumn(V,i,&w));
359 : }
360 1 : PetscCall(BVOrthogonalize(V,NULL));
361 : }
362 60 : for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
363 8 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
364 8 : PetscFunctionReturn(PETSC_SUCCESS);
365 : }
366 :
367 : /*@
368 : EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
369 : EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.
370 :
371 : Collective
372 :
373 : Input Parameters:
374 : + eps - eigensolver context
375 : - i - index of the solution
376 :
377 : Output Parameters:
378 : + eigr - real part of eigenvalue
379 : . eigi - imaginary part of eigenvalue
380 : . Vr - real part of eigenvector
381 : - Vi - imaginary part of eigenvector
382 :
383 : Notes:
384 : It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
385 : required. Otherwise, the caller must provide valid Vec objects, i.e.,
386 : they must be created by the calling program with e.g. MatCreateVecs().
387 :
388 : If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
389 : configured with complex scalars the eigenvalue is stored
390 : directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
391 : set to zero). In both cases, the user can pass NULL in eigi and Vi.
392 :
393 : The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
394 : Eigenpairs are indexed according to the ordering criterion established
395 : with EPSSetWhichEigenpairs().
396 :
397 : The 2-norm of the eigenvector is one unless the problem is generalized
398 : Hermitian. In this case the eigenvector is normalized with respect to the
399 : norm defined by the B matrix.
400 :
401 : Level: beginner
402 :
403 : .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
404 : EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
405 : @*/
406 5149 : PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
407 : {
408 5149 : PetscInt nconv;
409 :
410 5149 : PetscFunctionBegin;
411 5149 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
412 15447 : PetscValidLogicalCollectiveInt(eps,i,2);
413 5149 : EPSCheckSolved(eps,1);
414 5149 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
415 5149 : PetscCall(EPS_GetActualConverged(eps,&nconv));
416 5149 : PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
417 5149 : PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
418 5149 : if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
419 5149 : PetscFunctionReturn(PETSC_SUCCESS);
420 : }
421 :
422 : /*@
423 : EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().
424 :
425 : Not Collective
426 :
427 : Input Parameters:
428 : + eps - eigensolver context
429 : - i - index of the solution
430 :
431 : Output Parameters:
432 : + eigr - real part of eigenvalue
433 : - eigi - imaginary part of eigenvalue
434 :
435 : Notes:
436 : If the eigenvalue is real, then eigi is set to zero. If PETSc is
437 : configured with complex scalars the eigenvalue is stored
438 : directly in eigr (eigi is set to zero).
439 :
440 : The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
441 : Eigenpairs are indexed according to the ordering criterion established
442 : with EPSSetWhichEigenpairs().
443 :
444 : Level: beginner
445 :
446 : .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
447 : @*/
448 10481 : PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
449 : {
450 10481 : PetscInt k,nconv;
451 :
452 10481 : PetscFunctionBegin;
453 10481 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
454 10481 : EPSCheckSolved(eps,1);
455 10481 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
456 10481 : PetscCall(EPS_GetActualConverged(eps,&nconv));
457 10481 : PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
458 10481 : if (nconv==eps->nconv) {
459 10369 : k = eps->perm[i];
460 : #if defined(PETSC_USE_COMPLEX)
461 : if (eigr) *eigr = eps->eigr[k];
462 : if (eigi) *eigi = 0;
463 : #else
464 10369 : if (eigr) *eigr = eps->eigr[k];
465 10369 : if (eigi) *eigi = eps->eigi[k];
466 : #endif
467 : } else {
468 112 : PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
469 : /* BSE problem, even index is +lambda, odd index is -lambda */
470 112 : k = eps->perm[i/2];
471 : #if defined(PETSC_USE_COMPLEX)
472 : if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
473 : if (eigi) *eigi = 0;
474 : #else
475 112 : if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
476 112 : if (eigi) *eigi = eps->eigi[k];
477 : #endif
478 : }
479 10481 : PetscFunctionReturn(PETSC_SUCCESS);
480 : }
481 :
482 : /*@
483 : EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().
484 :
485 : Collective
486 :
487 : Input Parameters:
488 : + eps - eigensolver context
489 : - i - index of the solution
490 :
491 : Output Parameters:
492 : + Vr - real part of eigenvector
493 : - Vi - imaginary part of eigenvector
494 :
495 : Notes:
496 : The caller must provide valid Vec objects, i.e., they must be created
497 : by the calling program with e.g. MatCreateVecs().
498 :
499 : If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
500 : configured with complex scalars the eigenvector is stored
501 : directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
502 : or Vi if one of them is not required.
503 :
504 : The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
505 : Eigenpairs are indexed according to the ordering criterion established
506 : with EPSSetWhichEigenpairs().
507 :
508 : The 2-norm of the eigenvector is one unless the problem is generalized
509 : Hermitian. In this case the eigenvector is normalized with respect to the
510 : norm defined by the B matrix.
511 :
512 : Level: beginner
513 :
514 : .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
515 : @*/
516 6787 : PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
517 : {
518 6787 : PetscInt nconv;
519 :
520 6787 : PetscFunctionBegin;
521 6787 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
522 20361 : PetscValidLogicalCollectiveInt(eps,i,2);
523 6787 : if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Vr,3); }
524 6787 : if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Vi,4); }
525 6787 : EPSCheckSolved(eps,1);
526 6787 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
527 6787 : PetscCall(EPS_GetActualConverged(eps,&nconv));
528 6787 : PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
529 6787 : PetscCall(EPSComputeVectors(eps));
530 6787 : PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
531 6787 : PetscFunctionReturn(PETSC_SUCCESS);
532 : }
533 :
534 : /*@
535 : EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().
536 :
537 : Collective
538 :
539 : Input Parameters:
540 : + eps - eigensolver context
541 : - i - index of the solution
542 :
543 : Output Parameters:
544 : + Wr - real part of left eigenvector
545 : - Wi - imaginary part of left eigenvector
546 :
547 : Notes:
548 : The caller must provide valid Vec objects, i.e., they must be created
549 : by the calling program with e.g. MatCreateVecs().
550 :
551 : If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
552 : configured with complex scalars the eigenvector is stored directly in Wr
553 : (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
554 : one of them is not required.
555 :
556 : The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
557 : Eigensolutions are indexed according to the ordering criterion established
558 : with EPSSetWhichEigenpairs().
559 :
560 : Left eigenvectors are available only if the twosided flag was set, see
561 : EPSSetTwoSided().
562 :
563 : Level: intermediate
564 :
565 : .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
566 : @*/
567 234 : PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
568 : {
569 234 : PetscInt nconv;
570 234 : PetscBool trivial;
571 234 : Mat H;
572 234 : IS is[2];
573 234 : Vec v;
574 :
575 234 : PetscFunctionBegin;
576 234 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
577 702 : PetscValidLogicalCollectiveInt(eps,i,2);
578 234 : if (Wr) { PetscValidHeaderSpecific(Wr,VEC_CLASSID,3); PetscCheckSameComm(eps,1,Wr,3); }
579 234 : if (Wi) { PetscValidHeaderSpecific(Wi,VEC_CLASSID,4); PetscCheckSameComm(eps,1,Wi,4); }
580 234 : EPSCheckSolved(eps,1);
581 234 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
582 234 : PetscCall(EPS_GetActualConverged(eps,&nconv));
583 234 : PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
584 :
585 234 : trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE)? PETSC_TRUE: PETSC_FALSE;
586 234 : if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");
587 :
588 234 : PetscCall(EPSComputeVectors(eps));
589 234 : if (trivial) {
590 158 : PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
591 158 : if (eps->problem_type==EPS_BSE) { /* change sign of bottom part of the vector */
592 158 : PetscCall(STGetMatrix(eps->st,0,&H));
593 158 : PetscCall(MatNestGetISs(H,is,NULL));
594 158 : if (Wr) {
595 158 : PetscCall(VecGetSubVector(Wr,is[1],&v));
596 158 : PetscCall(VecScale(v,-1.0));
597 158 : PetscCall(VecRestoreSubVector(Wr,is[1],&v));
598 : }
599 : #if !defined(PETSC_USE_COMPLEX)
600 158 : if (Wi) {
601 0 : PetscCall(VecGetSubVector(Wi,is[1],&v));
602 0 : PetscCall(VecScale(v,-1.0));
603 0 : PetscCall(VecRestoreSubVector(Wi,is[1],&v));
604 : }
605 : #endif
606 : }
607 : } else {
608 76 : PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
609 : }
610 234 : PetscFunctionReturn(PETSC_SUCCESS);
611 : }
612 :
613 : /*@
614 : EPSGetErrorEstimate - Returns the error estimate associated to the i-th
615 : computed eigenpair.
616 :
617 : Not Collective
618 :
619 : Input Parameters:
620 : + eps - eigensolver context
621 : - i - index of eigenpair
622 :
623 : Output Parameter:
624 : . errest - the error estimate
625 :
626 : Notes:
627 : This is the error estimate used internally by the eigensolver. The actual
628 : error bound can be computed with EPSComputeError(). See also the users
629 : manual for details.
630 :
631 : Level: advanced
632 :
633 : .seealso: EPSComputeError()
634 : @*/
635 46 : PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
636 : {
637 46 : PetscInt nconv;
638 :
639 46 : PetscFunctionBegin;
640 46 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
641 46 : PetscAssertPointer(errest,3);
642 46 : EPSCheckSolved(eps,1);
643 46 : PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
644 46 : PetscCall(EPS_GetActualConverged(eps,&nconv));
645 46 : PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
646 46 : if (nconv==eps->nconv) {
647 46 : *errest = eps->errest[eps->perm[i]];
648 : } else {
649 0 : PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
650 : /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
651 0 : *errest = eps->errest[eps->perm[i/2]];
652 : }
653 46 : PetscFunctionReturn(PETSC_SUCCESS);
654 : }
655 :
656 : /*
657 : EPSComputeResidualNorm_Private - Computes the norm of the residual vector
658 : associated with an eigenpair.
659 :
660 : Input Parameters:
661 : trans - whether A' must be used instead of A
662 : kr,ki - eigenvalue
663 : xr,xi - eigenvector
664 : z - three work vectors (the second one not referenced in complex scalars)
665 : */
666 4016 : PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
667 : {
668 4016 : PetscInt nmat;
669 4016 : Mat A,B;
670 4016 : Vec u,w;
671 4016 : PetscScalar alpha;
672 : #if !defined(PETSC_USE_COMPLEX)
673 4016 : Vec v;
674 4016 : PetscReal ni,nr;
675 : #endif
676 4016 : PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
677 :
678 4016 : PetscFunctionBegin;
679 4016 : u = z[0]; w = z[2];
680 4016 : PetscCall(STGetNumMatrices(eps->st,&nmat));
681 4016 : PetscCall(STGetMatrix(eps->st,0,&A));
682 4016 : if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
683 :
684 : #if !defined(PETSC_USE_COMPLEX)
685 4016 : v = z[1];
686 4016 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
687 : #endif
688 3886 : PetscCall((*matmult)(A,xr,u)); /* u=A*x */
689 3886 : if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
690 3886 : if (nmat>1) PetscCall((*matmult)(B,xr,w));
691 1781 : else PetscCall(VecCopy(xr,w)); /* w=B*x */
692 3886 : alpha = trans? -PetscConj(kr): -kr;
693 3886 : PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
694 : }
695 3886 : PetscCall(VecNorm(u,NORM_2,norm));
696 : #if !defined(PETSC_USE_COMPLEX)
697 : } else {
698 130 : PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
699 130 : if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
700 130 : if (nmat>1) PetscCall((*matmult)(B,xr,v));
701 126 : else PetscCall(VecCopy(xr,v)); /* v=B*xr */
702 130 : PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
703 130 : if (nmat>1) PetscCall((*matmult)(B,xi,w));
704 126 : else PetscCall(VecCopy(xi,w)); /* w=B*xi */
705 130 : PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
706 : }
707 130 : PetscCall(VecNorm(u,NORM_2,&nr));
708 130 : PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
709 130 : if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
710 130 : PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
711 130 : PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
712 : }
713 130 : PetscCall(VecNorm(u,NORM_2,&ni));
714 130 : *norm = SlepcAbsEigenvalue(nr,ni);
715 : }
716 : #endif
717 4016 : PetscFunctionReturn(PETSC_SUCCESS);
718 : }
719 :
720 : /*@
721 : EPSComputeError - Computes the error (based on the residual norm) associated
722 : with the i-th computed eigenpair.
723 :
724 : Collective
725 :
726 : Input Parameters:
727 : + eps - the eigensolver context
728 : . i - the solution index
729 : - type - the type of error to compute
730 :
731 : Output Parameter:
732 : . error - the error
733 :
734 : Notes:
735 : The error can be computed in various ways, all of them based on the residual
736 : norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.
737 :
738 : Level: beginner
739 :
740 : .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
741 : @*/
742 3719 : PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
743 : {
744 3719 : Mat A,B;
745 3719 : Vec xr,xi,w[3];
746 3719 : PetscReal t,vecnorm=1.0,errorl;
747 3719 : PetscScalar kr,ki;
748 3719 : PetscBool flg;
749 :
750 3719 : PetscFunctionBegin;
751 3719 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
752 11157 : PetscValidLogicalCollectiveInt(eps,i,2);
753 11157 : PetscValidLogicalCollectiveEnum(eps,type,3);
754 3719 : PetscAssertPointer(error,4);
755 3719 : EPSCheckSolved(eps,1);
756 :
757 : /* allocate work vectors */
758 : #if defined(PETSC_USE_COMPLEX)
759 : PetscCall(EPSSetWorkVecs(eps,3));
760 : xi = NULL;
761 : w[1] = NULL;
762 : #else
763 3719 : PetscCall(EPSSetWorkVecs(eps,5));
764 3719 : xi = eps->work[3];
765 3719 : w[1] = eps->work[4];
766 : #endif
767 3719 : xr = eps->work[0];
768 3719 : w[0] = eps->work[1];
769 3719 : w[2] = eps->work[2];
770 :
771 : /* compute residual norm */
772 3719 : PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
773 3719 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));
774 :
775 : /* compute 2-norm of eigenvector */
776 3719 : if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));
777 :
778 : /* if two-sided, compute left residual norm and take the maximum */
779 3719 : if (eps->twosided) {
780 36 : PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
781 36 : PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
782 51 : *error = PetscMax(*error,errorl);
783 : }
784 :
785 : /* compute error */
786 3719 : switch (type) {
787 : case EPS_ERROR_ABSOLUTE:
788 : break;
789 2047 : case EPS_ERROR_RELATIVE:
790 2047 : *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
791 2047 : break;
792 1638 : case EPS_ERROR_BACKWARD:
793 : /* initialization of matrix norms */
794 1638 : if (!eps->nrma) {
795 34 : PetscCall(STGetMatrix(eps->st,0,&A));
796 34 : PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
797 34 : PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
798 34 : PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
799 : }
800 1638 : if (eps->isgeneralized) {
801 1632 : if (!eps->nrmb) {
802 33 : PetscCall(STGetMatrix(eps->st,1,&B));
803 33 : PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
804 33 : PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
805 33 : PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
806 : }
807 6 : } else eps->nrmb = 1.0;
808 1638 : t = SlepcAbsEigenvalue(kr,ki);
809 1638 : *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
810 1638 : break;
811 0 : default:
812 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
813 : }
814 3719 : PetscFunctionReturn(PETSC_SUCCESS);
815 : }
816 :
817 : /*
818 : EPSGetStartVector - Generate a suitable vector to be used as the starting vector
819 : for the recurrence that builds the right subspace.
820 :
821 : Collective
822 :
823 : Input Parameters:
824 : + eps - the eigensolver context
825 : - i - iteration number
826 :
827 : Output Parameters:
828 : . breakdown - flag indicating that a breakdown has occurred
829 :
830 : Notes:
831 : The start vector is computed from another vector: for the first step (i=0),
832 : the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
833 : vector is created. Then this vector is forced to be in the range of OP (only
834 : for generalized definite problems) and orthonormalized with respect to all
835 : V-vectors up to i-1. The resulting vector is placed in V[i].
836 :
837 : The flag breakdown is set to true if either i=0 and the vector belongs to the
838 : deflation space, or i>0 and the vector is linearly dependent with respect
839 : to the V-vectors.
840 : */
841 636 : PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
842 : {
843 636 : PetscReal norm;
844 636 : PetscBool lindep;
845 636 : Vec w,z;
846 :
847 636 : PetscFunctionBegin;
848 636 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
849 1908 : PetscValidLogicalCollectiveInt(eps,i,2);
850 :
851 : /* For the first step, use the first initial vector, otherwise a random one */
852 636 : if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));
853 :
854 : /* Force the vector to be in the range of OP for definite generalized problems */
855 636 : if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
856 112 : PetscCall(BVCreateVec(eps->V,&w));
857 112 : PetscCall(BVCopyVec(eps->V,i,w));
858 112 : PetscCall(BVGetColumn(eps->V,i,&z));
859 112 : PetscCall(STApply(eps->st,w,z));
860 112 : PetscCall(BVRestoreColumn(eps->V,i,&z));
861 112 : PetscCall(VecDestroy(&w));
862 : }
863 :
864 : /* Orthonormalize the vector with respect to previous vectors */
865 636 : PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
866 636 : if (breakdown) *breakdown = lindep;
867 599 : else if (lindep || norm == 0.0) {
868 0 : PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
869 0 : PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
870 : }
871 636 : PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
872 636 : PetscFunctionReturn(PETSC_SUCCESS);
873 : }
874 :
875 : /*
876 : EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
877 : vector for the recurrence that builds the left subspace. See EPSGetStartVector().
878 : */
879 24 : PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
880 : {
881 24 : PetscReal norm;
882 24 : PetscBool lindep;
883 :
884 24 : PetscFunctionBegin;
885 24 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
886 72 : PetscValidLogicalCollectiveInt(eps,i,2);
887 :
888 : /* For the first step, use the first initial vector, otherwise a random one */
889 24 : if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));
890 :
891 : /* Orthonormalize the vector with respect to previous vectors */
892 24 : PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
893 24 : if (breakdown) *breakdown = lindep;
894 18 : else if (lindep || norm == 0.0) {
895 0 : PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
896 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
897 : }
898 24 : PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
899 24 : PetscFunctionReturn(PETSC_SUCCESS);
900 : }
|
