Actual source code: pepsolve.c
slepc-3.18.2 2023-01-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: PEP routines related to the solution process
13: References:
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: */
20: #include <slepc/private/pepimpl.h>
21: #include <slepc/private/bvimpl.h>
22: #include <petscdraw.h>
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";
37: PetscErrorCode PEPComputeVectors(PEP pep)
38: {
39: PEPCheckSolved(pep,1);
40: if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,computevectors);
41: pep->state = PEP_STATE_EIGENVECTORS;
42: return 0;
43: }
45: PetscErrorCode PEPExtractVectors(PEP pep)
46: {
47: PEPCheckSolved(pep,1);
48: if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,extractvectors);
49: return 0;
50: }
52: /*@
53: PEPSolve - Solves the polynomial eigensystem.
55: Collective on pep
57: Input Parameter:
58: . pep - eigensolver context obtained from PEPCreate()
60: Options Database Keys:
61: + -pep_view - print information about the solver used
62: . -pep_view_matk - view the coefficient matrix Ak (replace k by an integer from 0 to nmat-1)
63: . -pep_view_vectors - view the computed eigenvectors
64: . -pep_view_values - view the computed eigenvalues
65: . -pep_converged_reason - print reason for convergence, and number of iterations
66: . -pep_error_absolute - print absolute errors of each eigenpair
67: . -pep_error_relative - print relative errors of each eigenpair
68: - -pep_error_backward - print backward errors of each eigenpair
70: Notes:
71: All the command-line options listed above admit an optional argument specifying
72: the viewer type and options. For instance, use '-pep_view_mat0 binary:amatrix.bin'
73: to save the A matrix to a binary file, '-pep_view_values draw' to draw the computed
74: eigenvalues graphically, or '-pep_error_relative :myerr.m:ascii_matlab' to save
75: the errors in a file that can be executed in Matlab.
77: Level: beginner
79: .seealso: PEPCreate(), PEPSetUp(), PEPDestroy(), PEPSetTolerances()
80: @*/
81: PetscErrorCode PEPSolve(PEP pep)
82: {
83: PetscInt i,k;
84: PetscBool flg,islinear;
85: char str[16];
88: if (pep->state>=PEP_STATE_SOLVED) return 0;
89: PetscLogEventBegin(PEP_Solve,pep,0,0,0);
91: /* call setup */
92: PEPSetUp(pep);
93: pep->nconv = 0;
94: pep->its = 0;
95: k = pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1);
96: for (i=0;i<k;i++) {
97: pep->eigr[i] = 0.0;
98: pep->eigi[i] = 0.0;
99: pep->errest[i] = 0.0;
100: pep->perm[i] = i;
101: }
102: PEPViewFromOptions(pep,NULL,"-pep_view_pre");
103: RGViewFromOptions(pep->rg,NULL,"-rg_view");
105: /* Call solver */
106: PetscUseTypeMethod(pep,solve);
108: pep->state = PEP_STATE_SOLVED;
110: /* Only the first nconv columns contain useful information */
111: BVSetActiveColumns(pep->V,0,pep->nconv);
113: PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear);
114: if (!islinear) {
115: STPostSolve(pep->st);
116: /* Map eigenvalues back to the original problem */
117: STGetTransform(pep->st,&flg);
118: if (flg) PetscTryTypeMethod(pep,backtransform);
119: }
121: #if !defined(PETSC_USE_COMPLEX)
122: /* reorder conjugate eigenvalues (positive imaginary first) */
123: for (i=0;i<pep->nconv-1;i++) {
124: if (pep->eigi[i] != 0) {
125: if (pep->eigi[i] < 0) {
126: pep->eigi[i] = -pep->eigi[i];
127: pep->eigi[i+1] = -pep->eigi[i+1];
128: /* the next correction only works with eigenvectors */
129: PEPComputeVectors(pep);
130: BVScaleColumn(pep->V,i+1,-1.0);
131: }
132: i++;
133: }
134: }
135: #endif
137: if (pep->refine!=PEP_REFINE_NONE) PetscCitationsRegister(citation,&cited);
139: if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
140: PEPComputeVectors(pep);
141: PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv);
142: }
144: /* sort eigenvalues according to pep->which parameter */
145: SlepcSortEigenvalues(pep->sc,pep->nconv,pep->eigr,pep->eigi,pep->perm);
146: PetscLogEventEnd(PEP_Solve,pep,0,0,0);
148: /* various viewers */
149: PEPViewFromOptions(pep,NULL,"-pep_view");
150: PEPConvergedReasonViewFromOptions(pep);
151: PEPErrorViewFromOptions(pep);
152: PEPValuesViewFromOptions(pep);
153: PEPVectorsViewFromOptions(pep);
154: for (i=0;i<pep->nmat;i++) {
155: PetscSNPrintf(str,sizeof(str),"-pep_view_mat%" PetscInt_FMT,i);
156: MatViewFromOptions(pep->A[i],(PetscObject)pep,str);
157: }
159: /* Remove the initial subspace */
160: pep->nini = 0;
161: return 0;
162: }
164: /*@
165: PEPGetIterationNumber - Gets the current iteration number. If the
166: call to PEPSolve() is complete, then it returns the number of iterations
167: carried out by the solution method.
169: Not Collective
171: Input Parameter:
172: . pep - the polynomial eigensolver context
174: Output Parameter:
175: . its - number of iterations
177: Note:
178: During the i-th iteration this call returns i-1. If PEPSolve() is
179: complete, then parameter "its" contains either the iteration number at
180: which convergence was successfully reached, or failure was detected.
181: Call PEPGetConvergedReason() to determine if the solver converged or
182: failed and why.
184: Level: intermediate
186: .seealso: PEPGetConvergedReason(), PEPSetTolerances()
187: @*/
188: PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
189: {
192: *its = pep->its;
193: return 0;
194: }
196: /*@
197: PEPGetConverged - Gets the number of converged eigenpairs.
199: Not Collective
201: Input Parameter:
202: . pep - the polynomial eigensolver context
204: Output Parameter:
205: . nconv - number of converged eigenpairs
207: Note:
208: This function should be called after PEPSolve() has finished.
210: Level: beginner
212: .seealso: PEPSetDimensions(), PEPSolve(), PEPGetEigenpair()
213: @*/
214: PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
215: {
218: PEPCheckSolved(pep,1);
219: *nconv = pep->nconv;
220: return 0;
221: }
223: /*@
224: PEPGetConvergedReason - Gets the reason why the PEPSolve() iteration was
225: stopped.
227: Not Collective
229: Input Parameter:
230: . pep - the polynomial eigensolver context
232: Output Parameter:
233: . reason - negative value indicates diverged, positive value converged
235: Options Database Key:
236: . -pep_converged_reason - print the reason to a viewer
238: Notes:
239: Possible values for reason are
240: + PEP_CONVERGED_TOL - converged up to tolerance
241: . PEP_CONVERGED_USER - converged due to a user-defined condition
242: . PEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
243: . PEP_DIVERGED_BREAKDOWN - generic breakdown in method
244: - PEP_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry
246: Can only be called after the call to PEPSolve() is complete.
248: Level: intermediate
250: .seealso: PEPSetTolerances(), PEPSolve(), PEPConvergedReason
251: @*/
252: PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
253: {
256: PEPCheckSolved(pep,1);
257: *reason = pep->reason;
258: return 0;
259: }
261: /*@C
262: PEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
263: PEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
265: Logically Collective on pep
267: Input Parameters:
268: + pep - polynomial eigensolver context
269: - i - index of the solution
271: Output Parameters:
272: + eigr - real part of eigenvalue
273: . eigi - imaginary part of eigenvalue
274: . Vr - real part of eigenvector
275: - Vi - imaginary part of eigenvector
277: Notes:
278: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
279: required. Otherwise, the caller must provide valid Vec objects, i.e.,
280: they must be created by the calling program with e.g. MatCreateVecs().
282: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
283: configured with complex scalars the eigenvalue is stored
284: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
285: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
286: them is not required.
288: The index i should be a value between 0 and nconv-1 (see PEPGetConverged()).
289: Eigenpairs are indexed according to the ordering criterion established
290: with PEPSetWhichEigenpairs().
292: Level: beginner
294: .seealso: PEPSolve(), PEPGetConverged(), PEPSetWhichEigenpairs()
295: @*/
296: PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
297: {
298: PetscInt k;
304: PEPCheckSolved(pep,1);
308: PEPComputeVectors(pep);
309: k = pep->perm[i];
311: /* eigenvalue */
312: #if defined(PETSC_USE_COMPLEX)
313: if (eigr) *eigr = pep->eigr[k];
314: if (eigi) *eigi = 0;
315: #else
316: if (eigr) *eigr = pep->eigr[k];
317: if (eigi) *eigi = pep->eigi[k];
318: #endif
320: /* eigenvector */
321: BV_GetEigenvector(pep->V,k,pep->eigi[k],Vr,Vi);
322: return 0;
323: }
325: /*@
326: PEPGetErrorEstimate - Returns the error estimate associated to the i-th
327: computed eigenpair.
329: Not Collective
331: Input Parameters:
332: + pep - polynomial eigensolver context
333: - i - index of eigenpair
335: Output Parameter:
336: . errest - the error estimate
338: Notes:
339: This is the error estimate used internally by the eigensolver. The actual
340: error bound can be computed with PEPComputeError(). See also the users
341: manual for details.
343: Level: advanced
345: .seealso: PEPComputeError()
346: @*/
347: PetscErrorCode PEPGetErrorEstimate(PEP pep,PetscInt i,PetscReal *errest)
348: {
351: PEPCheckSolved(pep,1);
354: *errest = pep->errest[pep->perm[i]];
355: return 0;
356: }
358: /*
359: PEPComputeResidualNorm_Private - Computes the norm of the residual vector
360: associated with an eigenpair.
362: Input Parameters:
363: kr,ki - eigenvalue
364: xr,xi - eigenvector
365: z - array of 4 work vectors (z[2],z[3] not referenced in complex scalars)
366: */
367: PetscErrorCode PEPComputeResidualNorm_Private(PEP pep,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
368: {
369: Mat *A=pep->A;
370: PetscInt i,nmat=pep->nmat;
371: PetscScalar t[20],*vals=t,*ivals=NULL;
372: Vec u,w;
373: #if !defined(PETSC_USE_COMPLEX)
374: Vec ui,wi;
375: PetscReal ni;
376: PetscBool imag;
377: PetscScalar it[20];
378: #endif
380: u = z[0]; w = z[1];
381: VecSet(u,0.0);
382: #if !defined(PETSC_USE_COMPLEX)
383: ui = z[2]; wi = z[3];
384: ivals = it;
385: #endif
386: if (nmat>20) {
387: PetscMalloc1(nmat,&vals);
388: #if !defined(PETSC_USE_COMPLEX)
389: PetscMalloc1(nmat,&ivals);
390: #endif
391: }
392: PEPEvaluateBasis(pep,kr,ki,vals,ivals);
393: #if !defined(PETSC_USE_COMPLEX)
394: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
395: imag = PETSC_FALSE;
396: else {
397: imag = PETSC_TRUE;
398: VecSet(ui,0.0);
399: }
400: #endif
401: for (i=0;i<nmat;i++) {
402: if (vals[i]!=0.0) {
403: MatMult(A[i],xr,w);
404: VecAXPY(u,vals[i],w);
405: }
406: #if !defined(PETSC_USE_COMPLEX)
407: if (imag) {
408: if (ivals[i]!=0 || vals[i]!=0) {
409: MatMult(A[i],xi,wi);
410: if (vals[i]==0) MatMult(A[i],xr,w);
411: }
412: if (ivals[i]!=0) {
413: VecAXPY(u,-ivals[i],wi);
414: VecAXPY(ui,ivals[i],w);
415: }
416: if (vals[i]!=0) VecAXPY(ui,vals[i],wi);
417: }
418: #endif
419: }
420: VecNorm(u,NORM_2,norm);
421: #if !defined(PETSC_USE_COMPLEX)
422: if (imag) {
423: VecNorm(ui,NORM_2,&ni);
424: *norm = SlepcAbsEigenvalue(*norm,ni);
425: }
426: #endif
427: if (nmat>20) {
428: PetscFree(vals);
429: #if !defined(PETSC_USE_COMPLEX)
430: PetscFree(ivals);
431: #endif
432: }
433: return 0;
434: }
436: /*@
437: PEPComputeError - Computes the error (based on the residual norm) associated
438: with the i-th computed eigenpair.
440: Collective on pep
442: Input Parameters:
443: + pep - the polynomial eigensolver context
444: . i - the solution index
445: - type - the type of error to compute
447: Output Parameter:
448: . error - the error
450: Notes:
451: The error can be computed in various ways, all of them based on the residual
452: norm ||P(l)x||_2 where l is the eigenvalue and x is the eigenvector.
453: See the users guide for additional details.
455: Level: beginner
457: .seealso: PEPErrorType, PEPSolve(), PEPGetErrorEstimate()
458: @*/
459: PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
460: {
461: Vec xr,xi,w[4];
462: PetscScalar kr,ki;
463: PetscReal t,z=0.0;
464: PetscInt j;
465: PetscBool flg;
471: PEPCheckSolved(pep,1);
473: /* allocate work vectors */
474: #if defined(PETSC_USE_COMPLEX)
475: PEPSetWorkVecs(pep,3);
476: xi = NULL;
477: w[2] = NULL;
478: w[3] = NULL;
479: #else
480: PEPSetWorkVecs(pep,6);
481: xi = pep->work[3];
482: w[2] = pep->work[4];
483: w[3] = pep->work[5];
484: #endif
485: xr = pep->work[0];
486: w[0] = pep->work[1];
487: w[1] = pep->work[2];
489: /* compute residual norms */
490: PEPGetEigenpair(pep,i,&kr,&ki,xr,xi);
491: PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error);
493: /* compute error */
494: switch (type) {
495: case PEP_ERROR_ABSOLUTE:
496: break;
497: case PEP_ERROR_RELATIVE:
498: *error /= SlepcAbsEigenvalue(kr,ki);
499: break;
500: case PEP_ERROR_BACKWARD:
501: /* initialization of matrix norms */
502: if (!pep->nrma[pep->nmat-1]) {
503: for (j=0;j<pep->nmat;j++) {
504: MatHasOperation(pep->A[j],MATOP_NORM,&flg);
506: MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]);
507: }
508: }
509: t = SlepcAbsEigenvalue(kr,ki);
510: for (j=pep->nmat-1;j>=0;j--) {
511: z = z*t+pep->nrma[j];
512: }
513: *error /= z;
514: break;
515: default:
516: SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
517: }
518: return 0;
519: }