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>

  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, ldz = 1;

 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;

 21:   PetscCall(PetscArrayzero(c, nmax));
 22:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
 23:   d  = cgP->d;
 24:   e  = cgP->e;
 25:   ee = cgP->ee;

 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:   }

 33:   PetscCall(PetscBLASIntCast(n, &bn));
 34:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 35:   PetscCallLAPACKInfo("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, ee, NULL, &ldz, NULL, &info));
 36:   PetscCall(PetscFPTrapPop());
 37:   PetscCall(PetscSortReal(n, r));
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 41: PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP ksp, PetscReal *emax, PetscReal *emin)
 42: {
 43:   KSP_CG      *cgP = (KSP_CG *)ksp->data;
 44:   PetscScalar *d, *e;
 45:   PetscReal   *dd, *ee;
 46:   PetscInt     n = ksp->its;
 47:   PetscBLASInt bn, ldz = 1;

 49:   PetscFunctionBegin;
 50:   if (!n) {
 51:     *emax = *emin = 1.0;
 52:     PetscFunctionReturn(PETSC_SUCCESS);
 53:   }
 54:   d  = cgP->d;
 55:   e  = cgP->e;
 56:   dd = cgP->dd;
 57:   ee = cgP->ee;

 59:   /* copy tridiagonal matrix to work space */
 60:   for (PetscInt j = 0; j < n; j++) {
 61:     dd[j] = PetscRealPart(d[j]);
 62:     ee[j] = PetscRealPart(e[j + 1]);
 63:   }

 65:   PetscCall(PetscBLASIntCast(n, &bn));
 66:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 67:   PetscCallLAPACKInfo("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, ee, NULL, &ldz, NULL, &info));
 68:   PetscCall(PetscFPTrapPop());
 69:   *emin = dd[0];
 70:   *emax = dd[n - 1];
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }