diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index d8b81ad..6b6c552 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -88,7 +88,7 @@ jobs:
           cd $SOURCES_DIR
           ./configure --prefix=$PETSC_INSTALL --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 \
              --download-mumps --download-scalapack --download-parmetis --download-metis \
-             --download-ptscotch --with-debugging --with-x=0 --with-shared=1 \
+             --download-fblaslapack --download-ptscotch --with-debugging --with-x=0 --with-shared=1 \
              --with-mpi=1 --with-64-bit-indices
           make
           make install
diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl
index ff806b3..849a059 100644
--- a/src/LinearSolvers/LinearSolvers.jl
+++ b/src/LinearSolvers/LinearSolvers.jl
@@ -28,6 +28,7 @@ export SchurComplementSolver
 export SchwarzLinearSolver
 
 export CallbackSolver
+export NullspaceSolver
 
 # Wrappers for IterativeSolvers.jl
 export IS_ConjugateGradientSolver
@@ -66,5 +67,6 @@ include("MatrixSolvers.jl")
 include("SchwarzLinearSolvers.jl")
 
 include("CallbackSolver.jl")
+include("NullspaceSolvers.jl")
 
 end
\ No newline at end of file
diff --git a/src/LinearSolvers/NullspaceSolvers.jl b/src/LinearSolvers/NullspaceSolvers.jl
new file mode 100644
index 0000000..3f8ee01
--- /dev/null
+++ b/src/LinearSolvers/NullspaceSolvers.jl
@@ -0,0 +1,122 @@
+
+"""
+    struct NullspaceSolver{A,B}
+      solver :: A <: LinearSolver
+      kernel :: B <: Union{AbstractMatrix, Vector{<:AbstractVector}}
+      constrain_matrix :: Bool = true
+    end
+
+Solver that computes the solution of a linear system `A⋅x = b` with a kernel constraint,
+i.e the returned solution is orthogonal to the provided kernel.
+
+We assume the kernel is provided as a matrix `K` of dimensions `(k,n)` with `n` the dimension
+of the original system and `k` the number of kernel vectors.
+
+Two modes of operation are supported:
+
+- If `constrain_matrix` is `true`, the solver will explicitly constrain the system matrix. 
+  I.e we consider the augmented system `Â⋅x̂ = b̂` where
+  -  `Ak = [A, K'; K, 0]`
+  -  `x̂ = [x; λ]`
+  -  `b̂ = [b; 0]`
+  This is often the only option for direct solvers, which require the system matrix to be
+  invertible. This is only supported in serial (its performance bottleneck in parallel).
+
+- If `constrain_matrix` is `false`, the solver preserve the original system and simply 
+  project the initial guess and the solution onto the orthogonal complement of the kernel. 
+  This option is more suitable for iterative solvers, which usually do not require the 
+  system matrix to be invertible (e.g. GMRES, BiCGStab, etc).
+"""
+struct NullspaceSolver{A,B}
+  solver    :: A
+  nullspace :: B
+  constrain_matrix :: Bool
+
+  function NullspaceSolver(
+    solver::LinearSolver,
+    nullspace::NullSpace; 
+    constrain_matrix::Bool = true
+  )
+    A, B = typeof(solver), typeof(nullspace)
+    new{A,B}(solver, nullspace, constrain_matrix)
+  end
+end
+
+struct NullspaceSolverSS{A} <: Algebra.SymbolicSetup
+  solver :: A
+end
+
+function Gridap.Algebra.symbolic_setup(solver::NullspaceSolver,A)
+  return NullspaceSolverSS(solver)
+end
+
+struct NullspaceSolverNS{S}
+  solver
+  ns
+  caches
+end
+
+function Algebra.numerical_setup(ss::NullspaceSolverSS, A::AbstractMatrix)
+  solver = ss.solver
+  N = solver.nullspace
+  nK, nV = size(N)
+  @assert size(A,1) == size(A,2) == nV
+  if solver.constrain_matrix
+    K = SolverInterfaces.matrix_representation(N)
+    mat = [A K; K' zeros(nK,nK)] # TODO: Explore reusing storage for A
+    display(mat)
+    display(A)
+    display(K)
+  else
+    mat = A
+  end
+  S = ifelse(solver.constrain_matrix, :constrained, :projected)
+  ns = numerical_setup(symbolic_setup(solver.solver, mat), mat)
+  caches = (allocate_in_domain(mat), allocate_in_domain(mat))
+  return NullspaceSolverNS{S}(solver, ns, caches)
+end
+
+function Algebra.numerical_setup!(ns::NullspaceSolverNS, A::AbstractMatrix)
+  solver = ns.solver
+  N = solver.nullspace
+  nK, nV = size(N)
+  @assert size(A,1) == size(A,2) == nV
+  if solver.constrain_matrix
+    K = SolverInterfaces.matrix_representation(N)
+    mat = [A K; K' zeros(nK,nK)] # TODO: Explore reusing storage for A
+  else
+    mat = A
+  end
+  numerical_setup!(ns.ns, mat)
+  return ns
+end
+
+function Algebra.solve!(x::AbstractVector, ns::NullspaceSolverNS{:constrained}, b::AbstractVector)
+  solver = ns.solver
+  A_ns, caches = ns.ns, ns.caches
+  N = solver.nullspace
+  w1, w2 = caches
+  nK, nV = size(N)
+  @assert length(x) == nV
+  @assert length(b) == nV
+  w1[1:nV] .= x
+  w1[nV+1:nV+nK] .= 0.0
+  w2[1:nV] .= b
+  w2[nV+1:nV+nK] .= 0.0
+  solve!(w1, A_ns, w2)
+  x .= w1[1:nV]
+  return x
+end
+
+function Algebra.solve!(x::AbstractVector, ns::NullspaceSolverNS{:projected}, b::AbstractVector)
+  solver = ns.solver
+  A_ns, caches = ns.ns, ns.caches
+  N = solver.nullspace
+  w1, w2 = caches
+  
+  w1, α = SolverInterfaces.project!(w1, N, x)
+  x .-= w1
+  solve!(x, A_ns, b)
+  x .+= w1
+  return x
+end
diff --git a/src/SolverInterfaces/NullSpaces.jl b/src/SolverInterfaces/NullSpaces.jl
new file mode 100644
index 0000000..1559411
--- /dev/null
+++ b/src/SolverInterfaces/NullSpaces.jl
@@ -0,0 +1,139 @@
+struct NullSpace{T,VT}
+  V :: VT
+  function NullSpace(
+    V::AbstractVector{<:AbstractVector{T}}
+  ) where T
+    n = length(first(V))
+    @assert all(length(v) == n for v in V)
+    VT = typeof(V)
+    new{T,VT}(V)
+  end
+end
+
+Base.size(N::NullSpace) = (length(N.V),length(first(N.V)))
+Base.size(N::NullSpace, i::Int) = size(N)[i]
+Base.merge(a::NullSpace{T}, b::NullSpace{T}) where T = NullSpace(vcat(a.V, b.V))
+
+function matrix_representation(N::NullSpace)
+  return stack(N.V)
+end
+
+NullSpace(v::AbstractVector{<:Number}) = NullSpace([v])
+
+function NullSpace(A::Matrix)
+  V = eachcol(nullspace(A))
+  return NullSpace(V)
+end
+
+function NullSpace(f::Function,X::FESpace,Λ::FESpace)
+  A = assemble_matrix(f, X, Λ)
+  return NullSpace(A)
+end
+
+function is_orthonormal(N::NullSpace; tol = 1.e-12)
+  for w in N.V
+    !(abs(norm(w) - 1.0) < tol) && return false
+  end
+  return is_orthogonal(N, tol = tol)
+end
+
+function is_orthogonal(N::NullSpace; tol = 1.e-12)
+  for (k,w) in enumerate(N.V)
+    for v in N.V[k+1:end]
+      !(abs(dot(w, v)) < tol) && return false
+    end
+  end
+  return true
+end
+
+function is_orthogonal(N::NullSpace, v::AbstractVector; tol = 1.e-12)
+  @assert length(v) == size(N,2)
+  for w in N.V
+    !(abs(dot(w, v)) < tol) && return false
+  end
+  return true
+end
+
+function is_orthogonal(N::NullSpace, A::AbstractMatrix; tol = 1.e-12)
+  @assert length(v) == size(N,2)
+  v = allocate_in_range(A)
+  for w in N.V
+    mul!(v, A, w)
+    !(abs(norm(v)) < tol) && return false
+  end
+  return true
+end
+
+function make_orthonormal!(N::NullSpace; method = :gram_schmidt)
+  if method == :gram_schmidt
+    gram_schmidt!(N.V)
+  elseif method == :modified_gram_schmidt
+    modified_gram_schmidt!(N.V)
+  else
+    error("Unknown method: $method")
+  end
+  return N
+end
+
+function gram_schmidt!(V::AbstractVector{<:AbstractVector{T}}) where T
+  n = length(V)
+  for j in 1:n
+    for i in 1:j-1
+      α = dot(V[j], V[i])
+      V[j] .-= α .* V[i]
+    end
+    V[j] ./= norm(V[j])
+  end
+  return V
+end
+
+function modified_gram_schmidt!(V::AbstractVector{<:AbstractVector{T}}) where T
+  n = length(V)
+  for j in 1:n
+    V[j] ./= norm(V[j])
+    for i in j+1:n
+      α = dot(V[j], V[i])
+      V[i] .-= α .* V[j]
+    end
+  end
+  return V
+end
+
+function project(N::NullSpace,v)
+  p = similar(v)
+  project!(p,N,v)
+end
+
+function project!(p,N::NullSpace{T},v) where T
+  @assert length(v) == size(N,2)
+  α = zeros(T, size(N,1))
+  fill!(p,0.0)
+  for (k,w) in enumerate(N.V)
+    α[k] = dot(v, w)
+    p .+= α[k] .* w
+  end
+  return p, α
+end
+
+function make_orthogonal!(N::NullSpace{T},v) where T
+  @assert length(v) == size(N,2)
+  α = zeros(T, size(N,1))
+  for (k,w) in enumerate(N.V)
+    α[k] = dot(v, w)
+    v .-= α[k] .* w
+  end
+  return v, α
+end
+
+function reconstruct(N::NullSpace,v,α)
+  w = copy(v)
+  reconstruct!(N,w,α)
+  return w
+end
+
+function reconstruct!(N::NullSpace,v,α)
+  for (k,w) in enumerate(N.V)
+    v .+= α[k] .* w
+  end
+  return v
+end
diff --git a/src/SolverInterfaces/SolverInterfaces.jl b/src/SolverInterfaces/SolverInterfaces.jl
index 4bcb3d3..9da52e0 100644
--- a/src/SolverInterfaces/SolverInterfaces.jl
+++ b/src/SolverInterfaces/SolverInterfaces.jl
@@ -1,5 +1,7 @@
 module SolverInterfaces
 
+using LinearAlgebra
+
 using Gridap
 using Gridap.Helpers
 using Gridap.Algebra
@@ -15,6 +17,7 @@ include("GridapExtras.jl")
 include("SolverTolerances.jl")
 include("ConvergenceLogs.jl")
 include("SolverInfos.jl")
+include("NullSpaces.jl")
 
 export SolverVerboseLevel, SolverConvergenceFlag
 export SolverTolerances, get_solver_tolerances, set_solver_tolerances!
@@ -22,4 +25,6 @@ export ConvergenceLog, init!, update!, finalize!, reset!, print_message, set_dep
 
 export SolverInfo
 
+export NullSpace
+
 end
\ No newline at end of file
diff --git a/test/_dev/nullspaces.jl b/test/_dev/nullspaces.jl
new file mode 100644
index 0000000..f980f46
--- /dev/null
+++ b/test/_dev/nullspaces.jl
@@ -0,0 +1,66 @@
+
+using Gridap
+using Gridap.Algebra
+
+using GridapSolvers
+using GridapSolvers.SolverInterfaces, GridapSolvers.LinearSolvers
+
+using GridapSolvers.SolverInterfaces: NullSpace, is_orthonormal, is_orthogonal, gram_schmidt!, modified_gram_schmidt!
+using GridapSolvers.SolverInterfaces: project, project!, make_orthogonal!, reconstruct, reconstruct!
+
+A = [1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
+N = NullSpace(A)
+
+is_orthonormal(N)
+is_orthogonal(N,[1.0,0.0,0.0])
+is_orthogonal(N,[0.0,1.0,0.0])
+
+v = [1.0,2.0,3.0]
+p, α = project(N,v)
+w, β = make_orthogonal!(N,copy(v))
+u = reconstruct(N,w,α)
+u ≈ v
+
+V1 = [[2.0,1.0,1.0],[1.0,2.0,1.0],[1.0,1.0,1.0]]
+gram_schmidt!(V1)
+is_orthonormal(NullSpace(V1))
+
+V2 = [[2.0,1.0,1.0],[1.0,2.0,1.0],[1.0,1.0,1.0]]
+modified_gram_schmidt!(V2)
+is_orthonormal(NullSpace(V2))
+V1 ≈ V2
+
+
+############################################################################################
+
+model = CartesianDiscreteModel((0,1,0,1),(4,4))
+Ω = Triangulation(model)
+Γ = BoundaryTriangulation(model)
+
+dΩ = Measure(Ω,4)
+dΓ = Measure(Γ,4)
+n = get_normal_vector(Γ)
+
+V = FESpace(model, ReferenceFE(lagrangian,Float64,1))
+
+u_exact(x) = x[1]*x[2]
+f(x) = Δ(u_exact)(x)
+a(u,v) = ∫(∇(u)⋅∇(v))*dΩ
+l(v) = ∫(f*v)*dΩ + ∫(v*(∇(u_exact)⋅n))*dΓ
+
+op = AffineFEOperator(a,l,V,V)
+A, b = get_matrix(op), get_vector(op)
+
+N = NullSpace(ones(size(A,2)))
+K = SolverInterfaces.matrix_representation(N)
+norm(A*K)
+
+s = GMRESSolver(10;verbose=true)
+s = LUSolver()
+solver = NullspaceSolver(s,N;constrain_matrix=true)
+ns = numerical_setup(symbolic_setup(solver,A),A)
+
+x = allocate_in_domain(A)
+fill!(x,0.0)
+solve!(x,ns,b)
+norm(A*x-b)