diff --git a/docs/src/api.md b/docs/src/api.md
index a43d10e..8362e2d 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -46,8 +46,12 @@ qrm_spposv
 qrm_spposv!
 qrm_least_squares
 qrm_least_squares!
+qrm_least_squares_semi_normal
+qrm_least_squares_semi_normal!
 qrm_min_norm
 qrm_min_norm!
+qrm_min_norm_semi_normal
+qrm_min_norm_semi_normal!
 qrm_refine
 qrm_refine!
 ```
diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index f9001a1..18e4905 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -50,7 +50,9 @@ export qrm_spmat, qrm_spfct,
     qrm_spbackslash!, qrm_spbackslash, \,
     qrm_spposv!, qrm_spposv,
     qrm_least_squares!, qrm_least_squares,
+    qrm_least_squares_semi_normal!, qrm_least_squares_semi_normal,
     qrm_min_norm!, qrm_min_norm,
+    qrm_min_norm_semi_normal!, qrm_min_norm_semi_normal,
     qrm_residual_norm!, qrm_residual_norm,
     qrm_residual_orth!, qrm_residual_orth,
     qrm_refine!, qrm_refine, qrm_set, qrm_get
@@ -348,6 +350,48 @@ function qrm_least_squares! end
 """
 function qrm_least_squares end
 
+@doc raw"""
+qrm_least_squares_semi_normal!(spmat, b, x, z, Δx, y; transp='n')
+
+This function can be used to solve a linear least squares problem
+
+```math
+\min \|Ax − b\|_2
+```
+
+in the case where the input matrix is square or overdetermined.
+Contrary to `qrm_least_squares!`, this function allows to solve the problem without storing the Q-factor in the QR factorization of A.
+
+It is a shortcut for the sequence
+
+    qrm_analyse!(spmat, spfct, transp  = 'n')
+    qrm_factorize!(spmat, spfct, transp = 'n')
+    qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = 't')
+    qrm_solve!(spfct, z, y, transp = 't')
+    qrm_solve!(spfct, y, x, transp = 'n')
+
+    qrm_refine!(spmat, spfct, x, z, Δx, y)
+
+Remark that the Q-factor is not used in this sequence but rather A and R. 
+
+#### Input Arguments :
+
+* `spmat`: the input matrix.
+* `spfct`: a sparse factorization object of type `qrm_spfct`.
+* `b`: the right-hand side(s).
+* `x`: the solution vector(s).
+* `Δx`: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x.
+* `z`: an auxiliary vector (or matrix if b and x are matrices) used to store the value z = Aᵀb, the size of this vector (resp. matrix) is the same as x and can be unitialized when the function is called.
+* `y`: an auxiliary vector (or matrix if b and x are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b.
+* `transp`: whether to use A, Aᵀ or Aᴴ. Can be either `'t'`, `'c'` or `'n'`.
+"""
+function qrm_least_squares_semi_normal! end
+
+@doc raw"""
+    x = qrm_least_squares_semi_normal(spmat, b)
+"""
+function qrm_least_squares_semi_normal end
+
 @doc raw"""
     qrm_min_norm!(spmat, b, x; transp='n')
 
@@ -379,6 +423,43 @@ function qrm_min_norm! end
 """
 function qrm_min_norm end
 
+@doc raw"""
+    qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y; transp='n')
+
+This function can be used to solve a linear minimum norm problem
+
+```math
+\min \|x\|_2 \quad s.t. \quad Ax = b
+```
+in the case where the input matrix is square or underdetermined.
+Contrary to `qrm_min_norm!`, this function allows to solve the problem without storing the Q-factor in the QR factorization of Aᵀ.
+It is a shortcut for the sequence
+
+    qrm_analyse!(spmat, spfct, transp = 't')
+    qrm_factorize!(spmat, spfct, transp = 't')
+    qrm_solve!(spfct, b, Δx, transp = 't')
+    qrm_solve!(spfct, Δx, y, transp = 'n')
+    qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = 't')
+
+Remark that the Q-factor is not used in this sequence but rather A and R. 
+
+#### Input Arguments :
+
+* `spmat`: the input matrix.
+* `spfct`: a sparse factorization object of type `qrm_spfct`.
+* `b`: the right-hand side(s).
+* `x`: the solution vector(s).
+* `Δx`: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as x.
+* `y`: an auxiliary vector (or matrix if x and b are matrices) used to compute the solution, the size of this vector (resp. matrix) is the same as b.
+* `transp`: whether to use A, Aᵀ or Aᴴ. Can be either `'t'`, `'c'` or `'n'`.
+"""
+function qrm_min_norm_semi_normal! end
+
+@doc raw"""
+    x = qrm_min_norm_semi_normal(spmat, b)
+"""
+function qrm_min_norm_semi_normal end
+
 @doc raw"""
     qrm_residual_norm!(spmat, b, x, nrm; transp='n')
 
@@ -427,10 +508,10 @@ Given an approximate solution x of the linear system RᵀRx ≈ z where R is the
 
 * `spmat`: the input matrix.
 * `spfct`: a sparse factorization object of type `qrm_spfct`.
-* `x`: the approximate solution vector, the size of this vector is n.
-* `z`: the RHS vector of the linear system, the size of this vector is n.
-* `Δx`: an auxiliary vector used to compute the refinement, the size of this vector is n.
-* `y`: an auxiliary vector used to compute the refinement, the size of this vector is m.
+* `x`: the approximate solution vector(s), the size of this vector is n (or n×k if there are multiple solutions).
+* `z`: the RHS vector(s) of the linear system, the size of this vector is n (or n×k if there are multiple RHS).
+* `Δx`: an auxiliary vector (or matrix if x and z are matrices) used to compute the refinement, the size of this vector (resp. matrix) is n (resp. n×k).
+* `y`: an auxiliary vector (or matrix if x and z are matrices) used to compute the refinement, the size of this vector (resp. matrix) is m (resp. m×k).
 """
 function qrm_refine! end
 
diff --git a/src/utils.jl b/src/utils.jl
index 7f37c95..40b9d86 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -7,19 +7,185 @@ function qrm_refine(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractV
   return x_refined
 end
 
+function qrm_refine(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}) where T
+  Δx = similar(x)
+  y = similar(x, (spfct.fct.m, size(x,2)))
+  x_refined = copy(x)
+
+  qrm_refine!(spmat, spfct, x_refined, z, Δx, y)
+  return x_refined
+end
+
 function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}) where T
   @assert length(x) == spfct.fct.n
   @assert length(z) == spfct.fct.n
   @assert length(Δx) == spfct.fct.n
   @assert length(y) == spfct.fct.m
 
-  transp = T <: Real ? 't' : 'c'
+  t = T <: Real ? 't' : 'c'
+  transp = spfct.fct.n == spmat.mat.n ? 'n' : t # Check whether it is A^T = QR or A = QR.
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
 
-  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = 'n')
-  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = ntransp)
   @. Δx = z - Δx 
 
-  qrm_solve!(spfct, Δx, y, transp = transp)
+  qrm_solve!(spfct, Δx, y, transp = t)
   qrm_solve!(spfct, y, Δx, transp = 'n')
   @. x = x + Δx
-end
\ No newline at end of file
+end
+
+function qrm_refine!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}) where T
+  @assert size(x, 1) == spfct.fct.n
+  @assert size(z, 1) == spfct.fct.n
+  @assert size(Δx, 1) == spfct.fct.n
+  @assert size(y, 1) == spfct.fct.m
+  @assert size(x, 2) == size(z, 2)
+  @assert size(Δx, 2) == size(z, 2)
+  @assert size(y, 2) == size(z, 2)
+
+  t = T <: Real ? 't' : 'c'
+  transp = spfct.fct.n == spmat.mat.n ? 'n' : t # Check whether it is A^T = QR or A = QR.
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+
+  qrm_spmat_mv!(spmat, T(1), x, T(0), y, transp = transp)
+  qrm_spmat_mv!(spmat, T(1), y, T(0), Δx, transp = ntransp)
+  @. Δx = z - Δx 
+
+  qrm_solve!(spfct, Δx, y, transp = t)
+  qrm_solve!(spfct, y, Δx, transp = 'n')
+  @. x = x + Δx
+end
+
+function qrm_min_norm_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}; transp :: Char = 'n') where T
+
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+
+  spfct = qrm_spfct_init(spmat)
+  qrm_set(spfct, "qrm_keeph", 0)
+
+  if typeof(b) <: AbstractVector{T}
+    x = similar(b, n)
+  else 
+    x = similar(b, (n, size(b, 2)))
+  end
+  Δx = similar(x)
+  y = similar(b)
+
+  qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y, transp = transp)
+  return x
+end
+
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}; transp :: Char = 'n') where T
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+  
+  @assert length(x) == n
+  @assert length(b) == m
+  @assert length(Δx) == n
+  @assert length(y) == m
+
+  qrm_analyse!(spmat, spfct, transp = ntransp)
+  qrm_factorize!(spmat, spfct, transp = ntransp)
+  qrm_solve!(spfct, b, Δx, transp = t)
+  qrm_solve!(spfct, Δx, y, transp = 'n')
+  #x = A^T y
+  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = ntransp)
+end
+
+function qrm_min_norm_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}; transp :: Char = 'n') where T
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+  
+  @assert size(x, 1) == n
+  @assert size(b, 1) == m
+  @assert size(Δx, 1) == n
+  @assert size(y, 1) == m
+  @assert size(x, 2) == size(b, 2)
+  @assert size(Δx, 2) == size(b, 2)
+  @assert size(y, 2) == size(b, 2)
+
+  qrm_analyse!(spmat, spfct, transp = ntransp)
+  qrm_factorize!(spmat, spfct, transp = ntransp)
+  qrm_solve!(spfct, b, Δx, transp = t)
+  qrm_solve!(spfct, Δx, y, transp = 'n')
+  #x = A^T y
+  qrm_spmat_mv!(spmat, T(1),  y, T(0), x, transp = ntransp)
+end
+
+function qrm_least_squares_semi_normal(spmat :: qrm_spmat{T}, b :: AbstractVecOrMat{T}; transp :: Char = 'n') where T
+
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+
+  spfct = qrm_spfct_init(spmat)
+  qrm_set(spfct, "qrm_keeph", 0)
+  if typeof(b) <: AbstractVector{T}
+    x = similar(b, n)
+  else 
+    x = similar(b, (n, size(b, 2)))
+  end
+  z = similar(x)
+  Δx = similar(x)
+  y = similar(b)
+
+  qrm_least_squares_semi_normal!(spmat, spfct, b, x, z, Δx, y, transp = transp)
+  return x
+end
+
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractVector{T}, x :: AbstractVector{T}, z :: AbstractVector{T}, Δx :: AbstractVector{T}, y :: AbstractVector{T}; transp :: Char = 'n') where T
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+
+  @assert length(x) == n
+  @assert length(b) == m
+  @assert length(z) == n
+  @assert length(Δx) == n
+  @assert length(y) == m
+
+  qrm_analyse!(spmat, spfct, transp  = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
+  # z = A^T b
+  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = ntransp)
+  #R^T y = z
+  qrm_solve!(spfct, z, y, transp = t)
+  qrm_solve!(spfct, y, x, transp = 'n')
+
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
+end
+
+function qrm_least_squares_semi_normal!(spmat :: qrm_spmat{T}, spfct :: qrm_spfct{T}, b :: AbstractMatrix{T}, x :: AbstractMatrix{T}, z :: AbstractMatrix{T}, Δx :: AbstractMatrix{T}, y :: AbstractMatrix{T}; transp :: Char = 'n') where T
+  
+  n = transp == 'n' ? spmat.mat.n : spmat.mat.m
+  m = transp == 'n' ? spmat.mat.m : spmat.mat.n
+  t = T <: Real ? 't' : 'c'
+  ntransp = transp == 't' || transp == 'c' ? 'n' : t
+
+  @assert size(x, 1) == n
+  @assert size(b, 1) == m
+  @assert size(z, 1) == n
+  @assert size(Δx, 1) == n
+  @assert size(y, 1) == m
+
+  @assert size(x, 2) == size(b, 2)
+  @assert size(z, 2) == size(b, 2)
+  @assert size(Δx, 2) == size(b, 2)
+  @assert size(y, 2) == size(b, 2)
+
+  qrm_analyse!(spmat, spfct, transp  = transp)
+  qrm_factorize!(spmat, spfct, transp = transp)
+  # z = A^T b
+  qrm_spmat_mv!(spmat, T(1),  b, T(0), z, transp = ntransp)
+  qrm_solve!(spfct, z, y, transp = t)
+  qrm_solve!(spfct, y, x, transp = 'n')
+
+  qrm_refine!(spmat, spfct, x, z, Δx, y)
+end
diff --git a/src/wrapper/qr_mumps_api.jl b/src/wrapper/qr_mumps_api.jl
index 2e1f12d..b114917 100644
--- a/src/wrapper/qr_mumps_api.jl
+++ b/src/wrapper/qr_mumps_api.jl
@@ -622,7 +622,9 @@ for (fname, elty) in ((:sqrm_least_squares_c, :Float32   ),
             return nothing
         end
 
-        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n')
+        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_least_squares_semi_normal(spmat, b, transp = transp)
+
             nrhs = 1
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n)
@@ -635,7 +637,9 @@ for (fname, elty) in ((:sqrm_least_squares_c, :Float32   ),
             return x
         end
 
-        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n')
+        function qrm_least_squares(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_least_squares_semi_normal(spmat, b, transp = transp)
+
             nrhs = size(b, 2)
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n, nrhs)
@@ -651,14 +655,14 @@ for (fname, elty) in ((:sqrm_least_squares_c, :Float32   ),
         @inline qrm_least_squares!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='t')
         @inline qrm_least_squares!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='t')
 
-        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_least_squares(spmat.parent, b, transp='t')
-        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_least_squares(spmat.parent, b, transp='t')
+        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='t', seminormal = false)
+        @inline qrm_least_squares(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='t', seminormal = false)
 
         @inline qrm_least_squares!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='c')
         @inline qrm_least_squares!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_least_squares!(spmat.parent, b, x, transp='c')
 
-        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_least_squares(spmat.parent, b, transp='c')
-        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_least_squares(spmat.parent, b, transp='c')
+        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='c', seminormal = false)
+        @inline qrm_least_squares(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_least_squares(spmat.parent, b, transp='c', seminormal = false)
 
     end
 end
@@ -692,7 +696,9 @@ for (fname, elty) in ((:sqrm_min_norm_c, :Float32   ),
             return nothing
         end
 
-        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n')
+        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Vector{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_min_norm_semi_normal(spmat, b, transp = transp)
+
             nrhs = 1
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n)
@@ -704,7 +710,9 @@ for (fname, elty) in ((:sqrm_min_norm_c, :Float32   ),
             return x
         end
 
-        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n')
+        function qrm_min_norm(spmat :: qrm_spmat{$elty}, b :: Matrix{$elty}; transp :: Char='n', seminormal :: Bool = false)
+            seminormal && return qrm_min_norm_semi_normal(spmat, b, transp = transp)
+
             nrhs = size(b, 2)
             if transp == 'n'
                 x = zeros($elty, spmat.mat.n, nrhs)
@@ -719,14 +727,14 @@ for (fname, elty) in ((:sqrm_min_norm_c, :Float32   ),
         @inline qrm_min_norm!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='t')
         @inline qrm_min_norm!(spmat :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='t')
 
-        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_min_norm(spmat.parent, b, transp='t')
-        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_min_norm(spmat.parent, b, transp='t')
+        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='t', seminormal = seminormal)
+        @inline qrm_min_norm(spmat  :: Transpose{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='t', seminormal = seminormal)
 
         @inline qrm_min_norm!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}, x :: Vector{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='c')
         @inline qrm_min_norm!(spmat :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}, x :: Matrix{$elty}) = qrm_min_norm!(spmat.parent, b, x, transp='c')
 
-        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}) = qrm_min_norm(spmat.parent, b, transp='c')
-        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}) = qrm_min_norm(spmat.parent, b, transp='c')
+        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Vector{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='c', seminormal = seminormal)
+        @inline qrm_min_norm(spmat  :: Adjoint{$elty,qrm_spmat{$elty}}, b :: Matrix{$elty}; seminormal :: Bool = false) = qrm_min_norm(spmat.parent, b, transp='c', seminormal = seminormal)
 
     end
 end
diff --git a/test/test_qrm.jl b/test/test_qrm.jl
index ac3fb6c..f57ee58 100644
--- a/test/test_qrm.jl
+++ b/test/test_qrm.jl
@@ -12,6 +12,10 @@ p = 5
     for I in (Int32 , Int64)
       A = sprand(T, m, n, 0.3)
       A = convert(SparseMatrixCSC{T,I}, A)
+
+      A_transp = sprand(T, n, m, 0.3)
+      A_transp = convert(SparseMatrixCSC{T,I}, A_transp)
+
       b = rand(T, m)
       B = rand(T, m, p)
 
@@ -43,6 +47,9 @@ p = 5
       spmat = qrm_spmat_init(T)
       qrm_spmat_init!(spmat, A)
 
+      spmat_transp = qrm_spmat_init(T)
+      qrm_spmat_init!(spmat_transp, A_transp)
+
       spfct = qrm_spfct_init(spmat)
       qrm_analyse!(spmat, spfct)
       qrm_factorize!(spmat, spfct)
@@ -83,20 +90,58 @@ p = 5
       r = b - A * x
       @test norm(A' * r) ≤ tol
 
+      x = qrm_least_squares(spmat, b, seminormal = true)
+      r = b - A * x
+      @test norm(A' * r) ≤ tol
+
+      x = qrm_least_squares_semi_normal(spmat, b)
+      r = b - A * x
+      @test norm(A' * r) ≤ tol
+
+      x = qrm_least_squares_semi_normal(spmat_transp, b, transp = transp)
+      r = b - A_transp' * x
+      @test norm(A_transp * r) ≤ tol
+
       X = qrm_least_squares(spmat, B)
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      X = qrm_least_squares(spmat, B, seminormal = true)
+      R = B - A * X
+      @test norm(A' * R) ≤ tol
+
+      X = qrm_least_squares_semi_normal(spmat, B)
+      R = B - A * X
+      @test norm(A' * R) ≤ tol
+
+      X = qrm_least_squares_semi_normal(spmat_transp, B, transp = transp)
+      R = B - A_transp' * X
+      @test norm(A_transp * r) ≤ tol
+
       bc = copy(b)
       qrm_least_squares!(spmat, bc, x)
       r = b - A * x
       @test norm(A' * r) ≤ tol
 
+      z = similar(x)
+      Δx = similar(x)
+      y = similar(b)
+      x = qrm_least_squares_semi_normal!(spmat, spfct, b, x, z, Δx, y)
+      r = b - A * x
+      @test norm(A' * r) ≤ tol
+
       Bc = copy(B)
       qrm_least_squares!(spmat, Bc, X)
       R = B - A * X
       @test norm(A' * R) ≤ tol
 
+      Z = similar(X)
+      ΔX = similar(X)
+      Y = similar(B)
+      qrm_least_squares_semi_normal!(spmat, spfct, B, X, Z, ΔX, Y)
+      R = B - A * X
+      @test norm(A' * R) ≤ tol
+
       x = qrm_spbackslash(spmat, b)
       r = b - A * x
       @test norm(A' * r) ≤ tol
@@ -161,10 +206,19 @@ end
     for I in (Int32 , Int64)
       A = sprand(T, n, m, 0.3)
       A = convert(SparseMatrixCSC{T,I}, A)
+
+      A_transp = sprand(T, m, n, 0.3)
+      A_transp = convert(SparseMatrixCSC{T,I}, A_transp)
+
       b = rand(T, n)
       B = rand(T, n, p)
+
       spmat = qrm_spmat_init(A)
       spfct = qrm_analyse(spmat, transp=transp)
+
+      spmat_transp = qrm_spmat_init(T)
+      qrm_spmat_init!(spmat_transp, A_transp)
+
       qrm_factorize!(spmat, spfct, transp=transp)
       spfct2 = (T <: Real) ? Transpose(spfct) : Adjoint(spfct)
 
@@ -227,18 +281,54 @@ end
       r = b - A * x
       @test norm(r) ≤ tol
 
+      x = qrm_min_norm(spmat, b, seminormal = true)
+      r = b - A * x
+      @test norm(r) ≤ tol
+
+      x = qrm_min_norm_semi_normal(spmat, b)
+      r = b - A * x
+      @test norm(r) ≤ tol
+
+      x = qrm_min_norm_semi_normal(spmat_transp, b, transp = transp)
+      r = b - A_transp' * x
+      @test norm(r) ≤ tol
+
       X = qrm_min_norm(spmat, B)
       R = B - A * X
       @test norm(R) ≤ tol
 
+      X = qrm_min_norm(spmat, B, seminormal = true)
+      R = B - A * X
+      @test norm(R) ≤ tol
+      
+      X = qrm_min_norm_semi_normal(spmat, B)
+      R = B - A * X
+      @test norm(R) ≤ tol
+
+      X = qrm_min_norm_semi_normal(spmat_transp, B, transp = transp)
+      R = B - A_transp' * X
+      @test norm(R) ≤ tol
+
       qrm_min_norm!(spmat, b, x)
       r = b - A * x
       @test norm(r) ≤ tol
+      
+      Δx = similar(x)
+      y = similar(b)
+      qrm_min_norm_semi_normal!(spmat, spfct, b, x, Δx, y)
+      r = b - A * x
+      @test norm(r) ≤ tol
 
       qrm_min_norm!(spmat, B, X)
       R = B - A * X
       @test norm(R) ≤ tol
 
+      ΔX = similar(X)
+      Y = similar(B)
+      qrm_min_norm_semi_normal!(spmat, spfct, B, X, ΔX, Y)
+      R = B - A * X
+      @test norm(R) ≤ tol
+
       x = qrm_spbackslash(spmat, b)
       r = b - A * x
       @test norm(r) ≤ tol
@@ -544,5 +634,6 @@ end
     x = qrm_solve(spfct, x₁, transp = 'n')
     x_refined = qrm_refine(spmat, spfct, x, b)
     @test norm(b - A'*(A*x)) ≥ norm(b - A'*(A*x_refined))
+
   end
 end