diff --git a/src/TestFunctions.jl b/src/TestFunctions.jl
index 6e444a769eb44973ba1642b1461cf7c50e9824b9..120542cadd252cdd557f87f5cd0e00694072f0f0 100644
--- a/src/TestFunctions.jl
+++ b/src/TestFunctions.jl
@@ -49,8 +49,8 @@ end
 Evaluates the Rosenbrock function.
 """
 function rosenbrock!(fval::Vector{Float64}, x::Vector{Float64})
-    fval[1] = ONE-x[1]
-    fval[2] = TEN*(x[2]-fval[1]*fval[1])
+    @inbounds fval[1] = ONE-x[1]
+    @inbounds fval[2] = TEN*(x[2]-fval[1]*fval[1])
 end
 
 """
@@ -59,10 +59,10 @@ end
 Evaluates the jacobian matrix of the Rosenbrock function.
 """
 function rosenbrock!(fjac::Matrix{Float64}, x::Vector{Float64})
-    fjac[1,1] = -ONE
-    fjac[1,2] = ZERO
-    fjac[2,1] = -TWENTY*x[1]
-    fjac[2,2] = TEN
+    @inbounds fjac[1,1] = -ONE
+    @inbounds fjac[1,2] = ZERO
+    @inbounds fjac[2,1] = -TWENTY*x[1]
+    @inbounds fjac[2,2] = TEN
 end
 
 """
@@ -75,24 +75,24 @@ function rosenbrock()
 end
 
 function powell1!(fval::Vector{Float64}, x::Vector{Float64})
-    fval[1] = x[1]+TEN*x[2]
-    fval[2] = sqrt(FIVE)*(x[3]-x[4])
-    tmp = x[2]-TWO*x[3]
-    fval[3] = tmp*tmp
-    tmp = x[1]-x[4]
-    fval[4] = sqrt(TEN)*tmp*tmp
+    @inbounds fval[1] = x[1]+TEN*x[2]
+    @inbounds fval[2] = sqrt(FIVE)*(x[3]-x[4])
+    @inbounds tmp = x[2]-TWO*x[3]
+    @inbounds fval[3] = tmp*tmp
+    @inbounds tmp = x[1]-x[4]
+    @inbounds fval[4] = sqrt(TEN)*tmp*tmp
 end
 
 function powell1!(fjac::Matrix{Float64}, x::Vector{Float64})
     fill!(fjac, 0.0)
-    fjac[1,1] = ONE
-    fjac[1,2] = TEN
-    fjac[2,3] = sqrt(FIVE)
-    fjac[2,4] = -fjac[2,3]
-    fjac[3,2] = TWO*(x[2]-TWO*x[3])
-    fjac[3,3] = -TWO*fjac[3,2]
-    fjac[4,1] = TWO*sqrt(TEN)*(x[1]-x[4])
-    fjac[4,4] = -fjac[4,1]
+    @inbounds fjac[1,1] = ONE
+    @inbounds fjac[1,2] = TEN
+    @inbounds fjac[2,3] = sqrt(FIVE)
+    @inbounds fjac[2,4] = -fjac[2,3]
+    @inbounds fjac[3,2] = TWO*(x[2]-TWO*x[3])
+    @inbounds fjac[3,3] = -TWO*fjac[3,2]
+    @inbounds fjac[4,1] = TWO*sqrt(TEN)*(x[1]-x[4])
+    @inbounds fjac[4,4] = -fjac[4,1]
 end
 
 function powell1()
@@ -100,15 +100,15 @@ function powell1()
 end
 
 function powell2!(fval::Vector{Float64}, x::Vector{Float64})
-    fval[1] = C1*x[1]*x[2]-ONE
-    fval[2] = exp(-x[1])+exp(-x[2])-C2
+    @inbounds fval[1] = C1*x[1]*x[2]-ONE
+    @inbounds fval[2] = exp(-x[1])+exp(-x[2])-C2
 end
 
 function powell2!(fjac::Matrix{Float64}, x::Vector{Float64})
-    fjac[1,1] = C1*x[2]
-    fjac[1,2] = C1*x[1]
-    fjac[2,1] = -exp(-x[1])
-    fjac[2,2] = -exp(-x[2])
+    @inbounds fjac[1,1] = C1*x[2]
+    @inbounds fjac[1,2] = C1*x[1]
+    @inbounds fjac[2,1] = -exp(-x[1])
+    @inbounds fjac[2,2] = -exp(-x[2])
 end
 
 function powell2()
@@ -116,28 +116,28 @@ function powell2()
 end
 
 function wood!(fval::Vector{Float64}, x::Vector{Float64})
-    tmp1 = x[2]-x[1]*x[1]
-    tmp2 = x[4]-x[3]*x[3]
-    fval[1] = -C3*x[1]*tmp1-(ONE-x[1])
-    fval[2] = C3*tmp1+C4*(x[2]-1.0)+C5*(x[4]-ONE)
-    fval[3] = -C6*x[3]*tmp2-(ONE-x[3])
-    fval[4] = C6*tmp2+C4*(x[4]-ONE)+C5*(x[2]-ONE)
+    @inbounds tmp1 = x[2]-x[1]*x[1]
+    @inbounds tmp2 = x[4]-x[3]*x[3]
+    @inbounds fval[1] = -C3*x[1]*tmp1-(ONE-x[1])
+    @inbounds fval[2] = C3*tmp1+C4*(x[2]-1.0)+C5*(x[4]-ONE)
+    @inbounds fval[3] = -C6*x[3]*tmp2-(ONE-x[3])
+    @inbounds fval[4] = C6*tmp2+C4*(x[4]-ONE)+C5*(x[2]-ONE)
 end
 
 function wood!(fjac::Matrix{Float64}, x::Vector{Float64})
     fill!(fjac, ZERO)
-    tmp1 = x[2]-THREE*x[1]*x[1]
-    tmp2 = x[4]-THREE*x[3]*x[3]
-    fjac[1,1] = -C3*tmp1+ONE
-    fjac[1,2] = -C3*x[1]
-    fjac[2,1] = -TWO*C3*x[1]
-    fjac[2,2] = C3+C4
-    fjac[2,4] = C5
-    fjac[3,3] = -C6*tmp2+ONE
-    fjac[3,4] = -C6*x[3]
-    fjac[4,2] = C5
-    fjac[4,3] = -TWO*C6*x[3]
-    fjac[4,4] = C6+C4
+    @inbounds tmp1 = x[2]-THREE*x[1]*x[1]
+    @inbounds tmp2 = x[4]-THREE*x[3]*x[3]
+    @inbounds fjac[1,1] = -C3*tmp1+ONE
+    @inbounds fjac[1,2] = -C3*x[1]
+    @inbounds fjac[2,1] = -TWO*C3*x[1]
+    @inbounds fjac[2,2] = C3+C4
+    @inbounds fjac[2,4] = C5
+    @inbounds fjac[3,3] = -C6*tmp2+ONE
+    @inbounds fjac[3,4] = -C6*x[3]
+    @inbounds fjac[4,2] = C5
+    @inbounds fjac[4,3] = -TWO*C6*x[3]
+    @inbounds fjac[4,4] = C6+C4
 end
 
 function wood()
@@ -145,32 +145,32 @@ function wood()
 end
 
 function helicalvalley!(fval::Vector{Float64}, x::Vector{Float64})
-    tmp1 = sign(C7, x[2])
+    @inbounds tmp1 = sign(C7, x[2])
     if x[1]>ZERO
-        tmp1 = atan(x[2]/x[1])/TPI
+        @inbounds tmp1 = atan(x[2]/x[1])/TPI
     end
     if x[1]<ZERO
-        tmp1 = atan(x[2]/x[1])/TPI+C8
+        @inbounds tmp1 = atan(x[2]/x[1])/TPI+C8
     end
-    tmp2 = sqrt(x[1]*x[1]+x[2]*x[2])
-    fval[1] = TEN*(x[3]-TEN*tmp1)
-    fval[2] = TEN*(tmp2-ONE)
-    fval[3] = x[3]
+    @inbounds tmp2 = sqrt(x[1]*x[1]+x[2]*x[2])
+    @inbounds fval[1] = TEN*(x[3]-TEN*tmp1)
+    @inbounds fval[2] = TEN*(tmp2-ONE)
+    @inbounds fval[3] = x[3]
 end
 
 function helicalvalley!(fjac::Matrix{Float64}, x::Vector{Float64})
-    tmp = x[1]*x[1]+x[2]*x[2]
+    @inbounds tmp = x[1]*x[1]+x[2]*x[2]
     tmp1 = TPI*tmp
     tmp2 = sqrt(tmp)
-    fjac[1,1] = HUNDRD*x[2]/tmp1
-    fjac[1,2] = -HUNDRD*x[1]/tmp1
-    fjac[1,3] = TEN
-    fjac[2,1] = TEN*x[1]/tmp2
-    fjac[2,2] = TEN*x[2]/tmp2
-    fjac[2,3] = ZERO
-    fjac[3,1] = ZERO
-    fjac[3,2] = ZERO
-    fjac[3,3] = ONE
+    @inbounds fjac[1,1] = HUNDRD*x[2]/tmp1
+    @inbounds fjac[1,2] = -HUNDRD*x[1]/tmp1
+    @inbounds fjac[1,3] = TEN
+    @inbounds fjac[2,1] = TEN*x[1]/tmp2
+    @inbounds fjac[2,2] = TEN*x[2]/tmp2
+    @inbounds fjac[2,3] = ZERO
+    @inbounds fjac[3,1] = ZERO
+    @inbounds fjac[3,2] = ZERO
+    @inbounds fjac[3,3] = ONE
 end
 
 function helicalvalley()
@@ -185,26 +185,26 @@ function watson!(fval::Vector{Float64}, x::Vector{Float64})
         sum1 = ZERO
         tmp = ONE
         for j=2:n
-            sum1 += Float64(j-1)*tmp*x[j]
+            @inbounds sum1 += Float64(j-1)*tmp*x[j]
             tmp *= ti
         end
         sum2 = ZERO
         tmp = ONE
         for j=1:n
-            sum2 += tmp*x[j]
+            @inbounds sum2 += tmp*x[j]
             tmp *= ti
         end
         tmp1 = sum1-sum2*sum2-ONE
         tmp2 = TWO*ti*sum2
         tmp = ONE/ti
         for k=1:n
-            fval[k] += tmp*(Float64(k-1)-tmp2)*tmp1
+            @inbounds fval[k] += tmp*(Float64(k-1)-tmp2)*tmp1
             tmp *= ti
         end
     end
-    tmp = x[2]-x[1]*x[1]-ONE
-    fval[1] += x[1]*(ONE-TWO*tmp)
-    fval[2] += tmp
+    @inbounds tmp = x[2]-x[1]*x[1]-ONE
+    @inbounds fval[1] += x[1]*(ONE-TWO*tmp)
+    @inbounds fval[2] += tmp
 end
 
 function watson!(fjac::Matrix{Float64}, x::Vector{Float64})
@@ -215,13 +215,13 @@ function watson!(fjac::Matrix{Float64}, x::Vector{Float64})
         sum1 = ZERO
         tmp = ONE
         for j=2:n
-            sum1 += Float64(j-1)*tmp*x[j]
+            @inbounds sum1 += Float64(j-1)*tmp*x[j]
             tmp *= ti
         end
         sum2 = ZERO
         tmp = ONE
         for j=1:n
-            sum2 += tmp*x[j]
+            @inbounds sum2 += tmp*x[j]
             tmp *= ti
         end
         tmp1 = TWO*(sum1-sum2*sum2-ONE)
@@ -231,18 +231,18 @@ function watson!(fjac::Matrix{Float64}, x::Vector{Float64})
         for k=1:n
             tj = tk
             for j=k:n
-                fjac[k,j] += tj*((Float64(k-1)/ti-tmp2)*(Float64(j-1)/ti-tmp2)-tmp1)
+                @inbounds fjac[k,j] += tj*((Float64(k-1)/ti-tmp2)*(Float64(j-1)/ti-tmp2)-tmp1)
                 tj *= ti
             end
             tk *= tmp
         end
     end
-    fjac[1,1] = fjac[1,1]+SIX*x[1]*x[1]-TWO*x[2]+THREE
-    fjac[1,2] = fjac[1,2]-TWO*x[1]
-    fjac[2,2] = fjac[2,2]+ONE
+    @inbounds fjac[1,1] = fjac[1,1]+SIX*x[1]*x[1]-TWO*x[2]+THREE
+    @inbounds fjac[1,2] = fjac[1,2]-TWO*x[1]
+    @inbounds fjac[2,2] = fjac[2,2]+ONE
     for k=1:n
         for j=k:n
-            fjac[j,k] = fjac[k,j]
+            @inbounds fjac[j,k] = fjac[k,j]
         end
     end
 end
@@ -256,10 +256,10 @@ function chebyquad!(fval::Vector{Float64}, x::Vector{Float64})
     fill!(fval, ZERO)
     for j=1:n
         tmp1 = ONE
-        tmp2 = TWO*x[j]-ONE
+        @inbounds tmp2 = TWO*x[j]-ONE
         tmp = TWO*tmp2
         for i=1:n
-            fval[i] += tmp2
+            @inbounds fval[i] += tmp2
             ti = tmp*tmp2-tmp1
             tmp1 = tmp2
             tmp2 = ti
@@ -268,9 +268,9 @@ function chebyquad!(fval::Vector{Float64}, x::Vector{Float64})
     tk = ONE/Float64(n)
     iev = -1
     for k=1:n
-        fval[k] *= tk
+        @inbounds fval[k] *= tk
         if iev>0
-            fval[k] += ONE/(Float64(k)^2-ONE)
+            @inbounds fval[k] += ONE/(Float64(k)^2-ONE)
         end
         iev = -iev
     end
@@ -281,12 +281,12 @@ function chebyquad!(fjac::Matrix{Float64}, x::Vector{Float64})
     tk = ONE/Float64(n)
     for j=1:n
         tmp1 = ONE
-        tmp2 = TWO*x[j]-ONE
+        @inbounds tmp2 = TWO*x[j]-ONE
         tmp = TWO*tmp2
         tmp3 = ZERO
         tmp4 = TWO
         for k=1:n
-            fjac[k,j] = tk*tmp4
+            @inbounds fjac[k,j] = tk*tmp4
             ti = FOUR*tmp2+tmp*tmp4-tmp3
             tmp3 = tmp4
             tmp4 = ti
@@ -301,7 +301,7 @@ function chebyquad(n::Int)
     h = ONE/Float64(n+1)
     x = Vector{Float64}(undef,n)
     for j=1:n
-        x[j] = Float64(j)*h
+        @inbounds x[j] = Float64(j)*h
     end
     return x
 end
@@ -310,38 +310,38 @@ function brown!(fval::Vector{Float64}, x::Vector{Float64})
     n = length(fval)
     sum = -Float64(n+1)
     prod = ONE
-    for j=1:n
+    @inbounds for j=1:n
         sum += x[j]
         prod *= x[j]
     end
     for k=1:n-1
-        fval[k]=x[k]+sum
+        @inbounds fval[k]=x[k]+sum
     end
-    fval[n] = prod-ONE
+    @inbounds fval[n] = prod-ONE
 end
 
 function brown!(fjac::Matrix{Float64}, x::Vector{Float64})
     n = length(x)
     prod = ONE
     for j=1:n
-        prod *= x[j]
+        @inbounds prod *= x[j]
         for k=1:n
-            fjac[k,j] = ONE
+            @inbounds fjac[k,j] = ONE
         end
-        fjac[j,j] = TWO
+        @inbounds fjac[j,j] = TWO
     end
     for j=1:n
-        tmp = x[j]
+        @inbounds tmp = x[j]
         if abs(tmp)<eps(Float64)
             tmp = ONE
             prod = ONE
             for k=1:n
                 if k!=j
-                    prod *= x[k]
+                    @inbounds prod *= x[k]
                 end
             end
         end
-        fjac[n,j] = prod/tmp
+        @inbounds fjac[n,j] = prod/tmp
     end
 end
 
@@ -359,13 +359,13 @@ function discreteboundaryvalue!(fval::Vector{Float64}, x::Vector{Float64})
         tmp = (x[k]+Float64(k)*h+ONE)^3
         tmp1 = ZERO
         if k!=1
-            tmp1 = x[k-1]
+            @inbounds tmp1 = x[k-1]
         end
         tmp2 = ZERO
         if k!=n
-            tmp2 = x[k+1]
+            @inbounds tmp2 = x[k+1]
         end
-        fval[k] = TWO*x[k]-tmp1-tmp2+tmp*hh/TWO
+        @inbounds fval[k] = TWO*x[k]-tmp1-tmp2+tmp*hh/TWO
     end
 end
 
@@ -375,13 +375,13 @@ function discreteboundaryvalue!(fjac::Matrix{Float64}, x::Vector{Float64})
     hh = h*h
     fill!(fjac, ZERO)
     for k=1:n
-        tmp = THREE*(x[k]+Float64(k)*h+ONE)^2
-        fjac[k,k] = TWO+tmp*hh/TWO
+        @inbounds tmp = THREE*(x[k]+Float64(k)*h+ONE)^2
+        @inbounds fjac[k,k] = TWO+tmp*hh/TWO
         if k!=1
-            fjac[k,k-1] = -ONE
+            @inbounds fjac[k,k-1] = -ONE
         end
         if k!=n
-            fjac[k,k+1] = -ONE
+            @inbounds fjac[k,k+1] = -ONE
         end
     end
 end
@@ -391,7 +391,7 @@ function discreteboundaryvalue(n::Int)
     x = Vector{Float64}(undef,n)
     for j=1:n
         tj = Float64(j)*h
-        x[j] = tj*(tj-ONE)
+        @inbounds x[j] = tj*(tj-ONE)
     end
     return x
 end
@@ -404,7 +404,7 @@ function discreteintegralequation!(fval::Vector{Float64}, x::Vector{Float64})
         sum1 = ZERO
         for j=1:k
             tj = Float64(j)*h
-            tmp = (x[j]+tj+ONE)^3
+            @inbounds tmp = (x[j]+tj+ONE)^3
             sum1 += tj*tmp
         end
         sum2 = ZERO
@@ -412,11 +412,11 @@ function discreteintegralequation!(fval::Vector{Float64}, x::Vector{Float64})
         if n>=kp1
             for j=kp1:n
                 tj = Float64(j)*h
-                tmp = (x[j]+tj+ONE)^3
+                @inbounds tmp = (x[j]+tj+ONE)^3
                 sum2 += (ONE-tj)*tmp
             end
         end
-        fval[k] = x[k] + h*((ONE-tk)*sum1+tk*sum2)/TWO
+        @inbounds fval[k] = x[k] + h*((ONE-tk)*sum1+tk*sum2)/TWO
     end
 end
 
@@ -427,10 +427,10 @@ function discreteintegralequation!(fjac::Matrix{Float64}, x::Vector{Float64})
         tk = Float64(k)*h
         for j=1:n
             tj = Float64(j)*h
-            tmp = THREE*(x[j]+tj+ONE)^2
-            fjac[k,j] = h*min(tj*(ONE-tk),tk*(ONE-tj))*tmp/TWO
+            @inbounds tmp = THREE*(x[j]+tj+ONE)^2
+            @inbounds fjac[k,j] = h*min(tj*(ONE-tk),tk*(ONE-tj))*tmp/TWO
         end
-        fjac[k,k] += ONE
+        @inbounds fjac[k,k] += ONE
     end
 end
 
@@ -440,22 +440,22 @@ function trigonometric!(fval::Vector{Float64}, x::Vector{Float64})
     n = length(fval)
     sum = ZERO
     for j=1:n
-        fval[j] = cos(x[j])
-        sum += fval[j]
+        @inbounds fval[j] = cos(x[j])
+        @inbounds sum += fval[j]
     end
     for k=1:n
-        fval[k] = Float64(n+k)-sin(x[k])-sum-Float64(k)*fval[k]
+        @inbounds fval[k] = Float64(n+k)-sin(x[k])-sum-Float64(k)*fval[k]
     end
 end
 
 function trigonometric!(fjac::Matrix{Float64}, x::Vector{Float64})
     n = length(x)
     for j=1:n
-        tmp = sin(x[j])
+        @inbounds tmp = sin(x[j])
         for k=1:n
-            fjac[k,j] = tmp
+            @inbounds fjac[k,j] = tmp
         end
-        fjac[j,j] = Float64(j+1)*tmp-cos(x[j])
+        @inbounds fjac[j,j] = Float64(j+1)*tmp-cos(x[j])
     end
 end
 
@@ -469,11 +469,11 @@ function variablydimensioned!(fval::Vector{Float64}, x::Vector{Float64})
     n = length(fval)
     sum = ZERO
     for j=1:n
-        sum += Float64(j)*(x[j]-ONE)
+        @inbounds sum += Float64(j)*(x[j]-ONE)
     end
     tmp = sum*(ONE+TWO*sum*sum)
     for k=1:n
-        fval[k] = x[k]-ONE+Float64(k)*tmp
+        @inbounds fval[k] = x[k]-ONE+Float64(k)*tmp
     end
 end
 
@@ -481,15 +481,15 @@ function variablydimensioned!(fjac::Matrix{Float64}, x::Vector{Float64})
     n = length(x)
     sum = ZERO
     for j=1:n
-        sum += Float64(j)*(x[j]-ONE)
+        @inbounds sum += Float64(j)*(x[j]-ONE)
     end
     tmp = ONE+SIX*sum*sum
     for k=1:n
         for j=1:n
-            fjac[k,j] = Float64(k*j)*tmp
-            fjac[j,k] = fjac[k,j]
+            @inbounds fjac[k,j] = Float64(k*j)*tmp
+            @inbounds fjac[j,k] = fjac[k,j]
         end
-        fjac[k,k] += ONE
+        @inbounds fjac[k,k] += ONE
     end
 end
 
@@ -497,7 +497,7 @@ function variablydimensioned(n::Int)
     h = ONE/Float64(n)
     x = Vector{Float64}(undef,n)
     for j=1:n
-        x[j] = ONE-Float64(j)*h
+        @inbounds x[j] = ONE-Float64(j)*h
     end
     return x
 end
@@ -505,16 +505,16 @@ end
 function broydentridiagonal!(fval::Vector{Float64}, x::Vector{Float64})
     n = length(fval)
     for k=1:n
-        tmp = (THREE-TWO*x[k])*x[k]
+        @inbounds tmp = (THREE-TWO*x[k])*x[k]
         tmp1 = ZERO
         if k!=1
-            tmp1 = x[k-1]
+            @inbounds tmp1 = x[k-1]
         end
         tmp2 = ZERO
         if k!=n
-            tmp2 = x[k+1]
+            @inbounds tmp2 = x[k+1]
         end
-        fval[k] = tmp-tmp1-TWO*tmp2+ONE
+        @inbounds fval[k] = tmp-tmp1-TWO*tmp2+ONE
     end
 end
 
@@ -522,12 +522,12 @@ function broydentridiagonal!(fjac::Matrix{Float64}, x::Vector{Float64})
     n = length(x)
     fill!(fjac, ZERO)
     for k=1:n
-        fjac[k,k] = THREE-FOUR*x[k]
+        @inbounds fjac[k,k] = THREE-FOUR*x[k]
         if k!=1
-            fjac[k,k-1] = -ONE
+            @inbounds fjac[k,k-1] = -ONE
         end
         if k!=n
-            fjac[k,k+1] = -TWO
+            @inbounds fjac[k,k+1] = -TWO
         end
     end
 end
@@ -548,10 +548,10 @@ function broydenbanded!(fval::Vector{Float64}, x::Vector{Float64})
         tmp = ZERO
         for j=k1:k2
             if j!=k
-                tmp += x[j]*(ONE+x[j])
+                @inbounds tmp += x[j]*(ONE+x[j])
             end
         end
-        fval[k] = x[k]*(TWO+FIVE*x[k]*x[k])+ONE-tmp
+        @inbounds fval[k] = x[k]*(TWO+FIVE*x[k]*x[k])+ONE-tmp
     end
 end
 
@@ -565,10 +565,10 @@ function broydenbanded!(fjac::Matrix{Float64}, x::Vector{Float64})
         k2 = min(k+mu,n)
         for j=k1:k2
             if j!=k
-                fjac[k,j] = -(ONE+TWO*x[j])
+                @inbounds fjac[k,j] = -(ONE+TWO*x[j])
             end
         end
-        fjac[k,k] = TWO+FIFTN*x[k]*x[k]
+        @inbounds fjac[k,k] = TWO+FIFTN*x[k]*x[k]
     end
 end