diff --git a/src/nlesolver_module.F90 b/src/nlesolver_module.F90
index f57b97f..807d78b 100755
--- a/src/nlesolver_module.F90
+++ b/src/nlesolver_module.F90
@@ -124,6 +124,15 @@ module nlesolver_module
         integer  :: nout    = 0     !! output unit for printing
         real(wp) :: damp    = zero  !! damp parameter for LSQR
 
+        integer :: localSize = 0 !! LSMR: Number of vectors for local reorthogonalization:
+                                 !!
+                                 !!  *  0    No reorthogonalization is performed.
+                                 !!  * >0    This many n-vectors "v" (the most recent ones)
+                                 !!          are saved for reorthogonalizing the next v.
+                                 !!
+                                 !! localSize need not be more than `min(m,n)`.
+                                 !! At most `min(m,n)` vectors will be allocated.
+
         ! LUSOL parameters:
         integer :: lusol_method = 0 !! * 0 => TPP: Threshold Partial   Pivoting.
                                     !! * 1 => TRP: Threshold Rook      Pivoting.
@@ -323,7 +332,7 @@ subroutine initialize_nlesolver_variables(me,&
                     verbose,iunit,n_uphill_max,n_intervals,&
                     sparsity_mode,irow,icol,&
                     atol,btol,conlim,damp,itnlim,nout,&
-                    lusol_method )
+                    lusol_method,localSize )
 
     implicit none
 
@@ -387,6 +396,14 @@ subroutine initialize_nlesolver_variables(me,&
                                                       !! * 0 => TPP: Threshold Partial   Pivoting.
                                                       !! * 1 => TRP: Threshold Rook      Pivoting.
                                                       !! * 2 => TCP: Threshold Complete  Pivoting.
+    integer,intent(in),optional :: localSize !! `LSMR`: Number of vectors for local reorthogonalization:
+                                             !!
+                                             !!  *  0    No reorthogonalization is performed.
+                                             !!  * >0    This many n-vectors "v" (the most recent ones)
+                                             !!          are saved for reorthogonalizing the next v.
+                                             !!
+                                             !! localSize need not be more than `min(m,n)`.
+                                             !! At most `min(m,n)` vectors will be allocated.
 
     logical :: status_ok !! true if there were no errors
 
@@ -503,6 +520,9 @@ subroutine initialize_nlesolver_variables(me,&
             ! LUSOL method
             if (present(nout)) me%lusol_method = lusol_method
 
+            ! LSMR optional inputs:
+            if (present(localSize)) me%localSize = localSize
+
         end if
     end if
 
@@ -553,7 +573,7 @@ subroutine nlesolver_solver(me,x)
     integer,dimension(:),allocatable :: idx, index_array !! for sparse indexing
     character(len=10) :: i_str !! string version of `i` for row string
     real(wp) :: normA, condA, normr, normAr, normx !! for LSMR
-    integer :: localSize, itn !! for LSMR
+    integer :: itn !! for LSMR
 
     if (me%istat<0) return ! class was not initialized properly
 
@@ -753,15 +773,10 @@ subroutine nlesolver_solver(me,x)
                            settings=lusol_options)
             case (4)
 
-                ! use LSMR solver
-
-                ! TODO this should be an input
-                localSize = 0
-                !localSize = min(me%m, me%n)
-
-                me%conlim = 1.0_wp/(100*epsilon(1.0_wp))
+                ! use LSMR solver:
+                !me%conlim = 1.0_wp/(100*epsilon(1.0_wp))
                 call lsmr_ez  ( me%m, me%n, me%irow, me%icol, fjac_sparse, rhs, me%damp, &
-                                me%atol, me%btol, me%conlim, me%itnlim, localSize, me%nout, &
+                                me%atol, me%btol, me%conlim, me%itnlim, me%localSize, me%nout, &
                                 p, info, itn, normA, condA, normr, normAr, normx )
 
                 ! check convergence:
@@ -778,6 +793,8 @@ subroutine nlesolver_solver(me,x)
                     info = 0
                 end select
 
+            case default
+                error stop 'invalid sparsity_mode'
             end select
 
             ! check for errors: