1 /*
2 Code for calculating extreme eigenvalues via the Lanczos method
3 running with CG. Note this only works for symmetric real and Hermitian
4 matrices (not complex matrices that are symmetric).
5 */
6 #include <../src/ksp/ksp/impls/cg/cgimpl.h>
7 #include <../include/petscblaslapack.h>
8
KSPComputeEigenvalues_CG(KSP ksp,PetscInt nmax,PetscReal * r,PetscReal * c,PetscInt * neig)9 PetscErrorCode KSPComputeEigenvalues_CG(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
10 {
11 KSP_CG *cgP = (KSP_CG *)ksp->data;
12 PetscScalar *d, *e;
13 PetscReal *ee;
14 PetscInt n = ksp->its;
15 PetscBLASInt bn, lierr = 0, ldz = 1;
16
17 PetscFunctionBegin;
18 PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
19 *neig = n;
20
21 PetscCall(PetscArrayzero(c, nmax));
22 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
23 d = cgP->d;
24 e = cgP->e;
25 ee = cgP->ee;
26
27 /* copy tridiagonal matrix to work space */
28 for (PetscInt j = 0; j < n; j++) {
29 r[j] = PetscRealPart(d[j]);
30 ee[j] = PetscRealPart(e[j + 1]);
31 }
32
33 PetscCall(PetscBLASIntCast(n, &bn));
34 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
35 PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, ee, NULL, &ldz, NULL, &lierr));
36 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
37 PetscCall(PetscFPTrapPop());
38 PetscCall(PetscSortReal(n, r));
39 PetscFunctionReturn(PETSC_SUCCESS);
40 }
41
KSPComputeExtremeSingularValues_CG(KSP ksp,PetscReal * emax,PetscReal * emin)42 PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP ksp, PetscReal *emax, PetscReal *emin)
43 {
44 KSP_CG *cgP = (KSP_CG *)ksp->data;
45 PetscScalar *d, *e;
46 PetscReal *dd, *ee;
47 PetscInt n = ksp->its;
48 PetscBLASInt bn, lierr = 0, ldz = 1;
49
50 PetscFunctionBegin;
51 if (!n) {
52 *emax = *emin = 1.0;
53 PetscFunctionReturn(PETSC_SUCCESS);
54 }
55 d = cgP->d;
56 e = cgP->e;
57 dd = cgP->dd;
58 ee = cgP->ee;
59
60 /* copy tridiagonal matrix to work space */
61 for (PetscInt j = 0; j < n; j++) {
62 dd[j] = PetscRealPart(d[j]);
63 ee[j] = PetscRealPart(e[j + 1]);
64 }
65
66 PetscCall(PetscBLASIntCast(n, &bn));
67 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
68 PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, ee, NULL, &ldz, NULL, &lierr));
69 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
70 PetscCall(PetscFPTrapPop());
71 *emin = dd[0];
72 *emax = dd[n - 1];
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75