diff --git a/README.md b/README.md
index 14edbcb..bd22dfa 100644
--- a/README.md
+++ b/README.md
@@ -38,6 +38,8 @@ A [FoBiS](https://github.com/szaghi/FoBiS) configuration file (`slsqp.fobis`) is
 
   To generate the documentation using [ford](https://github.com/cmacmackin/ford), run: ```FoBis.py rule --execute makedoc -f slsqp.fobis```
 
+  To run the test programs, run: ```FoBis.py rule --execute tests -f slsqp.fobis```
+
 ### Development
 
 [![Build Status](https://img.shields.io/travis/jacobwilliams/slsqp/master.svg?style=plastic)](https://travis-ci.org/jacobwilliams/slsqp)
diff --git a/slsqp.fobis b/slsqp.fobis
index cba6803..e995ddf 100644
--- a/slsqp.fobis
+++ b/slsqp.fobis
@@ -142,3 +142,15 @@ template        = template-tests
 [rule-makedoc]
 help   = Rule for building documentation from source files
 rule_1 = ford $FORD_FILE
+
+[rule-tests]
+help   = Rule for running the test programs
+rule_1 = (cd bin
+         GLOBIGNORE='*.*'
+         ls slsqp_test_* | sed 's/^\([^0-9]*\)\([0-9]*\)/\1 \2/' | sort -k2,2n | tr -d ' ' |
+         while read TEST; do
+             echo ""
+             echo "Running ${TEST}"
+             "./${TEST}"
+         done)
+
diff --git a/src/slsqp_module.f90 b/src/slsqp_module.f90
index f6a5331..5bc856f 100644
--- a/src/slsqp_module.f90
+++ b/src/slsqp_module.f90
@@ -291,28 +291,37 @@ subroutine slsqp_wrapper(me,x,istat,iterations,status_message)
                                                       !! string status message
                                                       !! corresponding to `istat`
 
-    !local variables:
-    real(wp)                               :: f        !! objective function
-    real(wp),dimension(max(1,me%m))        :: c        !! constraint vector
-    real(wp),dimension(max(1,me%m),me%n+1) :: a        !! a matrix for [[slsqp]]
-    real(wp),dimension(me%n+1)             :: g        !! g matrix for [[slsqp]]
-    real(wp),dimension(me%m)               :: cvec     !! constraint vector
-    real(wp),dimension(me%n)               :: dfdx     !! objective function partials
-    real(wp),dimension(me%m,me%n)          :: dcdx     !! constraint partials
-    integer                                :: i        !! iteration counter
-    integer                                :: mode     !! reverse communication flag for [[slsqp]]
-    integer                                :: la       !! input to [[slsqp]]
-    integer                                :: iter     !! in/out for [[slsqp]]
-    real(wp)                               :: acc      !! in/out for [[slsqp]]
-    integer                                :: ig       !! loop index to approximate gradient
-    real(wp),dimension(me%n)               :: delta    !! perturbation step size to approximate gradient
-    real(wp)                               :: fr       !! right function value to approximate objective function's gradient
-    real(wp)                               :: fl       !! left function value to approximate objective function's gradient
-    real(wp),dimension(me%m)               :: cvecr    !! right function value to approximate constraints vector's gradient
-    real(wp),dimension(me%m)               :: cvecl    !! left function value to approximate constraints vector's gradient
-    real(wp)                               :: fact     !! denominator factor for finite difference approximation
+    ! local variables:
+    real(wp),dimension(:),allocatable   :: c        !! constraint vector -- `dimension(max(1,me%m))`
+    real(wp),dimension(:,:),allocatable :: a        !! a matrix for [[slsqp]] -- `dimension(max(1,me%m),me%n+1)`
+    real(wp),dimension(:),allocatable   :: g        !! g matrix for [[slsqp]] -- `dimension(me%n+1)`
+    real(wp),dimension(:),allocatable   :: cvec     !! constraint vector -- `dimension(me%m)`
+    real(wp),dimension(:),allocatable   :: dfdx     !! objective function partials -- `dimension(me%n)`
+    real(wp),dimension(:,:),allocatable :: dcdx     !! constraint partials -- `dimension(me%m,me%n)`
+    real(wp),dimension(:),allocatable   :: delta    !! perturbation step size to approximate gradient -- `dimension(me%n)`
+    real(wp),dimension(:),allocatable   :: cvecr    !! right function value to approximate constraints vector's gradient -- `dimension(me%m)`
+    real(wp),dimension(:),allocatable   :: cvecl    !! left function value to approximate constraints vector's gradient -- `dimension(me%m)`
+    real(wp)                            :: f        !! objective function
+    integer                             :: i        !! iteration counter
+    integer                             :: mode     !! reverse communication flag for [[slsqp]]
+    integer                             :: la       !! input to [[slsqp]]
+    integer                             :: iter     !! in/out for [[slsqp]]
+    real(wp)                            :: acc      !! in/out for [[slsqp]]
+    integer                             :: ig       !! loop index to approximate gradient
+    real(wp)                            :: fr       !! right function value to approximate objective function's gradient
+    real(wp)                            :: fl       !! left function value to approximate objective function's gradient
+    real(wp)                            :: fact     !! denominator factor for finite difference approximation
 
     !initialize:
+    allocate(c(max(1,me%m))       )
+    allocate(a(max(1,me%m),me%n+1))
+    allocate(g(me%n+1)            )
+    allocate(cvec(me%m)           )
+    allocate(dfdx(me%n)           )
+    allocate(dcdx(me%m,me%n)      )
+    allocate(delta(me%n)          )
+    allocate(cvecr(me%m)          )
+    allocate(cvecl(me%m)          )
     i    = 0
     iter = me%max_iter
     la   = max(1,me%m)
@@ -376,7 +385,7 @@ subroutine slsqp_wrapper(me,x,istat,iterations,status_message)
 
         if (mode==0 .or. mode==1) then  !function evaluation (f&c)
             call me%f(x,f,cvec)
-            c(1:me%m)   = cvec
+            c(1:me%m) = cvec
         end if
 
         if (mode==0 .or. mode==-1) then  !gradient evaluation (g&a)