diff --git a/Project.toml b/Project.toml
index 60bf894..505e08c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,7 @@ authors = ["Benoît Legat <benoit.legat@gmail.com>"]
 version = "0.1.0"
 
 [deps]
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 SDPLR_jll = "3a057b76-36a0-51f0-a66f-6d580b8e8efd"
 
 [compat]
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
new file mode 100644
index 0000000..4e33b6b
--- /dev/null
+++ b/src/MOI_wrapper.jl
@@ -0,0 +1,442 @@
+# This file modifies code from SDPAFamily.jl (https://github.com/ericphanson/SDPAFamily.jl/), which is available under an MIT license (see LICENSE).
+
+import MathOptInterface as MOI
+
+const SupportedSets =
+    Union{MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle}
+
+mutable struct Optimizer <: MOI.AbstractOptimizer
+    objective_constant::Float64
+    objective_sign::Int
+    blksz::Vector{Cptrdiff_t}
+    blktype::Vector{Cchar}
+    varmap::Vector{Tuple{Int,Int,Int}} # Variable Index vi -> blk, i, j
+    b::Vector{Cdouble}
+    Cent::Vector{Cdouble}
+    Crow::Vector{Csize_t}
+    Ccol::Vector{Csize_t}
+    Cinfo_entptr::Vector{Csize_t}
+    Cinfo_type::Vector{Cchar}
+    Aent::Vector{Cdouble}
+    Arow::Vector{Csize_t}
+    Acol::Vector{Csize_t}
+    Ainfo_entptr::Vector{Vector{Csize_t}}
+    Ainfo_type::Vector{Vector{Cchar}}
+    silent::Bool
+    params::Parameters
+    Rmap::Vector{Int}
+    R::Union{Nothing,Vector{Cdouble}}
+    lambda::Vector{Cdouble}
+    ranks::Vector{Csize_t}
+    pieces::Union{Nothing,Vector{Cdouble}}
+    function Optimizer()
+        return new(
+            0.0,
+            1,
+            Cptrdiff_t[],
+            Cchar[],
+            Tuple{Int,Int,Int}[],
+            Cdouble[],
+            Cdouble[],
+            Csize_t[],
+            Csize_t[],
+            Csize_t[],
+            Cchar[],
+            Cdouble[],
+            Csize_t[],
+            Csize_t[],
+            Vector{Csize_t}[],
+            Vector{Cchar}[],
+            false,
+            Parameters(),
+            Int[],
+            nothing,
+            Cdouble[],
+            Csize_t[],
+            nothing,
+        )
+    end
+end
+
+function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute)
+    return hasfield(Parameters, Symbol(param.name))
+end
+function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
+    if !MOI.supports(optimizer, param)
+        throw(MOI.UnsupportedAttribute(param))
+    end
+    setfield!(optimizer.params, Symbol(param.name), value)
+    return
+end
+function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
+    if !MOI.supports(optimizer, param)
+        throw(MOI.UnsupportedAttribute(param))
+    end
+    getfield!(optimizer.params, Symbol(param.name))
+    return
+end
+
+MOI.supports(::Optimizer, ::MOI.Silent) = true
+function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
+    optimizer.silent = value
+    return
+end
+
+MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent
+
+MOI.get(::Optimizer, ::MOI.SolverName) = "SDPLR"
+
+function MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals})
+    return false
+end
+
+function MOI.supports_add_constrained_variables(
+    ::Optimizer,
+    ::Type{<:SupportedSets},
+)
+    return true
+end
+
+function _new_block(model::Optimizer, set::MOI.Nonnegatives)
+    push!(model.blksz, MOI.dimension(set))
+    push!(model.blktype, Cchar('d'))
+    blk = length(model.blksz)
+    for i in 1:MOI.dimension(set)
+        push!(model.varmap, (blk, i, i))
+    end
+    return
+end
+
+function _new_block(model::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
+    push!(model.blksz, set.side_dimension)
+    push!(model.blktype, Cchar('s'))
+    blk = length(model.blksz)
+    for j in 1:set.side_dimension
+        for i in 1:j
+            push!(model.varmap, (blk, i, j))
+        end
+    end
+    return
+end
+
+function MOI.add_constrained_variables(model::Optimizer, set::SupportedSets)
+    offset = length(model.varmap)
+    _new_block(model, set)
+    ci = MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(offset + 1)
+    for i in eachindex(model.Ainfo_entptr)
+        _fill_until(
+            model,
+            length(model.blksz),
+            model.Ainfo_entptr[i],
+            model.Ainfo_type[i],
+            _prev(model, i),
+        )
+    end
+    _fill_until(
+        model,
+        length(model.blksz),
+        model.Cinfo_entptr,
+        model.Cinfo_type,
+        length(model.Cent),
+    )
+    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
+end
+
+function MOI.supports_constraint(
+    ::Optimizer,
+    ::Type{MOI.ScalarAffineFunction{Cdouble}},
+    ::Type{MOI.EqualTo{Cdouble}},
+)
+    return true
+end
+
+function _isless(t1::MOI.VectorAffineTerm, t2::MOI.VectorAffineTerm)
+    if t1.scalar_term.variable.value == t2.scalar_term.variable.value
+        return isless(t1.output_index, t2.output_index)
+    else
+        return isless(
+            t1.scalar_term.variable.value,
+            t2.scalar_term.variable.value,
+        )
+    end
+end
+
+function _prev(model::Optimizer, i)
+    prev = 0
+    for j in i:-1:1
+        if !isempty(model.Ainfo_entptr[j])
+            prev = last(model.Ainfo_entptr[j])
+        end
+    end
+    return prev
+end
+
+function _fill_until(
+    model::Optimizer,
+    numblk,
+    entptr::Vector{Csize_t},
+    type::Vector{Cchar},
+    prev,
+)
+    @assert length(type) == length(entptr)
+    while length(type) < numblk
+        blk = length(type) + 1
+        push!(type, model.blktype[blk])
+        push!(entptr, prev)
+    end
+    @assert length(type) == numblk
+    @assert length(entptr) == numblk
+    return
+end
+
+function _fill!(
+    model,
+    ent,
+    row,
+    col,
+    entptr::Vector{Csize_t},
+    type::Vector{Cchar},
+    func,
+)
+    for t in MOI.Utilities.canonical(func).terms
+        blk, i, j = model.varmap[t.variable.value]
+        _fill_until(model, blk, entptr, type, length(ent))
+        coef = t.coefficient
+        if i != j
+            coef /= 2
+        end
+        push!(ent, coef)
+        push!(row, i)
+        push!(col, j)
+    end
+    _fill_until(model, length(model.blksz), entptr, type, length(ent))
+    @assert length(entptr) == length(model.blksz)
+    @assert length(type) == length(model.blksz)
+end
+
+function MOI.add_constraint(
+    model::Optimizer,
+    func::MOI.ScalarAffineFunction{Cdouble},
+    set::MOI.EqualTo{Cdouble},
+)
+    push!(model.Ainfo_entptr, Csize_t[])
+    push!(model.Ainfo_type, Cchar[])
+    _fill!(
+        model,
+        model.Aent,
+        model.Arow,
+        model.Acol,
+        model.Ainfo_entptr[end],
+        model.Ainfo_type[end],
+        func,
+    )
+    push!(model.b, MOI.constant(set) - MOI.constant(func))
+    return MOI.ConstraintIndex{typeof(func),typeof(set)}(length(model.b))
+end
+
+MOI.supports_incremental_interface(::Optimizer) = true
+
+function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
+    return MOI.Utilities.default_copy_to(dest, src)
+end
+
+function MOI.optimize!(model::Optimizer)
+    CAent = vcat(model.Cent, model.Aent)
+    CArow = vcat(model.Crow, model.Arow)
+    CAcol = vcat(model.Ccol, model.Acol)
+    @assert length(model.Cinfo_entptr) == length(model.blksz)
+    @assert length(model.Cinfo_type) == length(model.blksz)
+    for i in eachindex(model.Ainfo_entptr)
+        @assert length(model.Ainfo_entptr[i]) == length(model.blksz)
+        @assert length(model.Ainfo_type[i]) == length(model.blksz)
+    end
+    CAinfo_entptr = reduce(
+        vcat,
+        vcat([model.Cinfo_entptr], model.Ainfo_entptr),
+        init = Csize_t[],
+    )
+    for i in (length(model.Cinfo_entptr)+1):length(CAinfo_entptr)
+        CAinfo_entptr[i] += length(model.Cent)
+    end
+    push!(CAinfo_entptr, length(CAent))
+    CAinfo_type =
+        reduce(vcat, vcat([model.Cinfo_type], model.Ainfo_type), init = Cchar[])
+    maxranks = default_maxranks(
+        model.blktype,
+        model.blksz,
+        CAinfo_entptr,
+        length(model.b),
+    )
+    Rsizes = map(eachindex(model.blktype)) do k
+        if model.blktype[k] == Cchar('s')
+            return model.blksz[k] * maxranks[k]
+        else
+            @assert model.blktype[k] == Cchar('d')
+            return model.blksz[k]
+        end
+    end
+    model.Rmap = [0; cumsum(Rsizes)]
+    # In `main.c`, it does `(rand() / RAND_MAX) - (rand() - RAND_MAX)`` to take the difference between
+    # two numbers between 0 and 1. Here, Julia's `rand()`` is already between 0 and 1 so we don't have
+    # to divide by anything.
+    nr = last(model.Rmap)
+    R = rand(nr) - rand(nr)
+    _, model.R, model.lambda, model.ranks, model.pieces = solve(
+        model.blksz,
+        model.blktype,
+        model.b,
+        CAent,
+        CArow,
+        CAcol,
+        CAinfo_entptr,
+        CAinfo_type,
+        params = model.params,
+        maxranks = maxranks,
+        R = R,
+    )
+    return
+end
+
+function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
+    majiter, iter, λupdate, CG, curr_CG, totaltime, σ, overallsc =
+        optimizer.pieces
+    return "majiter = $majiter, iter = $iter, λupdate = $λupdate, CG = $CG, curr_CG = $curr_CG, totaltime = $totaltime, σ = $σ, overallsc = $overallsc"
+end
+function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec)
+    return optimizer.pieces[6]
+end
+
+function MOI.is_empty(optimizer::Optimizer)
+    return isempty(optimizer.b) && isempty(optimizer.varmap)
+end
+
+function MOI.empty!(optimizer::Optimizer)
+    optimizer.objective_constant = 0.0
+    optimizer.objective_sign = 1
+    empty!(optimizer.blksz)
+    empty!(optimizer.blktype)
+    empty!(optimizer.varmap)
+    empty!(optimizer.b)
+    empty!(optimizer.Cent)
+    empty!(optimizer.Crow)
+    empty!(optimizer.Ccol)
+    empty!(optimizer.Cinfo_entptr)
+    empty!(optimizer.Cinfo_type)
+    empty!(optimizer.Aent)
+    empty!(optimizer.Arow)
+    empty!(optimizer.Acol)
+    empty!(optimizer.Ainfo_entptr)
+    empty!(optimizer.Ainfo_type)
+    empty!(optimizer.Rmap)
+    optimizer.R = nothing
+    empty!(optimizer.lambda)
+    empty!(optimizer.ranks)
+    optimizer.pieces = nothing
+    return
+end
+
+function MOI.supports(
+    ::Optimizer,
+    ::Union{
+        MOI.ObjectiveSense,
+        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Cdouble}},
+    },
+)
+    return true
+end
+
+function MOI.set(
+    model::Optimizer,
+    ::MOI.ObjectiveSense,
+    sense::MOI.OptimizationSense,
+)
+    sign = sense == MOI.MAX_SENSE ? -1 : 1
+    if model.objective_sign != sign
+        rmul!(model.b, -1)
+        model.objective_sign = sign
+    end
+    return
+end
+
+function MOI.set(
+    model::Optimizer,
+    ::MOI.ObjectiveFunction,
+    func::MOI.ScalarAffineFunction{Cdouble},
+)
+    model.objective_constant = MOI.constant(func)
+    empty!(model.Cent)
+    empty!(model.Crow)
+    empty!(model.Ccol)
+    empty!(model.Cinfo_entptr)
+    empty!(model.Cinfo_type)
+    _fill!(
+        model,
+        model.Cent,
+        model.Crow,
+        model.Ccol,
+        model.Cinfo_entptr,
+        model.Cinfo_type,
+        func,
+    )
+    return
+end
+
+function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
+    if isnothing(model.pieces)
+        return MOI.OPTIMIZE_NOT_CALLED
+    else
+        return MOI.LOCALLY_SOLVED
+    end
+end
+
+function MOI.get(m::Optimizer, attr::MOI.PrimalStatus)
+    if attr.result_index > MOI.get(m, MOI.ResultCount())
+        return MOI.NO_SOLUTION
+    end
+    return MOI.FEASIBLE_POINT
+end
+
+function MOI.get(m::Optimizer, attr::MOI.DualStatus)
+    if attr.result_index > MOI.get(m, MOI.ResultCount())
+        return MOI.NO_SOLUTION
+    end
+    return MOI.FEASIBLE_POINT
+end
+
+function MOI.get(model::Optimizer, ::MOI.ResultCount)
+    if MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
+        return 0
+    else
+        return 1
+    end
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.VariablePrimal,
+    vi::MOI.VariableIndex,
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    blk, i, j = optimizer.varmap[vi.value]
+    I = (optimizer.Rmap[blk]+1):optimizer.Rmap[blk+1]
+    if optimizer.blktype[blk] == Cchar('s')
+        d = optimizer.blksz[blk]
+        U = reshape(optimizer.R[I], d, div(length(I), d))
+        return U[i, :]' * U[j, :]
+    else
+        @assert optimizer.blktype[blk] == Cchar('d')
+        return optimizer.R[I[i]]
+    end
+end
+
+function MOI.get(
+    optimizer::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex{
+        MOI.ScalarAffineFunction{Cdouble},
+        MOI.EqualTo{Cdouble},
+    },
+)
+    MOI.check_result_index_bounds(optimizer, attr)
+    return optimizer.lambda[ci.value]
+end
diff --git a/src/SDPLR.jl b/src/SDPLR.jl
index c78325f..e76ea5f 100644
--- a/src/SDPLR.jl
+++ b/src/SDPLR.jl
@@ -25,30 +25,26 @@ Base.@kwdef struct Parameters
 end
 
 # See `macros.h`
-datablockind(data, block, numblock) = ((data + 1) - 1) * numblock + block
+datablockind(data, block, numblock) = data * numblock + block
 
-import Random
 function default_R(blktype::Vector{Cchar}, blksz, maxranks)
     # See `getstorage` in `main.c`
     nr = sum(eachindex(blktype)) do k
         if blktype[k] == Cchar('s')
             return blksz[k] * maxranks[k]
-        elseif blktype[k] == Cchar('d')
-            return blksz[k]
         else
-            return 0
+            @assert blktype[k] == Cchar('d')
+            return blksz[k]
         end
     end
-    # In `main.c`, it does (rand() / RAND_MAX) - (rand() - RAND_MAX) to take the difference between
-    # two numbers between 0 and 1. Here, Julia's rand() is already between 0 and 1 so we don't have
-    # to divide.
-    Random.seed!(925)
+    # In `main.c`, it does `(rand() / RAND_MAX) - (rand() - RAND_MAX)`` to take the difference between
+    # two numbers between 0 and 1. Here, Julia's `rand()`` is already between 0 and 1 so we don't have
+    # to divide by anything.
     return rand(nr) - rand(nr)
 end
 
-function default_maxranks(blktype, blksz, CAinfo_entptr)
+function default_maxranks(blktype, blksz, CAinfo_entptr, m)
     numblk = length(blktype)
-    m = div(length(CAinfo_entptr) - 1, numblk) - 1
     # See `getstorage` in `main.c`
     return map(eachindex(blktype)) do k
         if blktype[k] == Cchar('s')
@@ -57,14 +53,27 @@ function default_maxranks(blktype, blksz, CAinfo_entptr)
                 return CAinfo_entptr[ind+1] > CAinfo_entptr[ind]
             end
             return Csize_t(min(isqrt(2cons) + 1, blksz[k]))
-        elseif blktype[k] == Cchar('d')
-            return Csize_t(1)
         else
-            return Csize_t(0)
+            @assert blktype[k] == Cchar('d')
+            return Csize_t(1)
         end
     end
 end
 
+"""
+SDPA format (see `MOI.FileFormats.SDPA.Model`) with
+matrices `C`, `A_1`, ..., `A_m`, `X` that are block
+diagonal with `numblk` blocks and `b` is a length-`m`
+vector.
+
+Each block `1 <= k <= numblk` is has dimension `blksz[k] × blksz[k]`.
+The `k`th block of `X` is computed as `R * R'` where `R` is of size
+`blksz[k] × maxranks[k]` if `blktype[k]` is `Cchar('s')` and
+`Diagonal(R)` where `R` is a vector of size `blksz[k]` if `blktype[k]`
+is `Cchar('d')`.
+
+The `CA...` arguments specify the `C` and `A_i` matrices.
+"""
 function solve(
     blksz::Vector{Cptrdiff_t},
     blktype::Vector{Cchar},
@@ -75,7 +84,12 @@ function solve(
     CAinfo_entptr::Vector{Csize_t},
     CAinfo_type::Vector{Cchar};
     params::Parameters = Parameters(),
-    maxranks::Vector{Csize_t} = default_maxranks(blktype, blksz, CAinfo_entptr),
+    maxranks::Vector{Csize_t} = default_maxranks(
+        blktype,
+        blksz,
+        CAinfo_entptr,
+        length(b),
+    ),
     ranks::Vector{Csize_t} = copy(maxranks),
     R::Vector{Cdouble} = default_R(blktype, blksz, maxranks),
     lambda::Vector{Cdouble} = zeros(length(b)),
@@ -91,37 +105,64 @@ function solve(
     @assert length(maxranks) == numblk
     @assert length(ranks) == numblk
     @assert length(pieces) == 8
-    ret = @ccall SDPLR.SDPLR_jll.libsdplr.sdplrlib(
-        m::Csize_t,
-        numblk::Csize_t,
-        blksz::Ptr{Cptrdiff_t},
-        blktype::Ptr{Cchar},
-        b::Ptr{Cdouble},
-        CAent::Ptr{Cdouble},
-        CArow::Ptr{Csize_t},
-        CAcol::Ptr{Csize_t},
-        CAinfo_entptr::Ptr{Csize_t},
-        CAinfo_type::Ptr{Cchar},
-        params.numbfgsvecs::Csize_t,
-        params.rho_f::Cdouble,
-        params.rho_c::Cdouble,
-        params.sigmafac::Cdouble,
-        params.rankreduce::Csize_t,
-        params.gaptol::Cdouble,
-        params.checkbd::Cptrdiff_t,
-        params.typebd::Csize_t,
-        params.dthresh_dim::Csize_t,
-        params.dthresh_dens::Cdouble,
-        params.timelim::Csize_t,
-        params.rankredtol::Cdouble,
-        params.printlevel::Csize_t,
-        R::Ptr{Cdouble},
-        lambda::Ptr{Cdouble},
-        maxranks::Ptr{Csize_t},
-        ranks::Ptr{Csize_t},
-        pieces::Ptr{Cdouble},
-    )::Csize_t
-    return ret, R, lambda, ranks
+    @assert CAinfo_entptr[1] == 0
+    @assert CAinfo_entptr[end] == length(CArow)
+    k = 0
+    for _ in eachindex(b)
+        for blk in eachindex(blksz)
+            k += 1
+            @assert CAinfo_entptr[k] <= CAinfo_entptr[k+1]
+            for j in ((CAinfo_entptr[k]+1):CAinfo_entptr[k+1])
+                @assert blktype[blk] == CAinfo_type[k]
+                @assert 1 <= CArow[j] <= blksz[blk]
+                @assert 1 <= CAcol[j] <= blksz[blk]
+                if CAinfo_type[k] == Cchar('s')
+                    @assert CArow[j] <= CAcol[j]
+                else
+                    @assert CAinfo_type[k] == Cchar('d')
+                    @assert CArow[j] == CAcol[j]
+                end
+            end
+        end
+    end
+    GC.@preserve blksz blktype b CAent CArow CAcol CAinfo_entptr CAinfo_type R lambda maxranks ranks pieces begin
+        ret = @ccall SDPLR.SDPLR_jll.libsdplr.sdplrlib(
+            m::Csize_t,
+            numblk::Csize_t,
+            blksz::Ptr{Cptrdiff_t},
+            blktype::Ptr{Cchar},
+            b::Ptr{Cdouble},
+            CAent::Ptr{Cdouble},
+            CArow::Ptr{Csize_t},
+            CAcol::Ptr{Csize_t},
+            CAinfo_entptr::Ptr{Csize_t},
+            CAinfo_type::Ptr{Cchar},
+            params.numbfgsvecs::Csize_t,
+            params.rho_f::Cdouble,
+            params.rho_c::Cdouble,
+            params.sigmafac::Cdouble,
+            params.rankreduce::Csize_t,
+            params.gaptol::Cdouble,
+            params.checkbd::Cptrdiff_t,
+            params.typebd::Csize_t,
+            params.dthresh_dim::Csize_t,
+            params.dthresh_dens::Cdouble,
+            params.timelim::Csize_t,
+            params.rankredtol::Cdouble,
+            params.printlevel::Csize_t,
+            # We can see in `source/main.c` that `R - 1` and `lambda - 1`
+            # are passed to `sdplrlib` so we also need to shift by `-1`
+            # by using `pointer(_, 0)`.
+            pointer(R, 0)::Ptr{Cdouble},
+            pointer(lambda, 0)::Ptr{Cdouble},
+            maxranks::Ptr{Csize_t},
+            ranks::Ptr{Csize_t},
+            pieces::Ptr{Cdouble},
+        )::Csize_t
+    end
+    return ret, R, lambda, ranks, pieces
 end
 
+include("MOI_wrapper.jl")
+
 end # module
diff --git a/test/Project.toml b/test/Project.toml
index 5024722..7adae17 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,3 +1,5 @@
 [deps]
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SDPLR = "56161740-ea4e-4253-9d15-43c62ff94d95"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index b0eef8e..1a3ab98 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1 +1,2 @@
 include("test_vibra.jl")
+include("test_simple.jl")
diff --git a/test/simple.jl b/test/simple.jl
new file mode 100644
index 0000000..9911ef5
--- /dev/null
+++ b/test/simple.jl
@@ -0,0 +1,8 @@
+blksz = Cptrdiff_t[2]
+blktype = Cchar['s']
+b = Cdouble[1]
+CAinfo_entptr = Csize_t[0, 2, 3]
+CAent = Cdouble[1, 1, 0.5]
+CArow = Csize_t[1, 2, 1]
+CAcol = Csize_t[1, 2, 2]
+CAinfo_type = Cchar['s', 's']
diff --git a/test/test_simple.jl b/test/test_simple.jl
new file mode 100644
index 0000000..01776c6
--- /dev/null
+++ b/test/test_simple.jl
@@ -0,0 +1,49 @@
+using Test
+import SDPLR
+import Random
+import MathOptInterface as MOI
+
+# This is `test_conic_PositiveSemidefiniteConeTriangle_VectorOfVariables`
+@testset "Solve simple with sdplrlib" begin
+    include("simple.jl")
+    # The `925` seed is taken from SDPLR's `main.c`
+    Random.seed!(925)
+    ret, R, lambda, ranks, pieces = SDPLR.solve(
+        blksz,
+        blktype,
+        b,
+        CAent,
+        CArow,
+        CAcol,
+        CAinfo_entptr,
+        CAinfo_type,
+    )
+    @test iszero(ret)
+    @test length(R) == 4
+    U = reshape(R, 2, 2)
+    @test U * U' ≈ ones(2, 2) rtol = 1e-2
+    @test lambda ≈ [2.0] atol = 1e-2
+    @test ranks == Csize_t[2]
+end
+
+@testset "MOI wrapper" begin
+    atol = rtol = 1e-2
+    model = SDPLR.Optimizer()
+    X, cX = MOI.add_constrained_variables(
+        model,
+        MOI.PositiveSemidefiniteConeTriangle(2),
+    )
+    c = MOI.add_constraint(model, 1.0 * X[2], MOI.EqualTo(1.0))
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    obj = 1.0 * X[1] + 1.0 * X[3]
+    MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED
+    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
+    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
+    Xv = ones(Float64, 3)
+    @test MOI.get(model, MOI.VariablePrimal(), X) ≈ Xv atol = atol rtol = rtol
+    attr = MOI.ConstraintDual()
+    @test MOI.get(model, attr, c) ≈ 2 atol = atol rtol = rtol
+end
diff --git a/test/test_vibra.jl b/test/test_vibra.jl
index 48d3cea..2d7a0e4 100644
--- a/test/test_vibra.jl
+++ b/test/test_vibra.jl
@@ -1,11 +1,14 @@
 using Test
 import SDPLR
+import Random
 @testset "Solve vibra with sdplr executable" begin
     SDPLR.solve_sdpa_file("vibra1.dat-s")
 end
 @testset "Solve vibra with sdplrlib" begin
     include("vibra.jl")
-    ret, R, lambda, ranks = SDPLR.solve(
+    # The `925` seed is taken from SDPLR's `main.c`
+    Random.seed!(925)
+    ret, R, lambda, ranks, pieces = SDPLR.solve(
         blksz,
         blktype,
         b,
@@ -17,6 +20,6 @@ end
     )
     @test iszero(ret)
     @test length(R) == 477
-    @test sum(lambda) ≈ -40.8133 rtol = 1e-3
+    @test sum(lambda) ≈ -40.8133 rtol = 1e-2
     @test ranks == Csize_t[9, 9, 1]
 end
diff --git a/test/vibra.jl b/test/vibra.jl
index 7753eaa..eaa93b5 100644
--- a/test/vibra.jl
+++ b/test/vibra.jl
@@ -1,5 +1,4 @@
 m = 36
-numblk = 3
 blksz = Cptrdiff_t[24, 25, 36]
 blktype = Cchar['s', 's', 'd']
 b = ones(m)