1: /*
  2:    Code for calculating extreme eigenvalues via the Lanczo 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, lierr = 0, 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]);
 31:   }

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

 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;

 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;

 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]);
 64:   }

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