diff --git a/REQUIRE b/REQUIRE
index 35f3490..445fa43 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -3,3 +3,4 @@ LightGraphs 1.2
 JuMP 0.18
 MathProgBase 0.7
 BlossomV 0.4
+Hungarian 0.2
diff --git a/src/LightGraphsMatching.jl b/src/LightGraphsMatching.jl
index aebf53d..e909ccd 100644
--- a/src/LightGraphsMatching.jl
+++ b/src/LightGraphsMatching.jl
@@ -7,8 +7,9 @@ using SparseArrays: spzeros
 using JuMP
 using MathProgBase: AbstractMathProgSolver
 import BlossomV # 'using BlossomV'  leads to naming conflicts with JuMP
+using Hungarian
 
-export MatchingResult, maximum_weight_matching, maximum_weight_maximal_matching, minimum_weight_perfect_matching
+export MatchingResult, maximum_weight_matching, maximum_weight_maximal_matching, minimum_weight_perfect_matching, HungarianAlgorithm, LPAlgorithm
 
 """
     struct MatchingResult{U}
@@ -31,6 +32,7 @@ end
 include("lp.jl")
 include("maximum_weight_matching.jl")
 include("blossomv.jl")
+include("hungarian.jl")
+include("maximum_weight_maximal_matching.jl")
 
 end # module
-
diff --git a/src/hungarian.jl b/src/hungarian.jl
new file mode 100644
index 0000000..5349fe3
--- /dev/null
+++ b/src/hungarian.jl
@@ -0,0 +1,65 @@
+function maximum_weight_maximal_matching_hungarian(g::Graph,
+          w::AbstractMatrix{T}=default_weights(g)) where {T <: Real}
+  edge_list = collect(edges(g))
+  n = nv(g)
+
+  # Determine the bipartition of the graph. 
+  bipartition = bipartite_map(g)
+  if length(bipartition) != n # Equivalent to !is_bipartite(g), but reuses the results of the previous function call.
+    error("The Hungarian algorithm only works for bipartite graphs; otherwise, prefer the Blossom algorithm (not yet available in LightGraphsMatching")
+  end
+  n_first = count(bipartition .== 1)
+  n_second = count(bipartition .== 2)
+
+  to_bipartition_1 = [count(bipartition[1:i] .== 1) for i in 1:n]
+  to_bipartition_2 = [count(bipartition[1:i] .== 2) for i in 1:n]
+
+  # hungarian() minimises the total cost, while this function is supposed to maximise the total weights.
+  wDual = maximum(w) .- w
+
+  # Remove weights that are not in the graph (Hungarian.jl considers all weights that are not missing values as real edges).
+  # Assume w is symmetric, so that the weight of matching i->j is the same as the one for j->i. 
+  weights = Matrix{Union{Missing, T}}(missing, n_first, n_second)
+
+  for i in 1:n
+    for j in 1:n
+      if Edge(i, j) ∈ edge_list || Edge(j, i) ∈ edge_list 
+        if bipartition[i] == 1 # and bipartition[j] == 2
+          idx_first = to_bipartition_1[i]
+          idx_second = to_bipartition_2[j]
+        else # bipartition[i] == 2 and bipartition[j] == 1
+          idx_first = to_bipartition_1[j]
+          idx_second = to_bipartition_2[i]
+        end
+
+        weight_to_add = (Edge(i, j) ∈ edge_list) ? wDual[i, j] : wDual[j, i]
+
+        weights[idx_first, idx_second] = weight_to_add
+      end
+    end
+  end
+
+  # Run the Hungarian algorithm.
+  assignment, _ = hungarian(weights)
+
+  # Convert the output format to match LGMatching's. 
+  pairs = Tuple{Int, Int}[]
+  mate = fill(-1, n) # Initialise to unmatched. 
+  for i in eachindex(assignment)
+    if assignment[i] != 0 # If matched: 
+      original_i = findfirst(to_bipartition_1 .== i)
+      original_j = findfirst(to_bipartition_2 .== assignment[i])
+
+      mate[original_i] = original_j
+      mate[original_j] = original_i
+
+      push!(pairs, (original_i, original_j))
+    end
+  end
+
+  # Compute the cost for this matching (as weights had to be changed for Hungarian.jl, the one returned by hungarian() makes no sense). 
+  cost = sum(w[p[1], p[2]] for p in pairs)
+
+  # Return the result.
+  return MatchingResult(cost, mate)
+end
diff --git a/src/lp.jl b/src/lp.jl
index 472ea73..7e1abc9 100644
--- a/src/lp.jl
+++ b/src/lp.jl
@@ -1,30 +1,8 @@
-"""
-    maximum_weight_maximal_matching(g, w::Dict{Edge,Real})
-    maximum_weight_maximal_matching(g, w::Dict{Edge,Real}, cutoff)
-
-Given a bipartite graph `g` and an edgemap `w` containing weights associated to edges,
-returns a matching with the maximum total weight among the ones containing the
-greatest number of edges.
-
-Edges in `g` not present in `w` will not be considered for the matching.
-
-The algorithm relies on a linear relaxation on of the matching problem, which is
-guaranteed to have integer solution on bipartite graps.
-
-Eventually a `cutoff` argument can be given, to reduce computational times
-excluding edges with weights lower than the cutoff.
-
-The package JuMP.jl and one of its supported solvers is required.
-
-The returned object is of type `MatchingResult`.
-"""
-function maximum_weight_maximal_matching end
-
-function maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{U}, cutoff::R) where {U<:Real, R<:Real}
-    return maximum_weight_maximal_matching(g, solver, cutoff_weights(w, cutoff))
+function maximum_weight_maximal_matching_lp(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
+    return maximum_weight_maximal_matching_lp(g, solver, cutoff_weights(w, cutoff))
 end
 
-function maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{U}) where {U<:Real}
+function maximum_weight_maximal_matching_lp(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}) where {T<:Real}
 # TODO support for graphs with zero degree nodes
 # TODO apply separately on each connected component
     bpmap = bipartite_map(g)
@@ -88,7 +66,7 @@ function maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolve
 
     mate = fill(-1, nv(g))
     for e in edges(g)
-        if w[src(e),dst(e)] > zero(U)
+        if w[src(e),dst(e)] > zero(T)
             inmatch = convert(Bool, sol[edgemap[e]])
             if inmatch
                 mate[src(e)] = dst(e)
@@ -99,18 +77,3 @@ function maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolve
 
     return MatchingResult(cost, mate)
 end
-
-"""
-    cutoff_weights copies the weight matrix with all elements below cutoff set to 0
-"""
-function cutoff_weights(w::AbstractMatrix{U}, cutoff::R) where {U<:Real, R<:Real}
-    wnew = copy(w)
-    for j in 1:size(w,2)
-        for i in 1:size(w,1)
-            if wnew[i,j] < cutoff
-                wnew[i,j] = zero(U)
-            end
-        end
-    end
-    wnew
-end
diff --git a/src/maximum_weight_maximal_matching.jl b/src/maximum_weight_maximal_matching.jl
new file mode 100644
index 0000000..b368864
--- /dev/null
+++ b/src/maximum_weight_maximal_matching.jl
@@ -0,0 +1,99 @@
+"""
+    AbstractMaximumWeightMaximalMatchingAlgorithm
+Abstract type that allows users to pass in their preferred algorithm
+"""
+abstract type AbstractMaximumWeightMaximalMatchingAlgorithm end
+
+"""
+    LPAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm
+Forces the maximum_weight_maximal_matching function to use a linear programming formulation.
+"""
+struct LPAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm end
+
+function maximum_weight_maximal_matching(
+    g::Graph, 
+    w::AbstractMatrix{T}, 
+    algorithm::LPAlgorithm, 
+    solver = nothing
+) where {T<:Real}
+    if ! isa(solver, AbstractMathProgSolver)
+        error("The keyword argument solver must be an AbstractMathProgSolver, as accepted by JuMP.")
+    end
+
+    return maximum_weight_maximal_matching_lp(g, solver, w)
+end
+
+"""
+    HungarianAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm
+Forces the maximum_weight_maximal_matching function to use the Hungarian algorithm.
+"""
+struct HungarianAlgorithm <: AbstractMaximumWeightMaximalMatchingAlgorithm end
+
+function maximum_weight_maximal_matching(
+    g::Graph, 
+    w::AbstractMatrix{T}, 
+    algorithm::HungarianAlgorithm, 
+    solver = nothing
+) where {T<:Real}
+    return maximum_weight_maximal_matching_hungarian(g, w)
+end
+
+"""
+    maximum_weight_maximal_matching{T<:Real}(g, w::Dict{Edge,T})
+
+Given a bipartite graph `g` and an edge map `w` containing weights associated to edges,
+returns a matching with the maximum total weight among the ones containing the
+greatest number of edges.
+
+Edges in `g` not present in `w` will not be considered for the matching.
+
+A `cutoff` keyword argument can be given, to reduce computational times
+excluding edges with weights lower than the cutoff.
+
+Finally, a specific algorithm can be chosen (`algorithm` keyword argument); 
+each algorithm has specific dependencies. For instance: 
+
+- If `algorithm=HungarianAlgorithm()` (the default), the package Hungarian.jl is used. 
+  This algorithm is always polynomial in time, with complexity O(n³). 
+- If `algorithm=LPAlgorithm()`, the package JuMP.jl and one of its supported solvers is required. 
+  In this case, the algorithm relies on a linear relaxation on of the matching problem, which is
+  guaranteed to have integer solution on bipartite graphs. A solver must be provided with 
+  the `solver` keyword parameter. 
+
+The returned object is of type `MatchingResult`.
+"""
+function maximum_weight_maximal_matching(
+        g::Graph, 
+        w::AbstractMatrix{T}; 
+        cutoff = nothing,
+        algorithm::AbstractMaximumWeightMaximalMatchingAlgorithm = HungarianAlgorithm(), 
+        solver = nothing
+    ) where {T<:Real}
+
+    if cutoff != nothing && ! isa(cutoff, Real)
+        error("The cutoff value must be of type Real or nothing.")
+    end
+
+    if cutoff != nothing
+        return maximum_weight_maximal_matching(g, cutoff_weights(w, cutoff), algorithm, solver)
+    else
+        return maximum_weight_maximal_matching(g, w, algorithm, solver)
+    end
+end
+
+"""
+    cutoff_weights copies the weight matrix with all elements below cutoff set to 0
+"""
+function cutoff_weights(w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
+    wnew = copy(w)
+    for j in 1:size(w,2)
+        for i in 1:size(w,1)
+            if wnew[i,j] < cutoff
+                wnew[i,j] = zero(T)
+            end
+        end
+    end
+    wnew
+end
+
+@deprecate maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real} maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), cutoff=cutoff, solver=solver)
diff --git a/test/runtests.jl b/test/runtests.jl
index 43fdcf1..fa0853f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,20 +6,6 @@ using LinearAlgebra: I
 
 @testset "LightGraphsMatching" begin
 
-g = CompleteGraph(4)
-w = LightGraphsMatching.default_weights(g)
-@test all((w + w') .≈ ones(4,4) - Matrix(I, 4,4))
-
-w1 = [
-      1 3
-      5 1
-     ]
-w0 = [
-      0 3
-      5 0
-     ]
-@test all(w0 .≈ LightGraphsMatching.cutoff_weights(w1, 2))
-
 @testset "maximum_weight_matching" begin
     g = CompleteGraph(3)
     w = [
@@ -105,7 +91,7 @@ end
     w[1,4] = 1.
     w[2,3] = 2.
     w[2,4] = 11.
-    match = maximum_weight_maximal_matching(g, CbcSolver(), w)
+    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=CbcSolver())
     @test match.weight ≈ 21
     @test match.mate[1] == 3
     @test match.mate[3] == 1
@@ -118,7 +104,7 @@ end
     w[1,4] = 0.5
     w[2,3] = 11
     w[2,4] = 1
-    match = maximum_weight_maximal_matching(g, CbcSolver(), w)
+    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=CbcSolver())
     @test match.weight ≈ 11.5
     @test match.mate[1] == 4
     @test match.mate[4] == 1
@@ -133,7 +119,7 @@ end
     w[2,4] = 1
     w[2,5] = -1
     w[2,6] = -1
-    match = maximum_weight_maximal_matching(g,CbcSolver(),w,0)
+    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=CbcSolver(), cutoff=0)
     @test match.weight ≈ 11.5
     @test match.mate[1] == 4
     @test match.mate[4] == 1
@@ -148,7 +134,7 @@ end
     w[1,6] = 1
     w[1,5] = -1
 
-    match = maximum_weight_maximal_matching(g,CbcSolver(),w,0)
+    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=CbcSolver(), cutoff=0)
     @test match.weight ≈ 12
     @test match.mate[1] == 6
     @test match.mate[2] == 5
@@ -156,6 +142,45 @@ end
     @test match.mate[4] == -1
     @test match.mate[5] == 2
     @test match.mate[6] == 1
+    
+    
+    g = CompleteBipartiteGraph(2, 2)
+    w = zeros(4, 4)
+    w[1, 3] = 10.
+    w[1, 4] = 1.
+    w[2, 3] = 2.
+    w[2, 4] = 11.
+    match = maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
+    @test match.weight ≈ 21
+    @test match.mate[1] == 3
+    @test match.mate[3] == 1
+    @test match.mate[2] == 4
+    @test match.mate[4] == 2
+
+    g = CompleteGraph(3)
+    w = zeros(3, 3)
+    @test ! is_bipartite(g)
+    @test_throws ErrorException maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
+
+    g = CompleteBipartiteGraph(2, 4)
+    w = zeros(6, 6)
+    w[1, 3] = 10
+    w[1, 4] = 0.5
+    w[2, 3] = 11
+    w[2, 4] = 1
+    match = maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
+    @test match.weight ≈ 11.5
+
+    g = Graph(4)
+    add_edge!(g, 1, 3)
+    add_edge!(g, 1, 4)
+    add_edge!(g, 2, 4)
+    w = zeros(4, 4)
+    w[1, 3] = 1
+    w[1, 4] = 3 
+    w[2, 4] = 1
+    match = maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
+    @test match.weight ≈ 2 
 
 end
 
@@ -214,8 +239,8 @@ end
     @test match.weight ≈ -11.5
 
 
-    g =CompleteGraph(4)
-    w =Dict{Edge,Float64}()
+    g = CompleteGraph(4)
+    w = Dict{Edge,Float64}()
     w[Edge(1,3)] = 10
     w[Edge(1,4)] = 0.5
     w[Edge(2,3)] = 11