diff --git a/src/QRMumps.jl b/src/QRMumps.jl
index f6f2d69..b9976b6 100644
--- a/src/QRMumps.jl
+++ b/src/QRMumps.jl
@@ -557,7 +557,7 @@ function qrm_update_shift_spmat! end
 function qrm_golub_riley end 
 
 @doc raw"""
-    qrm_golub_riley!(shifted_spmat, spfct, x, b, Δx, y; α = ϵm, max_iter = 50, tol = ϵm, transp = 'n')
+    qrm_golub_riley!(shifted_spmat, spfct, x, b, Δx, y, Δy; α = ϵm, max_iter = 3, tol = ϵm, transp = 'n')
 
 This method implements the Golub-Riley iteration.
 Given a (possibly ill-conditionned or rank deficient) system `Ax = b` where `A` can have any shape `m×n`, compute `x = A†b = Aᵀ(AAᵀ)†b` where `A†` is the Moore-Penrose pseudoinverse of `A`.
diff --git a/src/utils.jl b/src/utils.jl
index c8b8ced..119b899 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -214,7 +214,7 @@ function qrm_update_shift_spmat!(shifted_spmat :: qrm_shifted_spmat{T}, α :: T)
   shifted_spmat.spmat.val[shifted_spmat.spmat.mat.nz - shifted_spmat.spmat.mat.m + 1:end] .= sqrt(α)
 end
 
-function qrm_golub_riley(spmat :: qrm_spmat{T}, b :: AbstractVector{T}; α :: T = T(eps(real(T))), max_iter :: Int = 50, tol :: Real = eps(real(T)), transp :: Char = 'n') where T
+function qrm_golub_riley(spmat :: qrm_spmat{T}, b :: AbstractVector{T}; α :: T = T(eps(real(T))), max_iter :: Int = 3, tol :: Real = eps(real(T)), transp :: Char = 'n') where T
   shifted_spmat = qrm_shift_spmat(spmat, α)
   spfct = qrm_spfct_init(shifted_spmat.spmat)
   n = shifted_spmat.spmat.mat.n
@@ -250,7 +250,7 @@ function qrm_golub_riley!(
   y :: AbstractVector{T},
   Δy :: AbstractVector{T};  
   α :: T = T(eps(real(T))), 
-  max_iter :: Int = 50, 
+  max_iter :: Int = 3, 
   tol :: Real = eps(real(T)),
   transp :: Char = 'n'
   ) where T