How to use the slepc4py.SLEPc.EPS function in slepc4py

# Now solve eigenvalue problem for $F_u(u, \lambda)\phi = r\phi$
    # Want eigenmode phi with minimal eigenvalue r
    B = derivative(residual(th, lm, TestFunction(V)), th, TrialFunction(V))

    petsc_M = assemble(inner(TestFunction(V), TrialFunction(V))*dx, bcs=bcs).petscmat
    petsc_B = assemble(B, bcs=bcs).petscmat

    num_eigenvalues = 1

    opts = PETSc.Options()
    opts.setValue("eps_target_magnitude", None)
    opts.setValue("eps_target", 0)
    opts.setValue("st_type", "sinvert")

    es = SLEPc.EPS().create(comm=COMM_WORLD)
    es.setDimensions(num_eigenvalues)
    es.setOperators(petsc_B, petsc_M)
    es.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    es.setFromOptions()
    es.solve()

    ev_re, ev_im = petsc_B.getVecs()
    es.getEigenpair(0, ev_re, ev_im)
    eigenmode = Function(V)
    eigenmode.vector().set_local(ev_re)

    Z = MixedFunctionSpace([V, R, V])

    # Set initial guesses for state, parameter, null eigenmode
    z = Function(Z)
    z.split()[0].assign(th)

def eig_slepc(AA, BB, **kwargs):
    '''Slepc solver'''
    assert has_slepc4py()
    from slepc4py import SLEPc

    nev = 3  # Number of eigenvalues
    eigenvalues = np.array([])
    size = AA.size[0]
    info('System size: %i' % size)
    eps_type = kwargs.get('eps_type', SLEPc.EPS.Type.GD)
    for which in (SLEPc.EPS.Which.SMALLEST_MAGNITUDE, SLEPc.EPS.Which.LARGEST_MAGNITUDE):
        # Setup the eigensolver
        E = SLEPc.EPS().create()
        E.setOperators(AA ,BB)
        E.setType(E.Type.GD)
        E.setDimensions(nev, PETSc.DECIDE)
        E.setTolerances(1E-6, 8000)

        E.setWhichEigenpairs(which)
        E.setProblemType(SLEPc.EPS.ProblemType.GHEP)
        E.setFromOptions()

        # Solve the eigensystem
        E.solve()

        its = E.getIterationNumber()
        info('Number of iterations of the method: %i' % its)
        nconv = E.getConverged()
        info('Number of converged eigenpairs: %d' % nconv)

# 6 fields
        Nfields = 6
        self._A = PETSc.Mat()
        self._A.create(PETSc.COMM_WORLD)
        self._A.setSizes([Nfields*self._M*self._N, Nfields*self._M*self._N])
        self._A.setType('aij')
        self._A.setUp()

        self._B = PETSc.Mat()
        self._B.create(PETSc.COMM_WORLD)
        self._B.setSizes([Nfields*self._M*self._N, Nfields*self._M*self._N])
        self._B.setType('aij')
        self._B.setUp()

        # setup the solver
        self._solver = SLEPc.EPS()
        self._solver.create()

        # we need to set up the spectral transformation so that it doesnt try
        # to invert B
        st = self._solver.getST()
        st.setType('sinvert')

        # Let's use MUMPS for any system solving since it is fast
        ksp = st.getKSP()
        #ksp.setType('gmres')
        ksp.setType('preonly')
        pc = ksp.getPC()
        pc.setType('lu')

        # backwards compatibility
        try: