Actual source code: pepsolve.c
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 eigenproblem.
57: Collective
59: Input Parameter:
60: . pep - the polynomial eigensolver context
62: Options Database Keys:
63: + -pep_view - print information about the solver once the solve is complete
64: . -pep_view_pre - print information about the solver before the solve starts
65: . -pep_view_matk - view the coefficient matrix $A_k$ (replace `k` by an integer from 0 to `nmat`-1)
66: . -pep_view_vectors - view the computed eigenvectors
67: . -pep_view_values - view the computed eigenvalues
68: . -pep_converged_reason - print reason for convergence/divergence, and number of iterations
69: . -pep_error_absolute - print absolute errors of each eigenpair
70: . -pep_error_relative - print relative errors of each eigenpair
71: - -pep_error_backward - print backward errors of each eigenpair
73: Notes:
74: The problem matrices are specified with `PEPSetOperators()`.
76: `PEPSolve()` will return without generating an error regardless of whether
77: all requested solutions were computed or not. Call `PEPGetConverged()` to get the
78: actual number of computed solutions, and `PEPGetConvergedReason()` to determine if
79: the solver converged or failed and why.
81: All the command-line options listed above admit an optional argument specifying
82: the viewer type and options. For instance, use `-pep_view_mat0 binary:matrix0.bin`
83: to save the $A_0$ matrix to a binary file, `-pep_view_values draw` to draw the computed
84: eigenvalues graphically, or `-pep_error_relative :myerr.m:ascii_matlab` to save
85: the errors in a file that can be executed in Matlab.
86: See `PetscObjectViewFromOptions()` for more details.
88: Level: beginner
90: .seealso: [](ch:pep), `PEPCreate()`, `PEPSetUp()`, `PEPDestroy()`, `PEPSetTolerances()`, `PEPGetConverged()`, `PEPGetConvergedReason()`
91: @*/
92: PetscErrorCode PEPSolve(PEP pep)
93: {
94: PetscInt i,k;
95: PetscBool flg,islinear;
96: char str[16];
98: PetscFunctionBegin;
100: if (pep->state>=PEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
101: PetscCall(PetscLogEventBegin(PEP_Solve,pep,0,0,0));
103: /* call setup */
104: PetscCall(PEPSetUp(pep));
105: pep->nconv = 0;
106: pep->its = 0;
107: k = pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1);
108: for (i=0;i<k;i++) {
109: pep->eigr[i] = 0.0;
110: pep->eigi[i] = 0.0;
111: pep->errest[i] = 0.0;
112: pep->perm[i] = i;
113: }
114: PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view_pre"));
115: PetscCall(RGViewFromOptions(pep->rg,NULL,"-rg_view"));
117: /* Call solver */
118: PetscUseTypeMethod(pep,solve);
119: PetscCheck(pep->reason,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
120: pep->state = PEP_STATE_SOLVED;
122: /* Only the first nconv columns contain useful information */
123: PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
125: PetscCall(PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear));
126: if (!islinear) {
127: PetscCall(STPostSolve(pep->st));
128: /* Map eigenvalues back to the original problem */
129: PetscCall(STGetTransform(pep->st,&flg));
130: if (flg) PetscTryTypeMethod(pep,backtransform);
131: }
133: #if !defined(PETSC_USE_COMPLEX)
134: /* reorder conjugate eigenvalues (positive imaginary first) */
135: for (i=0;i<pep->nconv-1;i++) {
136: if (pep->eigi[i] != 0) {
137: if (pep->eigi[i] < 0) {
138: pep->eigi[i] = -pep->eigi[i];
139: pep->eigi[i+1] = -pep->eigi[i+1];
140: /* the next correction only works with eigenvectors */
141: PetscCall(PEPComputeVectors(pep));
142: PetscCall(BVScaleColumn(pep->V,i+1,-1.0));
143: }
144: i++;
145: }
146: }
147: #endif
149: if (pep->refine!=PEP_REFINE_NONE) PetscCall(PetscCitationsRegister(citation,&cited));
151: if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
152: PetscCall(PEPComputeVectors(pep));
153: PetscCall(PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv));
154: }
156: /* sort eigenvalues according to pep->which parameter */
157: PetscCall(SlepcSortEigenvalues(pep->sc,pep->nconv,pep->eigr,pep->eigi,pep->perm));
158: PetscCall(PetscLogEventEnd(PEP_Solve,pep,0,0,0));
160: /* various viewers */
161: PetscCall(PEPViewFromOptions(pep,NULL,"-pep_view"));
162: PetscCall(PEPConvergedReasonViewFromOptions(pep));
163: PetscCall(PEPErrorViewFromOptions(pep));
164: PetscCall(PEPValuesViewFromOptions(pep));
165: PetscCall(PEPVectorsViewFromOptions(pep));
166: for (i=0;i<pep->nmat;i++) {
167: PetscCall(PetscSNPrintf(str,sizeof(str),"-pep_view_mat%" PetscInt_FMT,i));
168: PetscCall(MatViewFromOptions(pep->A[i],(PetscObject)pep,str));
169: }
171: /* Remove the initial subspace */
172: pep->nini = 0;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*@
177: PEPGetIterationNumber - Gets the current iteration number. If the
178: call to `PEPSolve()` is complete, then it returns the number of iterations
179: carried out by the solution method.
181: Not Collective
183: Input Parameter:
184: . pep - the polynomial eigensolver context
186: Output Parameter:
187: . its - number of iterations
189: Note:
190: During the $i$-th iteration this call returns $i-1$. If `PEPSolve()` is
191: complete, then parameter `its` contains either the iteration number at
192: which convergence was successfully reached, or failure was detected.
193: Call `PEPGetConvergedReason()` to determine if the solver converged or
194: failed and why.
196: Level: intermediate
198: .seealso: [](ch:pep), `PEPGetConvergedReason()`, `PEPSetTolerances()`
199: @*/
200: PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
201: {
202: PetscFunctionBegin;
204: PetscAssertPointer(its,2);
205: *its = pep->its;
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*@
210: PEPGetConverged - Gets the number of converged eigenpairs.
212: Not Collective
214: Input Parameter:
215: . pep - the polynomial eigensolver context
217: Output Parameter:
218: . nconv - number of converged eigenpairs
220: Notes:
221: This function should be called after `PEPSolve()` has finished.
223: The value `nconv` may be different from the number of requested solutions
224: `nev`, but not larger than `ncv`, see `PEPSetDimensions()`.
226: Level: beginner
228: .seealso: [](ch:pep), `PEPSetDimensions()`, `PEPSolve()`, `PEPGetEigenpair()`
229: @*/
230: PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
231: {
232: PetscFunctionBegin;
234: PetscAssertPointer(nconv,2);
235: PEPCheckSolved(pep,1);
236: *nconv = pep->nconv;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@
241: PEPGetConvergedReason - Gets the reason why the `PEPSolve()` iteration was
242: stopped.
244: Not Collective
246: Input Parameter:
247: . pep - the polynomial eigensolver context
249: Output Parameter:
250: . reason - negative value indicates diverged, positive value converged, see
251: `PEPConvergedReason` for the possible values
253: Options Database Key:
254: . -pep_converged_reason - print reason for convergence/divergence, and number of iterations
256: Note:
257: If this routine is called before or doing the `PEPSolve()` the value of
258: `PEP_CONVERGED_ITERATING` is returned.
260: Level: intermediate
262: .seealso: [](ch:pep), `PEPSetTolerances()`, `PEPSolve()`, `PEPConvergedReason`
263: @*/
264: PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
265: {
266: PetscFunctionBegin;
268: PetscAssertPointer(reason,2);
269: PEPCheckSolved(pep,1);
270: *reason = pep->reason;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: PEPGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
276: `PEPSolve()`. The solution consists in both the eigenvalue and the eigenvector.
278: Collective
280: Input Parameters:
281: + pep - the polynomial eigensolver context
282: - i - index of the solution
284: Output Parameters:
285: + eigr - real part of eigenvalue
286: . eigi - imaginary part of eigenvalue
287: . Vr - real part of eigenvector
288: - Vi - imaginary part of eigenvector
290: Notes:
291: It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
292: required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
293: they must be created by the calling program with e.g. `MatCreateVecs()`.
295: If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
296: configured with complex scalars the eigenvalue is stored
297: directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi` is
298: set to zero). In any case, the user can pass `NULL` in `Vr` or `Vi` if one of
299: them is not required.
301: The index `i` should be a value between 0 and `nconv`-1 (see `PEPGetConverged()`).
302: Eigenpairs are indexed according to the ordering criterion established
303: with `PEPSetWhichEigenpairs()`.
305: The eigenvector is normalized to have unit norm.
307: Level: beginner
309: .seealso: [](ch:pep), `PEPSolve()`, `PEPGetConverged()`, `PEPSetWhichEigenpairs()`
310: @*/
311: PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
312: {
313: PetscInt k;
315: PetscFunctionBegin;
320: PEPCheckSolved(pep,1);
321: PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
322: PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
324: PetscCall(PEPComputeVectors(pep));
325: k = pep->perm[i];
327: /* eigenvalue */
328: #if defined(PETSC_USE_COMPLEX)
329: if (eigr) *eigr = pep->eigr[k];
330: if (eigi) *eigi = 0;
331: #else
332: if (eigr) *eigr = pep->eigr[k];
333: if (eigi) *eigi = pep->eigi[k];
334: #endif
336: /* eigenvector */
337: PetscCall(BV_GetEigenvector(pep->V,k,pep->eigi[k],Vr,Vi));
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: PEPGetErrorEstimate - Returns the error estimate associated to the `i`-th
343: computed eigenpair.
345: Not Collective
347: Input Parameters:
348: + pep - the polynomial eigensolver context
349: - i - index of eigenpair
351: Output Parameter:
352: . errest - the error estimate
354: Note:
355: This is the error estimate used internally by the eigensolver. The actual
356: error bound can be computed with `PEPComputeError()`.
358: Level: advanced
360: .seealso: [](ch:pep), `PEPComputeError()`
361: @*/
362: PetscErrorCode PEPGetErrorEstimate(PEP pep,PetscInt i,PetscReal *errest)
363: {
364: PetscFunctionBegin;
366: PetscAssertPointer(errest,3);
367: PEPCheckSolved(pep,1);
368: PetscCheck(i>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
369: PetscCheck(i<pep->nconv,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");
370: *errest = pep->errest[pep->perm[i]];
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*
375: PEPComputeResidualNorm_Private - Computes the norm of the residual vector
376: associated with an eigenpair.
378: Input Parameters:
379: kr,ki - eigenvalue
380: xr,xi - eigenvector
381: z - array of 4 work vectors (z[2],z[3] not referenced in complex scalars)
382: */
383: PetscErrorCode PEPComputeResidualNorm_Private(PEP pep,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
384: {
385: Mat *A=pep->A;
386: PetscInt i,nmat=pep->nmat;
387: PetscScalar t[20],*vals=t,*ivals=NULL;
388: Vec u,w;
389: #if !defined(PETSC_USE_COMPLEX)
390: Vec ui,wi;
391: PetscReal ni;
392: PetscBool imag;
393: PetscScalar it[20];
394: #endif
396: PetscFunctionBegin;
397: u = z[0]; w = z[1];
398: PetscCall(VecSet(u,0.0));
399: #if !defined(PETSC_USE_COMPLEX)
400: ui = z[2]; wi = z[3];
401: ivals = it;
402: #endif
403: if (nmat>20) {
404: PetscCall(PetscMalloc1(nmat,&vals));
405: #if !defined(PETSC_USE_COMPLEX)
406: PetscCall(PetscMalloc1(nmat,&ivals));
407: #endif
408: }
409: PetscCall(PEPEvaluateBasis(pep,kr,ki,vals,ivals));
410: #if !defined(PETSC_USE_COMPLEX)
411: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
412: imag = PETSC_FALSE;
413: else {
414: imag = PETSC_TRUE;
415: PetscCall(VecSet(ui,0.0));
416: }
417: #endif
418: for (i=0;i<nmat;i++) {
419: if (vals[i]!=0.0) {
420: PetscCall(MatMult(A[i],xr,w));
421: PetscCall(VecAXPY(u,vals[i],w));
422: }
423: #if !defined(PETSC_USE_COMPLEX)
424: if (imag) {
425: if (ivals[i]!=0 || vals[i]!=0) {
426: PetscCall(MatMult(A[i],xi,wi));
427: if (vals[i]==0) PetscCall(MatMult(A[i],xr,w));
428: }
429: if (ivals[i]!=0) {
430: PetscCall(VecAXPY(u,-ivals[i],wi));
431: PetscCall(VecAXPY(ui,ivals[i],w));
432: }
433: if (vals[i]!=0) PetscCall(VecAXPY(ui,vals[i],wi));
434: }
435: #endif
436: }
437: PetscCall(VecNorm(u,NORM_2,norm));
438: #if !defined(PETSC_USE_COMPLEX)
439: if (imag) {
440: PetscCall(VecNorm(ui,NORM_2,&ni));
441: *norm = SlepcAbsEigenvalue(*norm,ni);
442: }
443: #endif
444: if (nmat>20) {
445: PetscCall(PetscFree(vals));
446: #if !defined(PETSC_USE_COMPLEX)
447: PetscCall(PetscFree(ivals));
448: #endif
449: }
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@
454: PEPComputeError - Computes the error (based on the residual norm) associated
455: with the `i`-th computed eigenpair.
457: Collective
459: Input Parameters:
460: + pep - the polynomial eigensolver context
461: . i - the solution index
462: - type - the type of error to compute, see `PEPErrorType`
464: Output Parameter:
465: . error - the error
467: Note:
468: The error can be computed in various ways, all of them based on the residual
469: norm $\|P(\lambda)x\|_2$ where $(\lambda,x)$ is the approximate eigenpair.
471: Level: beginner
473: .seealso: [](ch:pep), `PEPErrorType`, `PEPSolve()`, `PEPGetErrorEstimate()`
474: @*/
475: PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
476: {
477: Vec xr,xi,w[4];
478: PetscScalar kr,ki;
479: PetscReal t,z=0.0;
480: PetscInt j;
481: PetscBool flg;
483: PetscFunctionBegin;
487: PetscAssertPointer(error,4);
488: PEPCheckSolved(pep,1);
490: /* allocate work vectors */
491: #if defined(PETSC_USE_COMPLEX)
492: PetscCall(PEPSetWorkVecs(pep,3));
493: xi = NULL;
494: w[2] = NULL;
495: w[3] = NULL;
496: #else
497: PetscCall(PEPSetWorkVecs(pep,6));
498: xi = pep->work[3];
499: w[2] = pep->work[4];
500: w[3] = pep->work[5];
501: #endif
502: xr = pep->work[0];
503: w[0] = pep->work[1];
504: w[1] = pep->work[2];
506: /* compute residual norms */
507: PetscCall(PEPGetEigenpair(pep,i,&kr,&ki,xr,xi));
508: PetscCall(PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error));
510: /* compute error */
511: switch (type) {
512: case PEP_ERROR_ABSOLUTE:
513: break;
514: case PEP_ERROR_RELATIVE:
515: *error /= SlepcAbsEigenvalue(kr,ki);
516: break;
517: case PEP_ERROR_BACKWARD:
518: /* initialization of matrix norms */
519: if (!pep->nrma[pep->nmat-1]) {
520: for (j=0;j<pep->nmat;j++) {
521: PetscCall(MatHasOperation(pep->A[j],MATOP_NORM,&flg));
522: PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
523: PetscCall(MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]));
524: }
525: }
526: t = SlepcAbsEigenvalue(kr,ki);
527: for (j=pep->nmat-1;j>=0;j--) {
528: z = z*t+pep->nrma[j];
529: }
530: *error /= z;
531: break;
532: default:
533: SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
534: }
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }