diff --git a/.gitignore b/.gitignore
index 63a7f65..8407416 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
 Manifest.toml
+build
 
 ## Core latex/pdflatex auxiliary files:
 *.aux
diff --git a/Project.toml b/Project.toml
index b962a11..74f9e98 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CatViews = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+Hiop = "f890763e-d4e5-4bf9-9dfa-7ed54d93ce96"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/build_hiop.sh b/build_hiop.sh
new file mode 100755
index 0000000..44a0bc2
--- /dev/null
+++ b/build_hiop.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+git submodule init
+git submodule update
+
+mkdir -p build
+cd build
+cmake -DHIOP_USE_MPI=OFF -DCMAKE_BUILD_TYPE=Debug -DHIOP_BUILD_SHARED=ON -DHIOP_BUILD_STATIC=ON -DCMAKE_INSTALL_PREFIX=$PWD ../hiop
+make -j
diff --git a/test/hiop.options b/test/hiop.options
new file mode 100755
index 0000000..1e98d61
--- /dev/null
+++ b/test/hiop.options
@@ -0,0 +1,13 @@
+dualsUpdateType linear
+dualsInitialization zero
+
+Hessian analytical_exact
+KKTLinsys auto
+compute_mode cpu
+
+verbosity_level 0
+mu0 1e-1
+fixed_var relax
+max_iter 3000
+kappa1 1e-8
+kappa2 1e-8
diff --git a/test/runtests.jl b/test/runtests.jl
index e3d6df6..2d74748 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -38,10 +38,14 @@ algparams = AlgParams()
 algparams.parallel = false #algparams.parallel = (nprocs() > 1)
 algparams.verbose = 0
 
+solver_list = ["Ipopt", "MadNLP"]
 if CUDA.has_cuda_gpu()
-    solver_list = ["Ipopt", "MadNLP", "MadNLPGPU"]
-else
-    solver_list = ["Ipopt", "MadNLP"]
+    push!(solver_list, "MadNLPGPU")
+end
+if isfile(joinpath(dirname(@__FILE__), "..", "build/libhiop.so"))
+    push!(solver_list, "Hiop")
+    ENV["JULIA_HIOP_LIBRARY_PATH"] = joinpath(dirname(@__FILE__), "..", "build") 
+    @info("Using Hiop at $(ENV["JULIA_HIOP_LIBRARY_PATH"])")
 end
 
 @testset "Test ProxAL on $(case)" begin
@@ -49,6 +53,7 @@ end
 
     for solver in solver_list
     @testset "$(solver)" begin
+        println("Testing using $(solver)")
         if solver == "Ipopt"
             using Ipopt
             algparams.optimizer =
@@ -99,7 +104,7 @@ end
                     result = solve_fullmodel(opfdata, rawdata, modelinfo, algparams)
                     @test isapprox(result["objective_value_nondecomposed"], 11.258316111585623, rtol = rtol)
                     @test isapprox(result["primal"].Pg[:], [0.8979870694509675, 1.3432060120295906, 0.9418738103137331, 0.9840203268625166, 1.448040098924617, 1.0149638876964715], rtol = rtol)
-                    if solver != "Ipopt"
+                    if solver == "MadNLP" || solver == "MadNLPGPU" || solver == "Hiop"
                         @test_broken isapprox(result["primal"].Zt[:], [0.0, 0.0, 0.0, 2.7859277234613066e-6, 2.3533760802049378e-6, 2.0234235436650152e-6], rtol = rtol)
                     else
                         @test isapprox(result["primal"].Zt[:], [0.0, 0.0, 0.0, 2.7859277234613066e-6, 2.3533760802049378e-6, 2.0234235436650152e-6], rtol = rtol)
@@ -148,19 +153,32 @@ end
                     algparams.mode = :nondecomposed
                     result = solve_fullmodel(opfdata, rawdata, modelinfo, algparams)
                     @test isapprox(result["objective_value_nondecomposed"], 11.258316111574212, rtol = rtol)
