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

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)
    z.split()[1].assign(lm)
    z.split()[2].assign(eigenmode)

ksp.getPC().setType("lu")
                    if value == "default":
                        value = wrapping.get_default_linear_solver()
                    if hasattr(ksp.getPC(), "setFactorSolverType"):  # PETSc >= 3.9
                        ksp.getPC().setFactorSolverType(value)
                    else:
                        ksp.getPC().setFactorSolverPackage(value)
                elif key == "maximum_iterations":
                    eps_tolerances[1] = value
                elif key == "problem_type":
                    if value == "hermitian":
                        self.eps.setProblemType(SLEPc.EPS.ProblemType.HEP)
                    elif value == "non_hermitian":
                        self.eps.setProblemType(SLEPc.EPS.ProblemType.NHEP)
                    elif value == "gen_hermitian":
                        self.eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
                    elif value == "gen_non_hermitian":
                        self.eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)
                    elif value == "pos_gen_non_hermitian":
                        self.eps.setProblemType(SLEPc.EPS.ProblemType.PGNHEP)
                    else:
                        raise RuntimeError("Invalid problem type")
                elif key == "solver":
                    if value == "power":
                        self.eps.setType(SLEPc.EPS.Type.POWER)
                    elif value == "subspace":
                        self.eps.setType(SLEPc.EPS.Type.SUBSPACE)
                    elif value == "arnoldi":
                        self.eps.setType(SLEPc.EPS.Type.ARNOLDI)
                    elif value == "lanczos":
                        self.eps.setType(SLEPc.EPS.Type.LANCZOS)
                    elif value == "krylov-schur":

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)
        assert nconv > 0
        eigenvalues = np.r_[eigenvalues, 
                            np.array([E.getEigenvalue(i).real
                                      for i in range(nconv)])]
    return size, eigenvalues