diff --git a/deps/build.jl b/deps/build.jl
index cc11739..64e376c 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -10,7 +10,7 @@ struct HSLVersion
     ext::String
 end
 
-getname(ver::HSLVersion) = ver.algname * "-" * ver.version * ver.ext
+getname(ver::HSLVersion, ext=true) = ver.algname * "-" * ver.version * (ext ? ver.ext : "")
 
 function checksha(version::HSLVersion, filepath)
     if isfile(filepath)
@@ -56,19 +56,29 @@ const hsl_ma57_patch = joinpath(hsl_ma57_path, "get_factors.patch")
 ##############################
 # MA97
 ##############################
-const hsl_ma97_verions = [
+const hsl_ma97_versions = [
     HSLVersion("hsl_ma97", "2.7.0", "ac3a081d3a28e9ecb8871ce769f4ced2a5ffa5a9c36defbd2c844ae3493ccb37", ".tar.gz"),
     HSLVersion("hsl_ma97", "2.7.0", "8221b607d96554d7a57cc60483c7305ef43a8785dc4171ac2e8da087900a1100", ".zip"),
     HSLVersion("hsl_ma97", "2.6.0", "be5fe822674be93e3d2e1a7d7ed6c5ad831b91cf8ca5150beb473f67af5fcb66", ".tar.gz"),
 ]
 const hsl_ma97_path = haskey(ENV, "HSL_MA97_PATH") ? ENV["HSL_MA97_PATH"] : joinpath(@__DIR__, "downloads")
-hsl_ma97_version = findversion(hsl_ma97_verions, hsl_ma97_path)
+hsl_ma97_version = findversion(hsl_ma97_versions, hsl_ma97_path)
 const hsl_ma97_archive = isnothing(hsl_ma97_version) ? "" : joinpath(hsl_ma97_path, getname(hsl_ma97_version))
 
+##############################
+# MA86
+##############################
+const hsl_ma86_versions = [
+  HSLVersion("hsl_ma86", "1.7.0",  "845c65bee0bf31507d7b99c87773997012b7b16434da349ad5c361ceb257191e", ".tar.gz"),
+]
+const hsl_ma86_path = haskey(ENV, "HSL_MA86_PATH") ? ENV["HSL_MA86_PATH"] : joinpath(@__DIR__, "downloads")
+hsl_ma86_version = findversion(hsl_ma86_versions, hsl_ma86_path)
+const hsl_ma86_archive = isnothing(hsl_ma86_version) ? "" : joinpath(hsl_ma86_path, getname(hsl_ma86_version))
+
 ##############################
 # Build
 ##############################
-const hsl_archives = [hsl_ma57_archive, hsl_ma97_archive]
+const hsl_archives = [hsl_ma57_archive, hsl_ma97_archive, hsl_ma86_archive]
 
 const HSL_FC = haskey(ENV, "HSL_FC") ? ENV["HSL_FC"] : "gfortran"
 const HSL_F77 = haskey(ENV, "HSL_F77") ? ENV["HSL_F77"] : HSL_FC
@@ -99,6 +109,12 @@ if any(isfile.(hsl_archives))
     include("build_hsl_ma57.jl")
   end
 
+  if isfile(hsl_ma86_archive)
+    @info "building ma86"
+    push!(products, FileProduct(prefix, "lib/libhsl_ma86.$so", :libhsl_ma86))
+    include("build_hsl_ma86.jl")
+  end
+
   if isfile(hsl_ma97_archive)
     @info "building ma97"
     push!(products, FileProduct(prefix, "lib/libhsl_ma97.$so", :libhsl_ma97))
diff --git a/deps/build_hsl_ma86.jl b/deps/build_hsl_ma86.jl
new file mode 100644
index 0000000..4b31b9e
--- /dev/null
+++ b/deps/build_hsl_ma86.jl
@@ -0,0 +1,12 @@
+version = hsl_ma86_version
+archive = hsl_ma86_archive
+if version.ext == ".tar.gz"
+    run(`tar -zxf $archive -C $builddir`)
+elseif hsl_ma97_version.ext == ".zip"
+    run(`unzip -o $archive -d $builddir`)
+end
+name = getname(version, false)
+cd("$builddir/$name")
+run(`./configure --prefix=$usrdir FC=$(HSL_FC) F77=$(HSL_F77) CC=$(HSL_CC) CFLAGS=-fPIC FFLAGS="-fPIC -fopenmp" FCFLAGS="-fPIC -fopenmp" --with-blas="-L$libopenblas_dir -lopenblas" --with-metis="-L$libmetis_dir -lmetis"`)
+run(`make install`)
+cd(@__DIR__)
diff --git a/src/HSL.jl b/src/HSL.jl
index bd25aa7..7249f89 100644
--- a/src/HSL.jl
+++ b/src/HSL.jl
@@ -24,6 +24,7 @@ const data_map = Dict{Type, Type}(Float32 => Cfloat,
                                   ComplexF32 => Cfloat,
                                   ComplexF64 => Cdouble)
 
+
 # package-specific definitions
 if (@isdefined libhsl_ma57) || haskey(ENV, "DOCUMENTER_KEY")
   include("hsl_ma57.jl")
@@ -34,5 +35,8 @@ end
 if (@isdefined libhsl_ma97) || haskey(ENV, "DOCUMENTER_KEY")
   include("hsl_ma97.jl")
 end
+if (@isdefined libhsl_ma86) || haskey(ENV, "DOCUMENTER_KEY")
+  include("hsl_ma86.jl")
+end
 
 end
diff --git a/src/hsl_ma86.jl b/src/hsl_ma86.jl
new file mode 100644
index 0000000..f54bbad
--- /dev/null
+++ b/src/hsl_ma86.jl
@@ -0,0 +1,431 @@
+export Ma86, Ma86_Control, Ma86_Info
+
+const Ma86Data = Union{Float32, Float64, ComplexF32, ComplexF64}
+const Ma86Real = Union{Cfloat, Cdouble}
+
+function appendtype(fname, T)
+  typesuffix = Dict(Float32=>"s", Float64=>"d", ComplexF32=>"c", ComplexF64=>"z")
+  return string(fname) * "_" * typesuffix[T]
+end
+
+"""
+    Ma86_Control
+
+A simple wrapper around the C control structure for C HSL Ma86 interface. Consult HSL 
+  documentation for more details.
+"""
+mutable struct Ma86_Control{S <: Ma86Real}
+  f_arrays::Cint
+  diagnostics_level::Cint
+  unit_diagnostics::Cint
+  unit_error::Cint
+  unit_warning::Cint
+  nemin::Cint
+  nb::Cint
+  action::Cint
+  nbi::Cint
+  bool_size::Cint
+  small_::S
+  static_::S
+  u::S
+  umin::S
+  scaling::Cint
+  function Ma86_Control(::Type{T}) where T
+    S = data_map[T]
+    control = new{S}(0,0,0,0,0,0,0,0,0,0,zero(S),zero(S),zero(S),zero(S),0)
+    ma86_default_control(T,control)
+    control.f_arrays = 1  # Use 1-based indexing for arrays, avoiding copies.
+    return control
+  end
+end
+Ma86_Control{S}() where S <: Ma86Real = Ma86_Control(S)
+
+"""
+    Ma86_Info
+
+A simple wrapper around the C info structure for C HSL Ma86 interface. Consult HSL 
+  documentation for more details.
+"""
+mutable struct Ma86_Info{S <: Ma86Real}
+  detlog::S
+  detsign::Cint
+  flag::Cint
+  matrix_rank::Cint
+  maxdepth::Cint
+  num_delay::Cint
+  num_factor::Clong
+  num_flops::Clong
+  num_neg::Cint
+  num_nodes::Cint
+  num_nothresh::Cint
+  num_perturbed::Cint
+  num_two::Cint
+  pool_size::Cint
+  stat::Cint
+  usmall::S
+  function Ma86_Info(::Type{T}) where T <: Ma86Data
+    S = data_map[T]
+    new{S}(
+      zero(S), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Clong), 
+      zero(Clong), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), 
+      zero(Cint), zero(S)
+    )
+  end
+end
+Ma86_Info{S}() where S <: Ma86Real = Ma86_Info(S)
+
+##############################
+# Convenience Wrapper Types 
+##############################
+struct Keep
+  ptr::Vector{Ptr{Cvoid}}
+end
+@inline Base.cconvert(::Type{Ref{Ptr{Cvoid}}}, akeep::Keep) = akeep.ptr
+@inline Base.unsafe_convert(::Type{Ptr{Ptr{Cvoid}}}, akeep::Keep) = 
+  Base.unsafe_convert(Ptr{Ptr{Cvoid}}, akeep.ptr)
+Keep() = Keep([C_NULL])
+isnull(akeep::Keep) = akeep.ptr[1] == C_NULL
+
+const matrix_types86 = (herm_indef = Cint(-4), cmpl_indef = Cint(-5))
+const jobs86 = (
+  A    = Cint(0),  # Solve Ax = b
+  PL   = Cint(1),  # Solve PLx = b
+  D    = Cint(2),  # Solve Dx = b
+  PLt  = Cint(3),  # Solve (PL)ᴴ x = b
+  DPLt = Cint(4),  # Solve D(PL)ᴴ x = b
+)
+
+Base.@kwdef mutable struct Ma86_Flags
+  isanalysisdone::Bool = false
+  isfactordone::Bool = false
+end
+
+"""
+    Ma86{T,S}
+
+Type for the HSL Ma86 solver, a multithreaded solver for symmetric indefinite or 
+symmetric positive definite systems.
+
+# Constructors
+
+    Ma86(colptr, rowval, nzval; [control, info, order, scale, doanalyse])
+
+Creates a new Ma86 solver with the lower-triangular matrix data stored in CSC format. 
+Both `colptr` and `rowval` will be converted to vectors of `Int32`. The element type of 
+the nonzero values `nzval` will determine the data type used by the solver. The data types
+of `control`, `info`, and `scale` should match the corresponding real data type.
+
+    Ma86(A::SparseMatrix; kwargs...)
+
+Extracts the CSC data from the sparse matrix `A`. If `A` is not lower triangular, the lower
+triangular data will be pulled out automatically. The keyword arguments are identical those
+of the previous constructor.
+
+# Usage
+With a constructed `Ma86` object, you can use any of the following methods to solve the 
+linear system `Ax = b` where `b` can have any number of columns. 
+
+## 1. Call HSL steps independently
+This mirrors the C interface of Ma86. Call the following methods in the specified order:
+
+    HSL.factor!(ma86)
+    HSL.solve!(ma86, b)
+
+Or
+
+    HSL.factorsolve!(ma86, b)
+
+Note the methods above will overwrite the data in `b` with the solution `x`. Alternatively,
+you can use [`solve`](@ref) and ][`factorsolve`](@ref) which will return the solution 
+vector `x` and leave `b` unchanged.
+
+By default, the constructor will run the [`analyse!`](@ref) method on the provided sparsity 
+structure. This step can be skipped by passing `doanalyse=false` to the constructor. The 
+analysis step must then be called manually before calling either `factor!` or 
+`factorsolve!`.
+
+## 2. Use backslash or `ldiv!`
+To return the solution matrix, use `x = ma86 \\ b`, or to overwrite the `b` matrix use 
+`ldiv!(ma86, b)`.
+"""
+mutable struct Ma86{T<:Ma86Data, S<:Ma86Real}
+  __keep::Keep
+  n::Cint
+  colptr::Vector{Cint}
+  rowval::Vector{Cint}
+  nzval::Vector{T}
+  control::Ma86_Control{S}
+  info::Ma86_Info{S}
+  order::Vector{Cint}
+  scale::Vector{S}
+  flags::Ma86_Flags
+  function Ma86(colptr::Vector{<:Integer}, rowval::Vector{<:Integer}, nzval::Vector{T}, 
+                control::Ma86_Control{S}, info::Ma86_Info{S}, order::Vector{<:Integer}, 
+                scale::Vector{<:Real}; doanalyse::Bool=true) where {T <: Ma86Data, S}
+    S == data_map[T] || throw(TypeError(:Ma86, "MA86{$T, $S}\n", data_map[T], T))
+    keep = Keep()
+    n = length(colptr) - 1
+    @assert length(order) == n
+    @assert length(scale) == n
+    flags = Ma86_Flags()
+    ma86 = new{T,S}(keep, n, colptr, rowval, nzval, control, info, order, scale, flags)
+    if doanalyse
+      analyse!(ma86)
+    end
+    finalizer(_finalize, ma86)
+  end
+end
+function Ma86(colptr::Vector{<:Integer}, rowval::Vector{<:Integer}, nzval::Vector{T}; 
+              control::Ma86_Control=Ma86_Control(T), 
+              info::Ma86_Info=Ma86_Info(T), 
+              order::Vector{<:Integer}=collect(1:length(colptr)-1),
+              scale::Vector{<:Real}=fill(one(data_map[T]), length(colptr)-1),
+              kwargs...) where T
+  # Assumes the data is already lower triangular
+  Ma86(colptr, rowval, nzval, control, info, order, scale; kwargs...)
+end
+function Ma86(A::SparseMatrixCSC; kwargs...)
+  if !istril(A)
+    A = tril(A)
+  end
+  Ma86(A.colptr, A.rowval, A.nzval; kwargs...)
+end
+datatype(::Ma86{T}) where T = T
+realtype(::Ma86{<:Any,S}) where S = S
+
+function _finalize(m::Ma86{T}) where T
+  if !isnull(m.__keep)
+    ma86_finalise(T, m.__keep, m.control)
+  end
+  m.flags.isanalysisdone = false
+  m.flags.isfactordone = false
+  m.__keep.ptr[1] = C_NULL
+  return m
+end
+
+##############################
+# Exception Handling
+##############################
+"""
+   Ma86Exception
+
+Exception thrown by HSL MA86 algorithm. The `flag` field contains the HSL flag 
+for the error. Consult the documentation for details of the exception.
+"""
+mutable struct Ma86Exception <: Exception
+  msg::AbstractString
+  flag::Int
+end
+
+function checkerror(m::Ma86, msg)
+  if m.info.flag != 0
+    _finalize(m)  # free memory in keep before throwing an error
+    throw(Ma86Exception(msg, m.info.flag))
+  end
+end
+
+##############################
+# Convenience Wrapper
+##############################
+
+function analyse!(m::Ma86{T}) where T
+  ma86_analyse(T, m.n, m.colptr, m.rowval, m.order, m.__keep, m.control, m.info)
+  checkerror(m, "MA86: Error during symbol analysis.")
+  m.flags.isanalysisdone = true
+  return nothing
+end
+
+function factor!(m::Ma86{<:Ma86Data}; matrix_type::Integer=matrix_types86.herm_indef)
+  if isnull(m.__keep) || !m.flags.isanalysisdone
+    error("Cannot call factor! before analyse!")
+    return nothing
+  end
+  ma86_factor(
+    Cint(matrix_type), m.n, m.colptr, m.rowval, m.nzval, m.order, m.__keep, m.control, m.info, 
+    m.scale
+  )
+  checkerror(m, "MA86: Error during factorization.")
+  m.flags.isfactordone = true
+  return nothing
+end
+
+function solve!(m::Ma86{T}, x::VecOrMat{T}; job::Integer=jobs86.A) where T
+  if !m.flags.isanalysisdone || !m.flags.isfactordone
+    error("Cannot call solve! before factor!")
+    return x
+  end
+  size(x,1) != m.n && throw(DimensionMismatch("x must have $(m.n) rows, got $(size(x,1))."))
+  ldx = Cint(stride(x,2))
+  nrhs = Cint(size(x,2))
+  ma86_solve(Cint(job), nrhs, ldx, x, m.order, m.__keep, m.control, m.info)
+  checkerror(m, "MA86: Error during solve.")
+  return x
+end
+solve(m::Ma86{T}, b::VecOrMat{T}; kwargs...) where T = solve!(m, copy(b); kwargs...)
+
+function factorsolve!(m::Ma86{T}, x::VecOrMat{T}; 
+  matrix_type::Integer=matrix_types86.herm_indef
+) where T <: Ma86Data
+  if isnull(m.__keep) || !m.flags.isanalysisdone
+    error("Cannot call factorsolve! before analyse!")
+    return x 
+  end
+  size(x,1) != m.n && throw(DimensionMismatch("x must have $(m.n) rows, got $(size(x,1))."))
+  ldx = Cint(stride(x,2))
+  nrhs = Cint(size(x,2))
+  ma86_factor_solve(
+    Cint(matrix_type), m.n, m.colptr, m.rowval, m.nzval, m.order, m.__keep, m.control, m.info, 
+    nrhs, ldx, x, m.scale
+  )
+  checkerror(m, "MA86: Error during factorsolve.")
+  m.flags.isfactordone = true
+  return x
+end
+
+# Out-of-place version
+factorsolve(m::Ma86{T}, b::VecOrMat{T}; kwargs...) where T <: Ma86Data = 
+  factorsolve!(m, copy(b); kwargs...)
+
+# LinearAlgebra wrappers
+function Base.:\(m::Ma86{T}, b::VecOrMat{T}) where T <: Ma86Data 
+  !m.flags.isanalysisdone && analyse!(m)
+  m.flags.isfactordone && return solve(m, b)
+  factorsolve(m, b)
+end
+
+function LinearAlgebra.ldiv!(m::Ma86{T}, x::VecOrMat{T}) where T <: Ma86Data
+  !m.flags.isanalysisdone && analyse!(m)
+  m.flags.isfactordone && return solve!(m, x)
+  factorsolve!(m, x)
+end
+
+##############################
+# C Wrapper
+##############################
+
+for T in (Float32, Float64, ComplexF32, ComplexF64) 
+  S = data_map[T]
+
+  @eval function ma86_default_control(::Type{$T}, control::Ma86_Control{$S})
+    ccall(($(appendtype(:ma86_default_control, T)), libhsl_ma86), Cvoid,
+      (Ref{Ma86_Control{$S}},),
+      control
+    )
+  end
+
+  # Add type as argument to avoid ambiguity
+  @eval function ma86_analyse(::Type{$T}, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+                              order::Vector{Cint}, keep::Keep, control::Ma86_Control, 
+                              info::Ma86_Info)
+    @assert length(ptr) == n+1
+    @assert ptr[end] == length(row) + 1
+    @assert length(order) == n
+    ccall(($(appendtype(:ma86_analyse, T)), libhsl_ma86), Cvoid,
+      (Cint, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}}, 
+      Ref{Ma86_Info{$S}}),
+      n, ptr, row, order, keep, control, info
+    )
+  end
+
+  if T <: Ma86Real
+    # Note that matrix_type argument is unused in these methods
+    # Included to make the signature match for all data types
+    @eval function ma86_factor(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
+                              order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                              info::Ma86_Info{$S}, scale::Vector{$S})
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      ccall(($(appendtype(:ma86_factor, T)), libhsl_ma86), Cvoid,
+        (Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Ref{$S}),
+        n, ptr, row, val, order, keep, control, info, scale
+      )
+    end
+
+    @eval function ma86_factor_solve(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+                                    val::Vector{$T}, order::Vector{Cint}, keep::Keep, 
+                                    control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
+                                    nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                                    scale::Vector{$S})
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      @assert nrhs >= 1
+      @assert ldx >= n
+      @assert length(scale) == n
+      @assert size(x,1) == n
+      @assert size(x,2) == nrhs
+      ccall(($(appendtype(:ma86_factor_solve, T)), libhsl_ma86), Cvoid,
+        (Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Cint, Cint, Ref{$T}, Ref{$S}),
+        n, ptr, row, val, order, keep, control, info, nrhs, ldx, x, scale
+      )
+    end
+  else
+    @eval function ma86_factor(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, 
+                               row::Vector{Cint}, val::Vector{$T}, 
+                               order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                               info::Ma86_Info{$S}, scale::Vector{$S})
+      @assert matrix_type ∈ (-4,-5)
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      ccall(($(appendtype(:ma86_factor, T)), libhsl_ma86), Cvoid,
+        (Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Ref{$S}),
+        matrix_type, n, ptr, row, val, order, keep, control, info, scale
+      )
+    end
+
+    @eval function ma86_factor_solve(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, 
+                                     row::Vector{Cint}, val::Vector{$T}, 
+                                     order::Vector{Cint}, keep::Keep, 
+                                     control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
+                                     nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                                     scale::Vector{$S})
+      @assert matrix_type ∈ (-4,-5)
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      @assert nrhs >= 1
+      @assert ldx >= n
+      @assert length(scale) == n
+      @assert size(x,1) == n
+      @assert size(x,2) == nrhs
+      ccall(($(appendtype(:ma86_factor_solve, T)), libhsl_ma86), Cvoid,
+        (Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Cint, Cint, Ref{$T}, Ref{$S}),
+        matrix_type, n, ptr, row, val, order, keep, control, info, nrhs, ldx, x, scale
+      )
+    end
+  end
+
+  @eval function ma86_solve(job::Cint, nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                          order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                          info::Ma86_Info{$S})
+    n = size(x,1)
+    @assert job ∈ 0:4
+    @assert nrhs >= 1
+    @assert ldx >= n
+    @assert size(x,2) == nrhs
+    ccall(($(appendtype(:ma86_solve, T)), libhsl_ma86), Cvoid,
+      (Cint, Cint, Cint, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}},
+      Ref{Ma86_Info{$S}}, Ptr{$S}),
+      job, nrhs, ldx, x, order, keep, control, info, C_NULL 
+    )
+  end
+
+  # Add type as argument to avoid ambiguity
+  @eval function ma86_finalise(::Type{$T}, keep::Keep, control::Ma86_Control{$S})
+    ccall(($(appendtype(:ma86_finalise, T)), libhsl_ma86), Cvoid,
+      (Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}}),
+      keep, control
+    )
+  end
+
+end
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..59b801f
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index 1fd4aee..667edea 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,14 +4,17 @@ using LinearAlgebra
 using SparseArrays
 using Test
 