-                    @test isapprox(result["primal"].Pg[:], [0.8979870693416382, 1.3432060108971793, 0.9418738115511179, 0.9055318507524525, 1.3522597485901564, 0.9500221754747974, 0.9840203265549852, 1.4480400977338292, 1.014963889201792, 0.9932006221514175, 1.459056452449548, 1.024878608445939], rtol = rtol)
-                    if solver != "Ipopt"
+                    if solver == "Hiop"
+                        @test_broken isapprox(result["primal"].Pg[:], [0.8979870693416382, 1.3432060108971793, 0.9418738115511179, 0.9055318507524525, 1.3522597485901564, 0.9500221754747974, 0.9840203265549852, 1.4480400977338292, 1.014963889201792, 0.9932006221514175, 1.459056452449548, 1.024878608445939], rtol = rtol)
+                    else
+                        @test isapprox(result["primal"].Pg[:], [0.8979870693416382, 1.3432060108971793, 0.9418738115511179, 0.9055318507524525, 1.3522597485901564, 0.9500221754747974, 0.9840203265549852, 1.4480400977338292, 1.014963889201792, 0.9932006221514175, 1.459056452449548, 1.024878608445939], rtol = rtol)
+                    end
+                    if solver == "MadNLP" || solver == "MadNLPGPU" || solver == "Hiop"
                         @test_broken isapprox(result["primal"].Zt[:], [0.0, 0.0, 0.0, 2.7857429391709934e-6, 2.353279608425683e-6, 2.023313001658965e-6], rtol = rtol)
                     else
                         @test isapprox(result["primal"].Zt[:], [0.0, 0.0, 0.0, 2.7857429391709934e-6, 2.353279608425683e-6, 2.023313001658965e-6], rtol = rtol)
                     end
-                    @test isapprox(result["primal"].ωt[:], [0.0, -0.00012071650257302939, 0.0, -0.00014688472954291597], rtol = rtol)
+                    if solver == "Hiop"
+                        @test_broken isapprox(result["primal"].ωt[:], [0.0, -0.00012071650257302939, 0.0, -0.00014688472954291597], rtol = rtol)
+                    else
+                        @test isapprox(result["primal"].ωt[:], [0.0, -0.00012071650257302939, 0.0, -0.00014688472954291597], rtol = rtol)
+                    end
                 end
                 @testset "Lyapunov bound" begin
                     algparams.mode = :lyapunov_bound
                     result = solve_fullmodel(opfdata, rawdata, modelinfo, algparams)
-                    @test isapprox(result["objective_value_lyapunov_bound"], 11.258316111574207)
-                    @test isapprox(result["primal"].Pg[:], [0.8979870693416395, 1.3432060108971777, 0.9418738115511173, 0.9055318507524539, 1.3522597485901549, 0.9500221754747968, 0.9840203265549855, 1.4480400977338292, 1.0149638892017916, 0.9932006221514178, 1.459056452449548, 1.0248786084459387], rtol = rtol)
+                    if solver == "Hiop"
+                        @test_broken isapprox(result["objective_value_lyapunov_bound"], 11.258316111574207)
+                        @test_broken isapprox(result["primal"].Pg[:], [0.8979870693416395, 1.3432060108971777, 0.9418738115511173, 0.9055318507524539, 1.3522597485901549, 0.9500221754747968, 0.9840203265549855, 1.4480400977338292, 1.0149638892017916, 0.9932006221514178, 1.459056452449548, 1.0248786084459387], rtol = rtol)
+                    else
+                        @test isapprox(result["objective_value_lyapunov_bound"], 11.258316111574207)
+                        @test isapprox(result["primal"].Pg[:], [0.8979870693416395, 1.3432060108971777, 0.9418738115511173, 0.9055318507524539, 1.3522597485901549, 0.9500221754747968, 0.9840203265549855, 1.4480400977338292, 1.0149638892017916, 0.9932006221514178, 1.459056452449548, 1.0248786084459387], rtol = rtol)
+                    end
                 end
 
                 @testset "ProxALM" begin
@@ -201,22 +219,34 @@ end
                     algparams.mode = :nondecomposed
                     result = solve_fullmodel(opfdata, rawdata, modelinfo, algparams)
                     @test isapprox(result["objective_value_nondecomposed"], 11.258316111574212, rtol = rtol)
-                    @test isapprox(result["primal"].Pg[:], [0.8979870693416382, 1.3432060108971793, 0.9418738115511179, 0.9055318507524525, 1.3522597485901564, 0.9500221754747974, 0.9840203265549852, 1.4480400977338292, 1.014963889201792, 0.9932006221514175, 1.459056452449548, 1.024878608445939], rtol = rtol)
-                    if solver != "Ipopt"
+                    if solver == "Hiop"
+                        @test_broken isapprox(result["primal"].Pg[:], [0.8979870693416382, 1.3432060108971793, 0.9418738115511179, 0.9055318507524525, 1.3522597485901564, 0.9500221754747974, 0.9840203265549852, 1.4480400977338292, 1.014963889201792, 0.9932006221514175, 1.459056452449548, 1.024878608445939], rtol = rtol)
+                    else
+                        @test isapprox(result["primal"].Pg[:], [0.8979870693416382, 1.3432060108971793, 0.9418738115511179, 0.9055318507524525, 1.3522597485901564, 0.9500221754747974, 0.9840203265549852, 1.4480400977338292, 1.014963889201792, 0.9932006221514175, 1.459056452449548, 1.024878608445939], rtol = rtol)
+                    end
+                    if solver == "MadNLP" || solver == "MadNLPGPU" || solver == "Hiop"
                         @test_broken isapprox(result["primal"].Zt[:], [0.0, 0.0, 0.0, 2.7857429391709934e-6, 2.353279608425683e-6, 2.023313001658965e-6], rtol = rtol)
                     else
                         @test isapprox(result["primal"].Zt[:], [0.0, 0.0, 0.0, 2.7857429391709934e-6, 2.353279608425683e-6, 2.023313001658965e-6], rtol = rtol)
                     end
-                    @test isapprox(result["primal"].ωt[:], [0.0, -0.00012071650257302939, 0.0, -0.00014688472954291597], rtol = rtol)
+                    if solver == "Hiop"
+                        @test_broken isapprox(result["primal"].ωt[:], [0.0, -0.00012071650257302939, 0.0, -0.00014688472954291597], rtol = rtol)
+                    else
+                        @test isapprox(result["primal"].ωt[:], [0.0, -0.00012071650257302939, 0.0, -0.00014688472954291597], rtol = rtol)
+                    end
                 end
                 @testset "Lyapunov bound" begin
                     algparams.mode = :lyapunov_bound
                     result = solve_fullmodel(opfdata, rawdata, modelinfo, algparams)
-                    @test isapprox(result["objective_value_lyapunov_bound"], 11.258316111574207)
+                    if solver == "Hiop"
+                        @test_broken isapprox(result["objective_value_lyapunov_bound"], 11.258316111574207)
+                    else
+                        @test isapprox(result["objective_value_lyapunov_bound"], 11.258316111574207)
+                    end
                     @test isapprox(result["primal"].Pg[:], [0.8979870693416395, 1.3432060108971777, 0.9418738115511173, 0.9055318507524539, 1.3522597485901549, 0.9500221754747968, 0.9840203265549855, 1.4480400977338292, 1.0149638892017916, 0.9932006221514178, 1.459056452449548, 1.0248786084459387], rtol = rtol)
                 end
                 # MadNLP gets stuck here
-                if solver == "Ipopt"
+                if solver == "Ipopt" || solver == "Hiop"
                 @testset "ProxALM" begin
                     algparams.mode = :coldstart
                     runinfo = run_proxALM(opfdata, rawdata, modelinfo, algparams)
@@ -224,10 +254,15 @@ end
                     @test isapprox(runinfo.x.ωt[:], [0.0, -6.872285646588893e-5, 0.0, -8.763210220007218e-5], rtol = rtol)
                     @test isapprox(runinfo.λ.ramping[:], [0.0, 0.0, 0.0, -0.02526641008086374, -0.03138993075214526, -0.022582337855672596], rtol = rtol)
                     @test isapprox(runinfo.λ.ctgs[:], [0.0, 0.0, 0.0, 0.029837851198013726, 0.0018123712232017784, -0.029641691338447657, 0.0, 0.0, 0.0, 0.03640088659004858, -0.0024938071819119163, -0.030934017357459445], rtol = rtol)
-                    @test isapprox(runinfo.maxviol_c[end], 9.297964943270377e-5)
-                    @test isapprox(runinfo.maxviol_t[end], 3.7014553733172306e-6, rtol = rtol)
+                    if solver == "Hiop"
+                        @test_broken isapprox(runinfo.maxviol_c[end], 9.297964943270377e-5)
+                        @test_broken isapprox(runinfo.maxviol_t[end], 3.7014553733172306e-6, rtol = rtol)
+                    else
+                        @test isapprox(runinfo.maxviol_c[end], 9.297964943270377e-5)
+                        @test isapprox(runinfo.maxviol_t[end], 3.7014553733172306e-6, rtol = rtol)
+                    end
                     @test isapprox(runinfo.maxviol_d[end], 1.2901722926388947e-6, rtol = rtol)
-                    if solver == Ipopt
+                    if solver == "Ipopt"
                         @test runinfo.iter == 81
                     end
                 end