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 : PEP routines related to the solution process
12 :
13 : References:
14 :
15 : [1] C. Campos and J.E. Roman, "Parallel iterative refinement in
16 : polynomial eigenvalue problems", Numer. Linear Algebra Appl.
17 : 23(4):730-745, 2016.
18 : */
19 :
20 : #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
21 : #include <slepc/private/bvimpl.h>
22 : #include <petscdraw.h>
23 :
24 : static PetscBool cited = PETSC_FALSE;
25 : static const char citation[] =
26 : "@Article{slepc-pep-refine,\n"
27 : " author = \"C. Campos and J. E. Roman\",\n"
28 : " title = \"Parallel iterative refinement in polynomial eigenvalue problems\",\n"
29 : " journal = \"Numer. Linear Algebra Appl.\",\n"
30 : " volume = \"23\",\n"
31 : " number = \"4\",\n"
32 : " pages = \"730--745\",\n"
33 : " year = \"2016,\"\n"
34 : " doi = \"https://doi.org/10.1002/nla.2052\"\n"
35 : "}\n";
36 :
37 1554 : PetscErrorCode PEPComputeVectors(PEP pep)
38 : {
39 1554 : PetscFunctionBegin;
40 1554 : PEPCheckSolved(pep,1);
41 1554 : if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,computevectors);
42 1554 : pep->state = PEP_STATE_EIGENVECTORS;
43 1554 : PetscFunctionReturn(PETSC_SUCCESS);
44 : }
45 :
46 160 : PetscErrorCode PEPExtractVectors(PEP pep)
47 : {
48 160 : PetscFunctionBegin;
49 160 : PEPCheckSolved(pep,1);
50 160 : if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,extractvectors);
51 160 : PetscFunctionReturn(PETSC_SUCCESS);
52 : }
53 :
54 : /*@
55 : PEPSolve - Solves the polynomial eigensystem.
56 :
57 : Collective
58 :
59 : Input Parameter:
60 : . pep - eigensolver context obtained from PEPCreate()
61 :
62 : Options Database Keys:
63 : + -pep_view - print information about the solver used
64 : . -pep_view_matk - view the coefficient matrix Ak (replace k by an integer from 0 to nmat-1)
65 : . -pep_view_vectors - view the computed eigenvectors
66 : . -pep_view_values - view the computed eigenvalues
67 : . -pep_converged_reason - print reason for convergence, and number of iterations
68 : . -pep_error_absolute - print absolute errors of each eigenpair
69 : . -pep_error_relative - print relative errors of each eigenpair
70 : - -pep_error_backward - print backward errors of each eigenpair
71 :
72 : Notes:
73 : All the command-line options listed above admit an optional argument specifying
74 : the viewer type and options. For instance, use '-pep_view_mat0 binary:amatrix.bin'
75 : to save the A matrix to a binary file, '-pep_view_values draw' to draw the computed
76 : eigenvalues graphically, or '-pep_error_relative :myerr.m:ascii_matlab' to save
77 : the errors in a file that can be executed in Matlab.
78 :
79 : Level: beginner
80 :
81 : .seealso: PEPCreate(), PEPSetUp(), PEPDestroy(), PEPSetTolerances()
82 : @*/
83 173 : PetscErrorCode PEPSolve(PEP pep)
84 : {
85 173 : PetscInt i,k;
86 173 : PetscBool flg,islinear;
87 173 : char str[16];
88 :
89 173 : PetscFunctionBegin;
90 173 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
91 173 : if (pep->state>=PEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
92 173 : PetscCall(PetscLogEventBegin(PEP_Solve,pep,0,0,0));
93 :
94 : /* call setup */
95 173 : PetscCall(PEPSetUp(pep));
96 173 : pep->nconv = 0;
97 173 : pep->its = 0;
98 173 : k = pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1);
99 3939 : for (i=0;i<k;i++) {
100 3766 : pep->eigr[i] = 0.0;
101 3766 : pep->eigi[i] = 0.0;
102 3766 : pep->errest[i] = 0.0;
103 3766 : pep->perm[i] = i;
104 : }
105 173 : PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view_pre"));
106 173 : PetscCall(RGViewFromOptions(pep->rg,NULL,"-rg_view"));
107 :
108 : /* Call solver */
109 173 : PetscUseTypeMethod(pep,solve);
110 173 : PetscCheck(pep->reason,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
111 173 : pep->state = PEP_STATE_SOLVED;
112 :
113 : /* Only the first nconv columns contain useful information */
114 173 : PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
115 :
116 173 : PetscCall(PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear));
117 173 : if (!islinear) {
118 139 : PetscCall(STPostSolve(pep->st));
119 : /* Map eigenvalues back to the original problem */
120 139 : PetscCall(STGetTransform(pep->st,&flg));
121 139 : if (flg) PetscTryTypeMethod(pep,backtransform);
122 : }
123 :
124 : #if !defined(PETSC_USE_COMPLEX)
125 : /* reorder conjugate eigenvalues (positive imaginary first) */
126 1035 : for (i=0;i<pep->nconv-1;i++) {
127 862 : if (pep->eigi[i] != 0) {
128 116 : if (pep->eigi[i] < 0) {
129 35 : pep->eigi[i] = -pep->eigi[i];
130 35 : pep->eigi[i+1] = -pep->eigi[i+1];
131 : /* the next correction only works with eigenvectors */
132 35 : PetscCall(PEPComputeVectors(pep));
133 35 : PetscCall(BVScaleColumn(pep->V,i+1,-1.0));
134 : }
135 116 : i++;
136 : }
137 : }
138 : #endif
139 :
140 173 : if (pep->refine!=PEP_REFINE_NONE) PetscCall(PetscCitationsRegister(citation,&cited));
141 :
142 173 : if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
143 15 : PetscCall(PEPComputeVectors(pep));
144 15 : PetscCall(PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv));
145 : }
146 :
147 : /* sort eigenvalues according to pep->which parameter */
148 173 : PetscCall(SlepcSortEigenvalues(pep->sc,pep->nconv,pep->eigr,pep->eigi,pep->perm));
149 173 : PetscCall(PetscLogEventEnd(PEP_Solve,pep,0,0,0));
150 :
151 : /* various viewers */
152 173 : PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view"));
153 173 : PetscCall(PEPConvergedReasonViewFromOptions(pep));
154 173 : PetscCall(PEPErrorViewFromOptions(pep));
155 173 : PetscCall(PEPValuesViewFromOptions(pep));
156 173 : PetscCall(PEPVectorsViewFromOptions(pep));
157 831 : for (i=0;i<pep->nmat;i++) {
158 658 : PetscCall(PetscSNPrintf(str,sizeof(str),"-pep_view_mat%" PetscInt_FMT,i));
159 658 : PetscCall(MatViewFromOptions(pep->A[i],(PetscObject)pep,str));
160 : }
161 :
162 : /* Remove the initial subspace */
163 173 : pep->nini = 0;
164 173 : PetscFunctionReturn(PETSC_SUCCESS);
165 : }
166 :
167 : /*@
168 : PEPGetIterationNumber - Gets the current iteration number. If the
169 : call to PEPSolve() is complete, then it returns the number of iterations
170 : carried out by the solution method.
171 :
172 : Not Collective
173 :
174 : Input Parameter:
175 : . pep - the polynomial eigensolver context
176 :
177 : Output Parameter:
178 : . its - number of iterations
179 :
180 : Note:
181 : During the i-th iteration this call returns i-1. If PEPSolve() is
182 : complete, then parameter "its" contains either the iteration number at
183 : which convergence was successfully reached, or failure was detected.
184 : Call PEPGetConvergedReason() to determine if the solver converged or
185 : failed and why.
186 :
187 : Level: intermediate
188 :
189 : .seealso: PEPGetConvergedReason(), PEPSetTolerances()
190 : @*/
191 30 : PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
192 : {
193 30 : PetscFunctionBegin;
194 30 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
195 30 : PetscAssertPointer(its,2);
196 30 : *its = pep->its;
197 30 : PetscFunctionReturn(PETSC_SUCCESS);
198 : }
199 :
200 : /*@
201 : PEPGetConverged - Gets the number of converged eigenpairs.
202 :
203 : Not Collective
204 :
205 : Input Parameter:
206 : . pep - the polynomial eigensolver context
207 :
208 : Output Parameter:
209 : . nconv - number of converged eigenpairs
210 :
211 : Note:
212 : This function should be called after PEPSolve() has finished.
213 :
214 : Level: beginner
215 :
216 : .seealso: PEPSetDimensions(), PEPSolve(), PEPGetEigenpair()
217 : @*/
218 85 : PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
219 : {
220 85 : PetscFunctionBegin;
221 85 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
222 85 : PetscAssertPointer(nconv,2);
223 85 : PEPCheckSolved(pep,1);
224 85 : *nconv = pep->nconv;
225 85 : PetscFunctionReturn(PETSC_SUCCESS);
226 : }
227 :
228 : /*@
229 : PEPGetConvergedReason - Gets the reason why the PEPSolve() iteration was
230 : stopped.
231 :
232 : Not Collective
233 :
234 : Input Parameter:
235 : . pep - the polynomial eigensolver context
236 :
237 : Output Parameter:
238 : . reason - negative value indicates diverged, positive value converged
239 :
240 : Options Database Key:
241 : . -pep_converged_reason - print the reason to a viewer
242 :
243 : Notes:
244 : Possible values for reason are
245 : + PEP_CONVERGED_TOL - converged up to tolerance
246 : . PEP_CONVERGED_USER - converged due to a user-defined condition
247 : . PEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
248 : . PEP_DIVERGED_BREAKDOWN - generic breakdown in method
249 : - PEP_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
250 :
251 : Can only be called after the call to PEPSolve() is complete.
252 :
253 : Level: intermediate
254 :
255 : .seealso: PEPSetTolerances(), PEPSolve(), PEPConvergedReason
256 : @*/
257 24 : PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
258 : {
259 24 : PetscFunctionBegin;
260 24 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
261 24 : PetscAssertPointer(reason,2);
262 24 : PEPCheckSolved(pep,1);
263 24 : *reason = pep->reason;
264 24 : PetscFunctionReturn(PETSC_SUCCESS);
265 : }
266 :
267 : /*@
268 : PEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
269 : PEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
270 :
271 : Collective
272 :
273 : Input Parameters:
274 : + pep - polynomial eigensolver context
275 : - i - index of the solution
276 :
277 : Output Parameters:
278 : + eigr - real part of eigenvalue
279 : . eigi - imaginary part of eigenvalue
280 : . Vr - real part of eigenvector
281 : - Vi - imaginary part of eigenvector
282 :
283 : Notes:
284 : It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
285 : required. Otherwise, the caller must provide valid Vec objects, i.e.,
286 : they must be created by the calling program with e.g. MatCreateVecs().
287 :
288 : If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
289 : configured with complex scalars the eigenvalue is stored
290 : directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
291 : set to zero). In any case, the user can pass NULL in Vr or Vi if one of
292 : them is not required.
293 :
294 : The index i should be a value between 0 and nconv-1 (see PEPGetConverged()).
295 : Eigenpairs are indexed according to the ordering criterion established
296 : with PEPSetWhichEigenpairs().
297 :
298 : Level: beginner
299 :
300 : .seealso: PEPSolve(), PEPGetConverged(), PEPSetWhichEigenpairs()
301 : @*/
302 1503 : PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
303 : {
304 1503 : PetscInt k;
305 :
306 1503 : PetscFunctionBegin;
307 1503 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
308 4509 : PetscValidLogicalCollectiveInt(pep,i,2);
309 1503 : if (Vr) { PetscValidHeaderSpecific(Vr,VEC_CLASSID,5); PetscCheckSameComm(pep,1,Vr,5); }
310 1503 : if (Vi) { PetscValidHeaderSpecific(Vi,VEC_CLASSID,6); PetscCheckSameComm(pep,1,Vi,6); }
311 1503 : PEPCheckSolved(pep,1);
312 1503 : PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
313 1503 : PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
314 :
315 1503 : PetscCall(PEPComputeVectors(pep));
316 1503 : k = pep->perm[i];
317 :
318 : /* eigenvalue */
319 : #if defined(PETSC_USE_COMPLEX)
320 : if (eigr) *eigr = pep->eigr[k];
321 : if (eigi) *eigi = 0;
322 : #else
323 1503 : if (eigr) *eigr = pep->eigr[k];
324 1503 : if (eigi) *eigi = pep->eigi[k];
325 : #endif
326 :
327 : /* eigenvector */
328 1503 : PetscCall(BV_GetEigenvector(pep->V,k,pep->eigi[k],Vr,Vi));
329 1503 : PetscFunctionReturn(PETSC_SUCCESS);
330 : }
331 :
332 : /*@
333 : PEPGetErrorEstimate - Returns the error estimate associated to the i-th
334 : computed eigenpair.
335 :
336 : Not Collective
337 :
338 : Input Parameters:
339 : + pep - polynomial eigensolver context
340 : - i - index of eigenpair
341 :
342 : Output Parameter:
343 : . errest - the error estimate
344 :
345 : Notes:
346 : This is the error estimate used internally by the eigensolver. The actual
347 : error bound can be computed with PEPComputeError(). See also the users
348 : manual for details.
349 :
350 : Level: advanced
351 :
352 : .seealso: PEPComputeError()
353 : @*/
354 4 : PetscErrorCode PEPGetErrorEstimate(PEP pep,PetscInt i,PetscReal *errest)
355 : {
356 4 : PetscFunctionBegin;
357 4 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
358 4 : PetscAssertPointer(errest,3);
359 4 : PEPCheckSolved(pep,1);
360 4 : PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
361 4 : PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
362 4 : *errest = pep->errest[pep->perm[i]];
363 4 : PetscFunctionReturn(PETSC_SUCCESS);
364 : }
365 :
366 : /*
367 : PEPComputeResidualNorm_Private - Computes the norm of the residual vector
368 : associated with an eigenpair.
369 :
370 : Input Parameters:
371 : kr,ki - eigenvalue
372 : xr,xi - eigenvector
373 : z - array of 4 work vectors (z[2],z[3] not referenced in complex scalars)
374 : */
375 1027 : PetscErrorCode PEPComputeResidualNorm_Private(PEP pep,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
376 : {
377 1027 : Mat *A=pep->A;
378 1027 : PetscInt i,nmat=pep->nmat;
379 1027 : PetscScalar t[20],*vals=t,*ivals=NULL;
380 1027 : Vec u,w;
381 : #if !defined(PETSC_USE_COMPLEX)
382 1027 : Vec ui,wi;
383 1027 : PetscReal ni;
384 1027 : PetscBool imag;
385 1027 : PetscScalar it[20];
386 : #endif
387 :
388 1027 : PetscFunctionBegin;
389 1027 : u = z[0]; w = z[1];
390 1027 : PetscCall(VecSet(u,0.0));
391 : #if !defined(PETSC_USE_COMPLEX)
392 1027 : ui = z[2]; wi = z[3];
393 1027 : ivals = it;
394 : #endif
395 1027 : if (nmat>20) {
396 0 : PetscCall(PetscMalloc1(nmat,&vals));
397 : #if !defined(PETSC_USE_COMPLEX)
398 0 : PetscCall(PetscMalloc1(nmat,&ivals));
399 : #endif
400 : }
401 1027 : PetscCall(PEPEvaluateBasis(pep,kr,ki,vals,ivals));
402 : #if !defined(PETSC_USE_COMPLEX)
403 1027 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
404 : imag = PETSC_FALSE;
405 : else {
406 201 : imag = PETSC_TRUE;
407 201 : PetscCall(VecSet(ui,0.0));
408 : }
409 : #endif
410 4173 : for (i=0;i<nmat;i++) {
411 3146 : if (vals[i]!=0.0) {
412 3146 : PetscCall(MatMult(A[i],xr,w));
413 3146 : PetscCall(VecAXPY(u,vals[i],w));
414 : }
415 : #if !defined(PETSC_USE_COMPLEX)
416 3146 : if (imag) {
417 659 : if (ivals[i]!=0 || vals[i]!=0) {
418 659 : PetscCall(MatMult(A[i],xi,wi));
419 659 : if (vals[i]==0) PetscCall(MatMult(A[i],xr,w));
420 : }
421 659 : if (ivals[i]!=0) {
422 458 : PetscCall(VecAXPY(u,-ivals[i],wi));
423 458 : PetscCall(VecAXPY(ui,ivals[i],w));
424 : }
425 3146 : if (vals[i]!=0) PetscCall(VecAXPY(ui,vals[i],wi));
426 : }
427 : #endif
428 : }
429 1027 : PetscCall(VecNorm(u,NORM_2,norm));
430 : #if !defined(PETSC_USE_COMPLEX)
431 1027 : if (imag) {
432 201 : PetscCall(VecNorm(ui,NORM_2,&ni));
433 201 : *norm = SlepcAbsEigenvalue(*norm,ni);
434 : }
435 : #endif
436 1027 : if (nmat>20) {
437 0 : PetscCall(PetscFree(vals));
438 : #if !defined(PETSC_USE_COMPLEX)
439 0 : PetscCall(PetscFree(ivals));
440 : #endif
441 : }
442 1027 : PetscFunctionReturn(PETSC_SUCCESS);
443 : }
444 :
445 : /*@
446 : PEPComputeError - Computes the error (based on the residual norm) associated
447 : with the i-th computed eigenpair.
448 :
449 : Collective
450 :
451 : Input Parameters:
452 : + pep - the polynomial eigensolver context
453 : . i - the solution index
454 : - type - the type of error to compute
455 :
456 : Output Parameter:
457 : . error - the error
458 :
459 : Notes:
460 : The error can be computed in various ways, all of them based on the residual
461 : norm ||P(l)x||_2 where l is the eigenvalue and x is the eigenvector.
462 : See the users guide for additional details.
463 :
464 : Level: beginner
465 :
466 : .seealso: PEPErrorType, PEPSolve(), PEPGetErrorEstimate()
467 : @*/
468 952 : PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
469 : {
470 952 : Vec xr,xi,w[4];
471 952 : PetscScalar kr,ki;
472 952 : PetscReal t,z=0.0;
473 952 : PetscInt j;
474 952 : PetscBool flg;
475 :
476 952 : PetscFunctionBegin;
477 952 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
478 2856 : PetscValidLogicalCollectiveInt(pep,i,2);
479 2856 : PetscValidLogicalCollectiveEnum(pep,type,3);
480 952 : PetscAssertPointer(error,4);
481 952 : PEPCheckSolved(pep,1);
482 :
483 : /* allocate work vectors */
484 : #if defined(PETSC_USE_COMPLEX)
485 : PetscCall(PEPSetWorkVecs(pep,3));
486 : xi = NULL;
487 : w[2] = NULL;
488 : w[3] = NULL;
489 : #else
490 952 : PetscCall(PEPSetWorkVecs(pep,6));
491 952 : xi = pep->work[3];
492 952 : w[2] = pep->work[4];
493 952 : w[3] = pep->work[5];
494 : #endif
495 952 : xr = pep->work[0];
496 952 : w[0] = pep->work[1];
497 952 : w[1] = pep->work[2];
498 :
499 : /* compute residual norms */
500 952 : PetscCall(PEPGetEigenpair(pep,i,&kr,&ki,xr,xi));
501 952 : PetscCall(PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error));
502 :
503 : /* compute error */
504 952 : switch (type) {
505 : case PEP_ERROR_ABSOLUTE:
506 : break;
507 34 : case PEP_ERROR_RELATIVE:
508 34 : *error /= SlepcAbsEigenvalue(kr,ki);
509 34 : break;
510 906 : case PEP_ERROR_BACKWARD:
511 : /* initialization of matrix norms */
512 906 : if (!pep->nrma[pep->nmat-1]) {
513 554 : for (j=0;j<pep->nmat;j++) {
514 419 : PetscCall(MatHasOperation(pep->A[j],MATOP_NORM,&flg));
515 419 : PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
516 419 : PetscCall(MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]));
517 : }
518 : }
519 906 : t = SlepcAbsEigenvalue(kr,ki);
520 3680 : for (j=pep->nmat-1;j>=0;j--) {
521 2774 : z = z*t+pep->nrma[j];
522 : }
523 906 : *error /= z;
524 906 : break;
525 0 : default:
526 0 : SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
527 : }
528 952 : PetscFunctionReturn(PETSC_SUCCESS);
529 : }
|