-include("../deps/deps.jl")
+# include("../deps/deps.jl")
 
-if @isdefined libhsl_ma57
+if isdefined(HSL, :libhsl_ma57)
   include("test_ma57.jl")
-  if @isdefined libhsl_ma57_patch
+  if isdefined(HSL, :libhsl_ma57_patch)
     include("test_ma57_patch.jl")
   end
 end
-if @isdefined libhsl_ma97
+if isdefined(HSL, :libhsl_ma97)
   include("test_ma97.jl")
 end
+if isdefined(HSL, :libhsl_ma86)
+  include("test_ma86.jl")
+end
diff --git a/test/test_ma86.jl b/test/test_ma86.jl
new file mode 100644
index 0000000..52b8e97
--- /dev/null
+++ b/test/test_ma86.jl
@@ -0,0 +1,178 @@
+using HSL
+using SparseArrays
+using LinearAlgebra
+using Test
+using Random
+
+function getallocs(A,b)
+    T = eltype(A)
+    n = size(A,1)
+    ma86 = Ma86(A)
+    allocs = @allocated HSL.analyse!(ma86)
+    allocs += @allocated HSL.factor!(ma86)
+    x = copy(b)
+    allocs += @allocated HSL.solve!(ma86, x)
+end
+
+T = Float64
+n = 10
+A = convert(T, 3) * sprand(T, n, n, 0.5)
+Aindef = A + A' + one(T)*I
+ma86 = Ma86(Aindef)
+ma86.flags.isanalysisdone
+HSL.factor!(ma86)
+b = rand(T,n)
+HSL.solve!(ma86, b)
+
+@testset "MA86 ($T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+    # @info "Testing hsl_ma86 with $T data"
+    n = 10
+    A = convert(T, 3) * sprand(T, n, n, 0.5)
+    Aindef = A + A' + one(T)*I
+
+    Apsd = A'A + one(T)*I
+    @test isposdef(Matrix(Apsd))
+
+
+    # Test symmetric indefinite A
+    @testset "$label" for (label,A) in (("Indefinite", Aindef), ("Positive Definite",Apsd))
+        n = size(A,1)
+        T = eltype(A)
+
+        # Test constructors
+        L = SparseMatrixCSC{T,Cint}(tril(A))
+        let ptr = L.colptr, row = L.rowval, val = L.nzval
+            # Check that new arrays aren't allocated if the input in Cint
+            let ma86 = Ma86(ptr, row, val)
+                @test ma86.colptr === ptr
+                @test ma86.rowval === row
+                @test ma86.nzval === val 
+                @test ma86.flags.isanalysisdone == true
+                @test !HSL.isnull(ma86.__keep)
+            end
+
+            # Check the keyword constructor            
+            control = Ma86_Control(T)
+            info = Ma86_Info(T)
+            S = HSL.data_map[T]
+            order = randperm(Cint(n))
+            scale = fill(one(S), n)
+            let ma86 = HSL.Ma86(ptr, row, val, 
+                    control=control, info=info, order=order, scale=scale)
+                @test ma86.order === order
+                @test ma86.scale === scale 
+                @test ma86.control === control
+                @test ma86.info === info 
+            end
+
+            # Pass in Int64 data
+            let ma86 = HSL.Ma86(Int64.(ptr), Int64.(row), val, 
+                    order=randperm(Int64(n)), scale=Float64.(scale))
+                @test ma86.colptr !== ptr
+                @test ma86.rowval !== row
+                @test ma86.order !== order
+                @test ma86.scale !== scale 
+            end
+
+            # Pass in lower trianglar matrix
+            # Make sure it doesn't reallocate new data
+            let ma86 = HSL.Ma86(L) 
+                @test ma86.colptr === L.colptr
+                @test ma86.rowval === L.rowval
+                @test ma86.nzval === L.nzval
+            end
+
+            # Check analyse skipping
+            let ma86 = HSL.Ma86(L, doanalyse=false)
+                @test ma86.flags.isanalysisdone == false
+            end
+        end
+
+        @testset "with $(size(b,2)) rhs vectors" for b in (rand(T, n), rand(T, n, 3))
+            n = size(A,1)
+            T = eltype(A)
+
+            # factor, then solve
+            ma86_0 = Ma86(A)
+            let ma86 = ma86_0
+                x = copy(b)
+                HSL.factor!(ma86)
+                @test ma86.info.flag == 0
+                HSL.solve!(ma86, x)
+                @test ma86.info.flag == 0
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+                @test HSL.solve(ma86, b) ≈ x  # out-of-place solve
+                @test ma86 \ b ≈ x            # backslash
+                x2 = copy(b)
+                ldiv!(ma86, x2)
+                @test x2 ≈ x
+            end
+
+            # in-place
+            ma86_1 = Ma86(A)
+            let ma86 = ma86_1
+                x = copy(b)
+                HSL.factorsolve!(ma86, x)
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+            end
+
+            # out-of-place
+            ma86_2 = Ma86(A)
+            let ma86 = ma86_2
+                x = HSL.factorsolve(ma86, b)
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+            end
+
+            # Test linear algebra
+            ma86_3 = Ma86(A)
+            x = ma86_3 \ b
+            @test A*x ≈ b atol=sqrt(eps(norm(x)))
+
+            ma86_4 = Ma86(A)
+            x = copy(b)
+            ldiv!(ma86_4, x)
+            @test A*x ≈ b atol=sqrt(eps(norm(x)))
+
+            # Test errors
+            ma86_err = Ma86(Apsd, doanalyse=false)
+            let ma86 = ma86_err
+                @test_throws ErrorException HSL.factor!(ma86)
+                b = rand(T, n)
+                @test_throws ErrorException HSL.solve!(ma86, b)
+                @test_throws ErrorException HSL.factorsolve!(ma86, b)
+
+                HSL.analyse!(ma86)
+                @test_throws ErrorException HSL.solve!(ma86, b) 
+
+                b = rand(T, n)
+                x = copy(b)
+                HSL.factorsolve!(ma86, x)
+                @test Apsd*x ≈ b atol=sqrt(eps(norm(x)))
+
+                x = copy(b)
+                HSL.solve!(ma86, x)
+                @test Apsd*x ≈ b atol=sqrt(eps(norm(x)))
+
+                b_bad = rand(T, 2n)
+                @test_throws DimensionMismatch HSL.solve!(ma86, b_bad)
+                @test_throws DimensionMismatch HSL.factorsolve!(ma86, b_bad)
+
+                # Test allocations
+                @test getallocs(A, b) == 0
+            end
+            
+            # Test HSL Errors
+            # Prints out annoying warning message from MC78.
+            # Not sure how to avoid this.
+            # Should we even allow for non-square matrices to be passed to the constructor?
+            # B = sprand(T, n, n, 0.1)
+            # ma86_hslerr = Ma86(B, doanalyse=false)
+            # ma86_hslerr.control.diagnostics_level = -1
+            # ma86_hslerr.control.unit_diagnostics = -1
+            # ma86_hslerr.control.unit_error = -1
+            # ma86_hslerr.control.unit_warning = -1
+            # @test_throws HSL.Ma86Exception HSL.analyse!(ma86_hslerr)
+        end
+    end
+end
+
diff --git a/test/test_ma97.jl b/test/test_ma97.jl
index 1c9d5a7..17f7818 100644
--- a/test/test_ma97.jl
+++ b/test/test_ma97.jl
@@ -1,3 +1,6 @@
+using Random
+Random.seed!(1)
+
 function getallocs(A)
     T = eltype(A)
     n = size(A,1)