Actual source code: pepsolve.c
slepc-3.22.1 2024-10-28
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: PetscFunctionBegin;
40: PEPCheckSolved(pep,1);
41: if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,computevectors);
42: pep->state = PEP_STATE_EIGENVECTORS;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: PetscErrorCode PEPExtractVectors(PEP pep)
47: {
48: PetscFunctionBegin;
49: PEPCheckSolved(pep,1);
50: if (pep->state==PEP_STATE_SOLVED) PetscTryTypeMethod(pep,extractvectors);
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: /*@
55: PEPSolve - Solves the polynomial eigensystem.
57: Collective
59: Input Parameter:
60: . pep - eigensolver context obtained from PEPCreate()
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
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.
79: Level: beginner
81: .seealso: PEPCreate(), PEPSetUp(), PEPDestroy(), PEPSetTolerances()
82: @*/
83: PetscErrorCode PEPSolve(PEP pep)
84: {
85: PetscInt i,k;
86: PetscBool flg,islinear;
87: char str[16];
89: PetscFunctionBegin;
91: if (pep->state>=PEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
92: PetscCall(PetscLogEventBegin(PEP_Solve,pep,0,0,0));
94: /* call setup */
95: PetscCall(PEPSetUp(pep));
96: pep->nconv = 0;
97: pep->its = 0;
98: k = pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1);
99: for (i=0;i<k;i++) {
100: pep->eigr[i] = 0.0;
101: pep->eigi[i] = 0.0;
102: pep->errest[i] = 0.0;
103: pep->perm[i] = i;
104: }
105: PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view_pre"));
106: PetscCall(RGViewFromOptions(pep->rg,NULL,"-rg_view"));
108: /* Call solver */
109: PetscUseTypeMethod(pep,solve);
110: PetscCheck(pep->reason,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
111: pep->state = PEP_STATE_SOLVED;
113: /* Only the first nconv columns contain useful information */
114: PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
116: PetscCall(PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear));
117: if (!islinear) {
118: PetscCall(STPostSolve(pep->st));
119: /* Map eigenvalues back to the original problem */
120: PetscCall(STGetTransform(pep->st,&flg));
121: if (flg) PetscTryTypeMethod(pep,backtransform);
122: }
124: #if !defined(PETSC_USE_COMPLEX)
125: /* reorder conjugate eigenvalues (positive imaginary first) */
126: for (i=0;i<pep->nconv-1;i++) {
127: if (pep->eigi[i] != 0) {
128: if (pep->eigi[i] < 0) {
129: pep->eigi[i] = -pep->eigi[i];
130: pep->eigi[i+1] = -pep->eigi[i+1];
131: /* the next correction only works with eigenvectors */
132: PetscCall(PEPComputeVectors(pep));
133: PetscCall(BVScaleColumn(pep->V,i+1,-1.0));
134: }
135: i++;
136: }
137: }
138: #endif
140: if (pep->refine!=PEP_REFINE_NONE) PetscCall(PetscCitationsRegister(citation,&cited));
142: if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
143: PetscCall(PEPComputeVectors(pep));
144: PetscCall(PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv));
145: }
147: /* sort eigenvalues according to pep->which parameter */
148: PetscCall(SlepcSortEigenvalues(pep->sc,pep->nconv,pep->eigr,pep->eigi,pep->perm));
149: PetscCall(PetscLogEventEnd(PEP_Solve,pep,0,0,0));
151: /* various viewers */
152: PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view"));
153: PetscCall(PEPConvergedReasonViewFromOptions(pep));
154: PetscCall(PEPErrorViewFromOptions(pep));
155: PetscCall(PEPValuesViewFromOptions(pep));
156: PetscCall(PEPVectorsViewFromOptions(pep));
157: for (i=0;i<pep->nmat;i++) {
158: PetscCall(PetscSNPrintf(str,sizeof(str),"-pep_view_mat%" PetscInt_FMT,i));
159: PetscCall(MatViewFromOptions(pep->A[i],(PetscObject)pep,str));
160: }
162: /* Remove the initial subspace */
163: pep->nini = 0;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
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.
172: Not Collective
174: Input Parameter:
175: . pep - the polynomial eigensolver context
177: Output Parameter:
178: . its - number of iterations
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.
187: Level: intermediate
189: .seealso: PEPGetConvergedReason(), PEPSetTolerances()
190: @*/
191: PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
192: {
193: PetscFunctionBegin;
195: PetscAssertPointer(its,2);
196: *its = pep->its;
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /*@
201: PEPGetConverged - Gets the number of converged eigenpairs.
203: Not Collective
205: Input Parameter:
206: . pep - the polynomial eigensolver context
208: Output Parameter:
209: . nconv - number of converged eigenpairs
211: Note:
212: This function should be called after PEPSolve() has finished.
214: Level: beginner
216: .seealso: PEPSetDimensions(), PEPSolve(), PEPGetEigenpair()
217: @*/
218: PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
219: {
220: PetscFunctionBegin;
222: PetscAssertPointer(nconv,2);
223: PEPCheckSolved(pep,1);
224: *nconv = pep->nconv;
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /*@
229: PEPGetConvergedReason - Gets the reason why the PEPSolve() iteration was
230: stopped.
232: Not Collective
234: Input Parameter:
235: . pep - the polynomial eigensolver context
237: Output Parameter:
238: . reason - negative value indicates diverged, positive value converged
240: Options Database Key:
241: . -pep_converged_reason - print the reason to a viewer
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
251: Can only be called after the call to PEPSolve() is complete.
253: Level: intermediate
255: .seealso: PEPSetTolerances(), PEPSolve(), PEPConvergedReason
256: @*/
257: PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
258: {
259: PetscFunctionBegin;
261: PetscAssertPointer(reason,2);
262: PEPCheckSolved(pep,1);
263: *reason = pep->reason;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
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.
271: Collective
273: Input Parameters:
274: + pep - polynomial eigensolver context
275: - i - index of the solution
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
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().
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.
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().
298: Level: beginner
300: .seealso: PEPSolve(), PEPGetConverged(), PEPSetWhichEigenpairs()
301: @*/
302: PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
303: {
304: PetscInt k;
306: PetscFunctionBegin;
311: PEPCheckSolved(pep,1);
312: PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
313: PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
315: PetscCall(PEPComputeVectors(pep));
316: k = pep->perm[i];
318: /* eigenvalue */
319: #if defined(PETSC_USE_COMPLEX)
320: if (eigr) *eigr = pep->eigr[k];
321: if (eigi) *eigi = 0;
322: #else
323: if (eigr) *eigr = pep->eigr[k];
324: if (eigi) *eigi = pep->eigi[k];
325: #endif
327: /* eigenvector */
328: PetscCall(BV_GetEigenvector(pep->V,k,pep->eigi[k],Vr,Vi));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /*@
333: PEPGetErrorEstimate - Returns the error estimate associated to the i-th
334: computed eigenpair.
336: Not Collective
338: Input Parameters:
339: + pep - polynomial eigensolver context
340: - i - index of eigenpair
342: Output Parameter:
343: . errest - the error estimate
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.
350: Level: advanced
352: .seealso: PEPComputeError()
353: @*/
354: PetscErrorCode PEPGetErrorEstimate(PEP pep,PetscInt i,PetscReal *errest)
355: {
356: PetscFunctionBegin;
358: PetscAssertPointer(errest,3);
359: PEPCheckSolved(pep,1);
360: PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
361: PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
362: *errest = pep->errest[pep->perm[i]];
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*
367: PEPComputeResidualNorm_Private - Computes the norm of the residual vector
368: associated with an eigenpair.
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: PetscErrorCode PEPComputeResidualNorm_Private(PEP pep,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
376: {
377: Mat *A=pep->A;
378: PetscInt i,nmat=pep->nmat;
379: PetscScalar t[20],*vals=t,*ivals=NULL;
380: Vec u,w;
381: #if !defined(PETSC_USE_COMPLEX)
382: Vec ui,wi;
383: PetscReal ni;
384: PetscBool imag;
385: PetscScalar it[20];
386: #endif
388: PetscFunctionBegin;
389: u = z[0]; w = z[1];
390: PetscCall(VecSet(u,0.0));
391: #if !defined(PETSC_USE_COMPLEX)
392: ui = z[2]; wi = z[3];
393: ivals = it;
394: #endif
395: if (nmat>20) {
396: PetscCall(PetscMalloc1(nmat,&vals));
397: #if !defined(PETSC_USE_COMPLEX)
398: PetscCall(PetscMalloc1(nmat,&ivals));
399: #endif
400: }
401: PetscCall(PEPEvaluateBasis(pep,kr,ki,vals,ivals));
402: #if !defined(PETSC_USE_COMPLEX)
403: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
404: imag = PETSC_FALSE;
405: else {
406: imag = PETSC_TRUE;
407: PetscCall(VecSet(ui,0.0));
408: }
409: #endif
410: for (i=0;i<nmat;i++) {
411: if (vals[i]!=0.0) {
412: PetscCall(MatMult(A[i],xr,w));
413: PetscCall(VecAXPY(u,vals[i],w));
414: }
415: #if !defined(PETSC_USE_COMPLEX)
416: if (imag) {
417: if (ivals[i]!=0 || vals[i]!=0) {
418: PetscCall(MatMult(A[i],xi,wi));
419: if (vals[i]==0) PetscCall(MatMult(A[i],xr,w));
420: }
421: if (ivals[i]!=0) {
422: PetscCall(VecAXPY(u,-ivals[i],wi));
423: PetscCall(VecAXPY(ui,ivals[i],w));
424: }
425: if (vals[i]!=0) PetscCall(VecAXPY(ui,vals[i],wi));
426: }
427: #endif
428: }
429: PetscCall(VecNorm(u,NORM_2,norm));
430: #if !defined(PETSC_USE_COMPLEX)
431: if (imag) {
432: PetscCall(VecNorm(ui,NORM_2,&ni));
433: *norm = SlepcAbsEigenvalue(*norm,ni);
434: }
435: #endif
436: if (nmat>20) {
437: PetscCall(PetscFree(vals));
438: #if !defined(PETSC_USE_COMPLEX)
439: PetscCall(PetscFree(ivals));
440: #endif
441: }
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: /*@
446: PEPComputeError - Computes the error (based on the residual norm) associated
447: with the i-th computed eigenpair.
449: Collective
451: Input Parameters:
452: + pep - the polynomial eigensolver context
453: . i - the solution index
454: - type - the type of error to compute
456: Output Parameter:
457: . error - the error
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.
464: Level: beginner
466: .seealso: PEPErrorType, PEPSolve(), PEPGetErrorEstimate()
467: @*/
468: PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
469: {
470: Vec xr,xi,w[4];
471: PetscScalar kr,ki;
472: PetscReal t,z=0.0;
473: PetscInt j;
474: PetscBool flg;
476: PetscFunctionBegin;
480: PetscAssertPointer(error,4);
481: PEPCheckSolved(pep,1);
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: PetscCall(PEPSetWorkVecs(pep,6));
491: xi = pep->work[3];
492: w[2] = pep->work[4];
493: w[3] = pep->work[5];
494: #endif
495: xr = pep->work[0];
496: w[0] = pep->work[1];
497: w[1] = pep->work[2];
499: /* compute residual norms */
500: PetscCall(PEPGetEigenpair(pep,i,&kr,&ki,xr,xi));
501: PetscCall(PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error));
503: /* compute error */
504: switch (type) {
505: case PEP_ERROR_ABSOLUTE:
506: break;
507: case PEP_ERROR_RELATIVE:
508: *error /= SlepcAbsEigenvalue(kr,ki);
509: break;
510: case PEP_ERROR_BACKWARD:
511: /* initialization of matrix norms */
512: if (!pep->nrma[pep->nmat-1]) {
513: for (j=0;j<pep->nmat;j++) {
514: PetscCall(MatHasOperation(pep->A[j],MATOP_NORM,&flg));
515: PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
516: PetscCall(MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]));
517: }
518: }
519: t = SlepcAbsEigenvalue(kr,ki);
520: for (j=pep->nmat-1;j>=0;j--) {
521: z = z*t+pep->nrma[j];
522: }
523: *error /= z;
524: break;
525: default:
526: SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
527: }
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }