diff --git a/fpm.toml b/fpm.toml
index 6d478fc..14159b8 100644
--- a/fpm.toml
+++ b/fpm.toml
@@ -71,3 +71,8 @@ main = "example_lmdif1.f90"
 name = "example_primes"
 source-dir = "examples"
 main = "example_primes.f90"
+
+[[test]]
+name = "c-test"
+source-dir = "test/api"
+main = "tester.c"
diff --git a/include/minpack.h b/include/minpack.h
new file mode 100644
index 0000000..f78a16b
--- /dev/null
+++ b/include/minpack.h
@@ -0,0 +1,363 @@
+#pragma once
+
+#ifdef __cplusplus
+#define MINPACK_EXTERN extern "C"
+#else
+#define MINPACK_EXTERN extern
+#endif
+
+// In case we have to add some attributes to all functions (e.g. for DLL exports)
+#define MINPACK_CALL
+
+MINPACK_EXTERN double MINPACK_CALL
+minpack_dpmpar(int /* i */);
+
+typedef void (*minpack_func)(
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    int* /* iflag */,
+    void* /* udata */);
+
+/*
+ * the purpose of hybrd is to find a zero of a system of
+ * n nonlinear functions in n variables by a modification
+ * of the powell hybrid method. the user must provide a
+ * subroutine which calculates the functions. the jacobian is
+ * then calculated by a forward-difference approximation.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_hybrd(
+    minpack_func /* fcn */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double /* xtol */,
+    int /* maxfev */,
+    int /* ml */,
+    int /* mu */,
+    double /* epsfcn */,
+    double* /* diag */,
+    int /* mode */,
+    double /* factor */,
+    int /* nprint */,
+    int* /* info */,
+    int* /* nfev */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double* /* r */,
+    int /* lr */,
+    double* /* qtf */,
+    double* /* wa1 */,
+    double* /* wa2 */,
+    double* /* wa3 */,
+    double* /* wa4 */,
+    void* /* udata */);
+
+/*
+ * the purpose of hybrd1 is to find a zero of a system of
+ * n nonlinear functions in n variables by a modification
+ * of the powell hybrid method. this is done by using the
+ * more general nonlinear equation solver hybrd. the user
+ * must provide a subroutine which calculates the functions.
+ * the jacobian is then calculated by a forward-difference
+ * approximation.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_hybrd1(
+    minpack_func /* fcn */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double /* tol */,
+    int* /* info */,
+    double* /* wa */,
+    int /* lwa */,
+    void* /* udata */);
+
+typedef void (*minpack_fcn_hybrj)(
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    int* /* iflag */,
+    void* /* udata */);
+
+/*
+ * the purpose of hybrj is to find a zero of a system of
+ * n nonlinear functions in n variables by a modification
+ * of the powell hybrid method. the user must provide a
+ * subroutine which calculates the functions and the jacobian.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_hybrj(
+    minpack_fcn_hybrj /* fcn */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double /* xtol */,
+    int /* maxfev */,
+    double* /* diag */,
+    int /* mode */,
+    double /* factor */,
+    int /* nprint */,
+    int* /* info */,
+    int* /* nfev */,
+    int* /* njev */,
+    double* /* r */,
+    int /* lr */,
+    double* /* qtf */,
+    double* /* wa1 */,
+    double* /* wa2 */,
+    double* /* wa3 */,
+    double* /* wa4 */,
+    void* /* udata */);
+
+/*
+ * The purpose of hybrj1 is to find a zero of a system of
+ * n nonlinear functions in n variables by a modification
+ * of the powell hybrid method. this is done by using the
+ * more general nonlinear equation solver hybrj. the user
+ * must provide a subroutine which calculates the functions
+ * and the jacobian.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_hybrj1(
+    minpack_fcn_hybrj /* fcn */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double /* tol */,
+    int* /* info */,
+    double* /* wa */,
+    int /* lwa */,
+    void* /* udata */);
+
+typedef void (*minpack_fcn_lmder)(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    int* /* iflag */,
+    void* /* udata */);
+
+/*
+ * the purpose of lmder is to minimize the sum of the squares of
+ * m nonlinear functions in n variables by a modification of
+ * the levenberg-marquardt algorithm. the user must provide a
+ * subroutine which calculates the functions and the jacobian.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_lmder(
+    minpack_fcn_lmder /* fcn */,
+    int /* m */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double /* ftol */,
+    double /* xtol */,
+    double /* gtol */,
+    int /* maxfev */,
+    double* /* diag */,
+    int /* mode */,
+    double /* factor */,
+    int /* nprint */,
+    int* /* info */,
+    int* /* nfev */,
+    int* /* njev */,
+    int* /* ipvt */,
+    double* /* qtf */,
+    double* /* wa1 */,
+    double* /* wa2 */,
+    double* /* wa3 */,
+    double* /* wa4 */,
+    void* /* udata */);
+
+/*
+ * the purpose of lmder1 is to minimize the sum of the squares of
+ * m nonlinear functions in n variables by a modification of the
+ * levenberg-marquardt algorithm. this is done by using the more
+ * general least-squares solver lmder. the user must provide a
+ * subroutine which calculates the functions and the jacobian.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_lmder1(
+    minpack_fcn_lmder /* fcn */,
+    int /* m */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double /* tol */,
+    int* /* info */,
+    int* /* ipvt */,
+    double* /* wa */,
+    int /* lwa */,
+    void* /* udata */);
+
+typedef void (*minpack_func2)(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    int* /* iflag */,
+    void* /* udata */);
+
+/*
+ * the purpose of lmdif is to minimize the sum of the squares of
+ * m nonlinear functions in n variables by a modification of
+ * the levenberg-marquardt algorithm. the user must provide a
+ * subroutine which calculates the functions. the jacobian is
+ * then calculated by a forward-difference approximation.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_lmdif(
+    minpack_func2 /* fcn */,
+    int /* m */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double /* ftol */,
+    double /* xtol */,
+    double /* gtol */,
+    int /* maxfev */,
+    double /* epsfcn */,
+    double* /* diag */,
+    int /* mode */,
+    double /* factor */,
+    int /* nprint */,
+    int* /* info */,
+    int* /* nfev */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    int* /* ipvt */,
+    double* /* qtf */,
+    double* /* wa1 */,
+    double* /* wa2 */,
+    double* /* wa3 */,
+    double* /* wa4 */,
+    void* /* udata */);
+
+/*
+ * the purpose of lmdif1 is to minimize the sum of the squares of
+ * m nonlinear functions in n variables by a modification of the
+ * levenberg-marquardt algorithm. this is done by using the more
+ * general least-squares solver lmdif. the user must provide a
+ * subroutine which calculates the functions. the jacobian is
+ * then calculated by a forward-difference approximation.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_lmdif1(
+    minpack_func2 /* fcn */,
+    int /* m */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double /* tol */,
+    int* /* info */,
+    int* /* iwa */,
+    double* /* wa */,
+    int /* lwa */,
+    void* /* udata */);
+
+typedef void (*minpack_fcn_lmstr)(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    double* /* fvec */,
+    double* /* fjrow */,
+    int* /* iflag */,
+    void* /* udata */);
+
+/*
+ * the purpose of lmstr is to minimize the sum of the squares of
+ * m nonlinear functions in n variables by a modification of
+ * the levenberg-marquardt algorithm which uses minimal storage.
+ * the user must provide a subroutine which calculates the
+ * functions and the rows of the jacobian.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_lmstr(
+    minpack_fcn_lmstr /* fcn */,
+    int /* m */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double /* ftol */,
+    double /* xtol */,
+    double /* gtol */,
+    int /* maxfev */,
+    double* /* diag */,
+    int /* mode */,
+    double /* factor */,
+    int /* nprint */,
+    int* /* info */,
+    int* /* nfev */,
+    int* /* njev */,
+    int* /* ipvt */,
+    double* /* qtf */,
+    double* /* wa1 */,
+    double* /* wa2 */,
+    double* /* wa3 */,
+    double* /* wa4 */,
+    void* /* udata */);
+
+/*
+ * the purpose of lmstr1 is to minimize the sum of the squares of
+ * m nonlinear functions in n variables by a modification of
+ * the levenberg-marquardt algorithm which uses minimal storage.
+ * this is done by using the more general least-squares solver
+ * lmstr. the user must provide a subroutine which calculates
+ * the functions and the rows of the jacobian.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_lmstr1(
+    minpack_fcn_lmstr /* fcn */,
+    int /* m */,
+    int /* n */,
+    double* /* x */,
+    double* /* fvec */,
+    double* /* fjac */,
+    int /* ldfjac */,
+    double /* tol */,
+    int* /* info */,
+    int* /* ipvt */,
+    double* /* wa */,
+    int /* lwa */,
+    void* /* udata */);
+
+/*
+ * this subroutine checks the gradients of m nonlinear functions
+ * in n variables, evaluated at a point x, for consistency with
+ * the functions themselves.
+ *
+ * the subroutine does not perform reliably if cancellation or
+ * rounding errors cause a severe loss of significance in the
+ * evaluation of a function. therefore, none of the components
+ * of x should be unusually small (in particular, zero) or any
+ * other value which may cause loss of significance.
+ */
+MINPACK_EXTERN void MINPACK_CALL
+minpack_chkder(
+    int /* m */,
+    int /* n */,
+    const double* /* x */,
+    const double* /* fvec */,
+    const double* /* fjac */,
+    int /* ldfjac */,
+    double* /* xp */,
+    const double* /* fvecp */,
+    int /* mode */,
+    double* /* err */);
diff --git a/src/minpack.f90 b/src/minpack.f90
index 7cf47d3..562d886 100644
--- a/src/minpack.f90
+++ b/src/minpack.f90
@@ -97,9 +97,9 @@ subroutine fcn_lmstr(m, n, x, fvec, fjrow, iflag)
                                        !! the value of iflag should not be changed by fcn unless
                                        !! the user wants to terminate execution of lmstr.
                                        !! in this case set iflag to a negative integer.
-            real(wp) :: x(n) !! independant variable vector
-            real(wp) :: fvec(m) !! value of function at `x`
-            real(wp) :: fjrow(n) !! jacobian row
+            real(wp), intent(in) :: x(n) !! independent variable vector
+            real(wp), intent(inout) :: fvec(m) !! value of function at `x`
+            real(wp), intent(inout) :: fjrow(n) !! jacobian row
         end subroutine fcn_lmstr
 
     end interface
@@ -2725,10 +2725,10 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, &
                                           !! multiplicative scale factors for the variables.
         real(wp), intent(out) :: Qtf(n) !! an output array of length n which contains
                                        !! the first n elements of the vector (q transpose)*fvec.
-        real(wp) :: Wa1(n) !! work array of length n.
-        real(wp) :: Wa2(n) !! work array of length n.
-        real(wp) :: Wa3(n) !! work array of length n.
-        real(wp) :: Wa4(m) !! work array of length m.
+        real(wp), intent(inout) :: Wa1(n) !! work array of length n.
+        real(wp), intent(inout) :: Wa2(n) !! work array of length n.
+        real(wp), intent(inout) :: Wa3(n) !! work array of length n.
+        real(wp), intent(inout) :: Wa4(m) !! work array of length m.
 
         integer :: i, iflag, iter, j, l
         real(wp) :: actred, delta, dirder, fnorm, &
@@ -3808,4 +3808,4 @@ end subroutine rwupdt
 
 !*****************************************************************************************
 end module minpack_module
-!*****************************************************************************************
\ No newline at end of file
+!*****************************************************************************************
diff --git a/src/minpack_capi.f90 b/src/minpack_capi.f90
new file mode 100644
index 0000000..cd91f90
--- /dev/null
+++ b/src/minpack_capi.f90
@@ -0,0 +1,472 @@
+module minpack_capi
+    use, intrinsic :: iso_c_binding, only : c_int, c_double, c_funptr, c_ptr, c_f_procpointer
+    use minpack_module
+    implicit none
+    private
+
+    public :: minpack_hybrd,minpack_hybrd1, minpack_hybrj, minpack_hybrj1, &
+        & minpack_lmdif, minpack_lmdif1, minpack_lmder, minpack_lmder1, &
+        & minpack_chkder
+
+    public :: minpack_dpmpar
+
+    public :: minpack_func, minpack_func2, minpack_fcn_hybrj, minpack_fcn_lmder, &
+        & minpack_fcn_lmstr
+
+    abstract interface
+        subroutine minpack_func(n, x, fvec, iflag, udata) bind(c)
+            import :: c_int, c_double, c_ptr
+            implicit none
+            integer(c_int), value :: n
+            real(c_double), intent(in) :: x(n)
+            real(c_double), intent(out) :: fvec(n)
+            integer(c_int), intent(inout) :: iflag
+            type(c_ptr), value :: udata
+        end subroutine minpack_func
+
+        subroutine minpack_func2(m, n, x, fvec, iflag, udata) bind(c)
+            import :: c_int, c_double, c_ptr
+            implicit none
+            integer(c_int), value :: m
+            integer(c_int), value :: n
+            real(c_double), intent(in) :: x(n)
+            real(c_double), intent(out) :: fvec(m)
+            integer(c_int), intent(inout) :: iflag
+            type(c_ptr), value :: udata
+        end subroutine minpack_func2
+
+        subroutine minpack_fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag, udata) bind(c)
+            import :: c_int, c_double, c_ptr
+            implicit none
+            integer(c_int), value :: n
+            real(c_double), intent(in) :: x(n)
+            integer(c_int), value :: ldfjac
+            real(c_double), intent(out) :: fvec(n)
+            real(c_double), intent(out) :: fjac(ldfjac, n)
+            integer(c_int), intent(inout) :: iflag
+            type(c_ptr), value :: udata
+        end subroutine minpack_fcn_hybrj
+
+        subroutine minpack_fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag, udata) bind(c)
+            import :: c_int, c_double, c_ptr
+            implicit none
+            integer(c_int), value :: m
+            integer(c_int), value :: n
+            integer(c_int), value :: ldfjac
+            integer(c_int), intent(inout) :: iflag
+            real(c_double), intent(in) :: x(n)
+            real(c_double), intent(inout) :: fvec(m)
+            real(c_double), intent(inout) :: fjac(ldfjac, n)
+            type(c_ptr), value :: udata
+        end subroutine minpack_fcn_lmder
+
+        subroutine minpack_fcn_lmstr(m, n, x, fvec, fjrow, iflag, udata) bind(c)
+            import :: c_int, c_double, c_ptr
+            implicit none
+            integer(c_int), value :: m
+            integer(c_int), value :: n
+            integer(c_int), intent(inout) :: iflag
+            real(c_double), intent(in) :: x(n)
+            real(c_double), intent(inout) :: fvec(m)
+            real(c_double), intent(inout) :: fjrow(n)
+            type(c_ptr), value :: udata
+        end subroutine minpack_fcn_lmstr
+    end interface
+
+contains
+
+    function minpack_dpmpar(i) result(par) bind(c)
+        integer(c_int), value :: i
+        real(c_double) :: par
+        if (i > 0_c_int .and. i <= 3_c_int) then
+            par = dpmpar(i)
+        end if
+    end function minpack_dpmpar
+
+    subroutine minpack_hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
+            & factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4, &
+            & udata) &
+            & bind(c)
+        procedure(minpack_func) :: fcn
+        integer(c_int), value :: n
+        integer(c_int), value :: maxfev
+        integer(c_int), value :: ml
+        integer(c_int), value :: mu
+        integer(c_int), value :: mode
+        integer(c_int), value :: nprint
+        integer(c_int), intent(out) :: info
+        integer(c_int), intent(out) :: nfev
+        integer(c_int), value :: ldfjac
+        integer(c_int), value :: lr
+        real(c_double), value :: xtol
+        real(c_double), value :: epsfcn
+        real(c_double), value :: factor
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(n)
+        real(c_double), intent(inout) :: diag(n)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(out) :: r(lr)
+        real(c_double), intent(out) :: qtf(n)
+        real(c_double), intent(inout) :: wa1(n)
+        real(c_double), intent(inout) :: wa2(n)
+        real(c_double), intent(inout) :: wa3(n)
+        real(c_double), intent(inout) :: wa4(n)
+        type(c_ptr), value :: udata
+
+        call hybrd(wrap_fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, &
+            & factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4)
+
+    contains
+        subroutine wrap_fcn(n, x, fvec, iflag)
+            integer, intent(in) :: n
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(out) :: fvec(n)
+            integer, intent(inout) :: iflag
+
+            call fcn(n, x, fvec, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_hybrd
+
+    subroutine minpack_hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa, udata) &
+          & bind(c)
+        procedure(minpack_func) :: fcn
+        integer(c_int), value :: n
+        integer(c_int), intent(out) :: info
+        real(c_double), value :: tol
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(n)
+        integer(c_int), value :: Lwa
+        real(c_double), intent(inout) :: Wa(Lwa)
+        type(c_ptr), value :: udata
+
+        call hybrd1(wrap_fcn, n, x, fvec, tol, info, Wa, Lwa)
+
+    contains
+        subroutine wrap_fcn(n, x, fvec, iflag)
+            integer, intent(in) :: n
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(out) :: fvec(n)
+            integer, intent(inout) :: iflag
+
+            call fcn(n, x, fvec, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_hybrd1
+
+    subroutine minpack_hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, &
+            & factor, nprint, info, nfev, njev, r, lr, qtf, wa1, wa2, wa3, wa4, udata) &
+            & bind(c)
+        procedure(minpack_fcn_hybrj) :: fcn
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), value :: maxfev
+        integer(c_int), value :: mode
+        integer(c_int), value :: nprint
+        integer(c_int), intent(out) :: info
+        integer(c_int), intent(out) :: nfev
+        integer(c_int), intent(out) :: njev
+        integer(c_int), value :: lr
+        real(c_double), value :: xtol
+        real(c_double), value :: factor
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(n)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(inout) :: diag(n)
+        real(c_double), intent(out) :: r(lr)
+        real(c_double), intent(out) :: qtf(n)
+        real(c_double), intent(inout) :: wa1(n)
+        real(c_double), intent(inout) :: wa2(n)
+        real(c_double), intent(inout) :: wa3(n)
+        real(c_double), intent(inout) :: wa4(n)
+        type(c_ptr), value :: udata
+
+        call hybrj(wrap_fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, &
+            & factor, nprint, info, nfev, njev, r, lr, qtf, wa1, wa2, wa3, wa4)
+
+    contains
+        subroutine wrap_fcn(n, x, fvec, fjac, ldfjac, iflag)
+            integer, intent(in) :: n
+            real(wp), intent(in) :: x(n)
+            integer, intent(in) :: ldfjac
+            real(wp), intent(out) :: fvec(n)
+            real(wp), intent(out) :: fjac(ldfjac, n)
+            integer, intent(inout) :: iflag
+
+            call fcn(n, x, fvec, fjac, ldfjac, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_hybrj
+
+    subroutine minpack_hybrj1(fcn, n, x, fvec, fjac, ldfjac, tol, info, wa, lwa, udata) &
+            & bind(c)
+        procedure(minpack_fcn_hybrj) :: fcn
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), intent(out) :: info
+        integer(c_int), value :: lwa
+        real(c_double), value :: tol
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(n)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(inout) :: wa(lwa)
+        type(c_ptr), value :: udata
+
+        call hybrj1(wrap_fcn, n, x, fvec, fjac, ldfjac, tol, info, wa, lwa)
+
+    contains
+        subroutine wrap_fcn(n, x, fvec, fjac, ldfjac, iflag)
+            integer, intent(in) :: n
+            real(wp), intent(in) :: x(n)
+            integer, intent(in) :: ldfjac
+            real(wp), intent(out) :: fvec(n)
+            real(wp), intent(out) :: fjac(ldfjac, n)
+            integer, intent(inout) :: iflag
+
+            call fcn(n, x, fvec, fjac, ldfjac, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_hybrj1
+
+    subroutine minpack_lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, &
+            & mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4, &
+            & udata) &
+            & bind(c)
+        procedure(minpack_func2) :: fcn
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), value :: maxfev
+        integer(c_int), value :: mode
+        integer(c_int), value :: nprint
+        integer(c_int), intent(out) :: info
+        integer(c_int), intent(out) :: nfev
+        integer(c_int), value :: ldfjac
+        integer(c_int), intent(out) :: ipvt(n)
+        real(c_double), value :: ftol
+        real(c_double), value :: xtol
+        real(c_double), value :: gtol
+        real(c_double), value :: epsfcn
+        real(c_double), value :: factor
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(m)
+        real(c_double), intent(inout) :: diag(n)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(out) :: qtf(n)
+        real(c_double), intent(inout) :: wa1(n)
+        real(c_double), intent(inout) :: wa2(n)
+        real(c_double), intent(inout) :: wa3(n)
+        real(c_double), intent(inout) :: wa4(m)
+        type(c_ptr), value :: udata
+
+        call lmdif(wrap_fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, &
+            & mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4)
+
+    contains
+        subroutine wrap_fcn(m, n, x, fvec, iflag)
+            integer, intent(in) :: m
+            integer, intent(in) :: n
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(out) :: fvec(m)
+            integer, intent(inout) :: iflag
+
+            call fcn(m, n, x, fvec, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_lmdif
+
+    subroutine minpack_lmdif1(fcn, m, n, x, fvec, tol, info, iwa, wa, lwa, udata) &
+            & bind(c)
+        procedure(minpack_func2) :: fcn
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), intent(out) :: info
+        integer(c_int), value :: lwa
+        integer(c_int), intent(inout) :: iwa(n)
+        real(c_double), value :: tol
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(inout) :: fvec(m)
+        real(c_double), intent(inout) :: wa(lwa)
+        type(c_ptr), value :: udata
+
+        call lmdif1(wrap_fcn, m, n, x, fvec, tol, info, iwa, wa, lwa)
+
+    contains
+        subroutine wrap_fcn(m, n, x, fvec, iflag)
+            integer, intent(in) :: m
+            integer, intent(in) :: n
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(out) :: fvec(m)
+            integer, intent(inout) :: iflag
+
+            call fcn(m, n, x, fvec, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_lmdif1
+
+    subroutine minpack_lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
+            & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4, &
+            & udata) &
+            & bind(c)
+        procedure(minpack_fcn_lmder) :: fcn
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), value :: maxfev
+        integer(c_int), value :: mode
+        integer(c_int), value :: nprint
+        integer(c_int), intent(out) :: info
+        integer(c_int), intent(out) :: nfev
+        integer(c_int), intent(out) :: njev
+        integer(c_int), intent(out) :: ipvt(n)
+        real(c_double), value :: ftol
+        real(c_double), value :: xtol
+        real(c_double), value :: gtol
+        real(c_double), value :: factor
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(m)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(inout) :: diag(n)
+        real(c_double), intent(out) :: qtf(n)
+        real(c_double), intent(inout) :: wa1(n)
+        real(c_double), intent(inout) :: wa2(n)
+        real(c_double), intent(inout) :: wa3(n)
+        real(c_double), intent(inout) :: wa4(m)
+        type(c_ptr), value :: udata
+
+        call lmder(wrap_fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
+            & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4)
+
+    contains
+        subroutine wrap_fcn(m, n, x, fvec, fjac, ldfjac, iflag)
+            integer, intent(in) :: m
+            integer, intent(in) :: n
+            integer, intent(in) :: ldfjac
+            integer, intent(inout) :: iflag
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(inout) :: fvec(m)
+            real(wp), intent(inout) :: fjac(ldfjac, n)
+
+            call fcn(m, n, x, fvec, fjac, ldfjac, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_lmder
+
+    subroutine minpack_lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, &
+            & udata) &
+            & bind(c)
+        procedure(minpack_fcn_lmder) :: fcn
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), intent(out) :: info
+        integer(c_int), value :: lwa
+        integer(c_int), intent(out) :: ipvt(n)
+        real(c_double), value :: tol
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(m)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(inout) :: wa(lwa)
+        type(c_ptr), value :: udata
+
+        call lmder1(wrap_fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa)
+
+    contains
+        subroutine wrap_fcn(m, n, x, fvec, fjac, ldfjac, iflag)
+            integer, intent(in) :: m
+            integer, intent(in) :: n
+            integer, intent(in) :: ldfjac
+            integer, intent(inout) :: iflag
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(inout) :: fvec(m)
+            real(wp), intent(inout) :: fjac(ldfjac, n)
+
+            call fcn(m, n, x, fvec, fjac, ldfjac, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_lmder1
+
+    subroutine minpack_lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
+            & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4, &
+            & udata) &
+            & bind(c)
+        procedure(minpack_fcn_lmstr) :: fcn
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), value :: maxfev
+        integer(c_int), value :: mode
+        integer(c_int), value :: nprint
+        integer(c_int), intent(out) :: info
+        integer(c_int), intent(out) :: nfev
+        integer(c_int), intent(out) :: njev
+        integer(c_int), intent(out) :: ipvt(n)
+        real(c_double), value :: ftol
+        real(c_double), value :: xtol
+        real(c_double), value :: gtol
+        real(c_double), value :: factor
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(m)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(inout) :: diag(n)
+        real(c_double), intent(out) :: qtf(n)
+        real(c_double), intent(inout) :: wa1(n)
+        real(c_double), intent(inout) :: wa2(n)
+        real(c_double), intent(inout) :: wa3(n)
+        real(c_double), intent(inout) :: wa4(m)
+        type(c_ptr), value :: udata
+
+        call lmstr(wrap_fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, &
+            & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4)
+    contains
+        subroutine wrap_fcn(m, n, x, fvec, fjrow, iflag)
+            integer, intent(in) :: m
+            integer, intent(in) :: n
+            integer, intent(inout) :: iflag
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(inout) :: fvec(m)
+            real(wp), intent(inout) :: fjrow(n)
+
+            call fcn(m, n, x, fvec, fjrow, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_lmstr
+
+    subroutine minpack_lmstr1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info, ipvt, wa, lwa, &
+            & udata) &
+            & bind(c)
+        procedure(minpack_fcn_lmstr) :: fcn
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), intent(out) :: info
+        integer(c_int), value :: lwa
+        integer(c_int), intent(out) :: ipvt(n)
+        real(c_double), value :: tol
+        real(c_double), intent(inout) :: x(n)
+        real(c_double), intent(out) :: fvec(m)
+        real(c_double), intent(out) :: fjac(ldfjac, n)
+        real(c_double), intent(inout) :: wa(lwa)
+        type(c_ptr), value :: udata
+
+        call lmstr1(wrap_fcn, m, n, x, fvec, fjac, ldfjac, tol, info, ipvt, wa, lwa)
+    contains
+        subroutine wrap_fcn(m, n, x, fvec, fjrow, iflag)
+            integer, intent(in) :: m
+            integer, intent(in) :: n
+            integer, intent(inout) :: iflag
+            real(wp), intent(in) :: x(n)
+            real(wp), intent(inout) :: fvec(m)
+            real(wp), intent(inout) :: fjrow(n)
+
+            call fcn(m, n, x, fvec, fjrow, iflag, udata)
+        end subroutine wrap_fcn
+    end subroutine minpack_lmstr1
+
+    subroutine minpack_chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err) &
+         & bind(c)
+        integer(c_int), value :: m
+        integer(c_int), value :: n
+        integer(c_int), value :: ldfjac
+        integer(c_int), value :: mode
+        real(c_double), intent(in) :: x(n)
+        real(c_double), intent(in) :: fvec(m)
+        real(c_double), intent(in) :: fjac(ldfjac, n)
+        real(c_double), intent(out) :: xp(n)
+        real(c_double), intent(in) :: fvecp(m)
+        real(c_double), intent(out) :: err(m)
+
+        call chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err)
+    end subroutine minpack_chkder
+
+end module minpack_capi
diff --git a/test/api/tester.c b/test/api/tester.c
new file mode 100644
index 0000000..03fb536
--- /dev/null
+++ b/test/api/tester.c
@@ -0,0 +1,414 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdbool.h>
+#include <math.h>
+#include <assert.h>
+
+#include "testsuite.h"
+#include "minpack.h"
+
+static double
+enorm(const int n, const double* x)
+{
+    double norm = 0.0;
+    for(int i = 0; i < n; i++) {
+        norm += x[i]*x[i];
+    }
+    return sqrt(norm);
+}
+
+void
+trial_hybrd_fcn(int n, const double* x, double* fvec, int* iflag, void* udata) {
+    assert(!udata);  // we always pass NULL in the test
+    if (*iflag == 0) return;
+    for(int k = 0; k < n; k++) {
+        double temp = (3.0 - 2.0*x[k])*x[k];
+        double temp1 = k != 0 ? x[k - 1] : 0.0;
+        double temp2 = k != n-1 ? x[k + 1] : 0.0;
+        fvec[k] = temp - temp1 - 2.0*temp2 + 1.0;
+    }
+}
+
+int
+test_hybrd1 (void)
+{
+    int n = 9;
+    int info = 0;
+    int lwa = 180;
+    double x[9] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
+    double fvec[9];
+    double wa[180];
+    double tol = sqrt(minpack_dpmpar(1));
+    double reference[9] = {
+        -0.5706545, -0.6816283, -0.7017325,
+        -0.7042129, -0.7013690, -0.6918656,
+        -0.6657920, -0.5960342, -0.4164121};
+
+    minpack_hybrd1(trial_hybrd_fcn, n, x, fvec, tol, &info, wa, lwa, NULL);
+
+    if (!check(info, 1, "Unexpected info value")) return 1;
+
+    if (!check(enorm(n, fvec), 0.0, tol, "Unexpected residual")) return 1;
+
+    for(int i = 0; i < 9; i++) {
+        if (!check(x[i], reference[i], 10*tol, "Unexpected solution")) return 1;
+    }
+
+    return 0;
+}
+
+int
+test_hybrd (void)
+{
+    int n = 9;
+    int info = 0, nfev = 0;
+    double x[9] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
+    double diag[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
+    double fvec[9];
+    double fjac[9*9];
+    double r[45], qtf[9], wa1[9], wa2[9], wa3[9], wa4[9];
+    double tol = sqrt(minpack_dpmpar(1));
+    double reference[9] = {
+        -0.5706545, -0.6816283, -0.7017325,
+        -0.7042129, -0.7013690, -0.6918656,
+        -0.6657920, -0.5960342, -0.4164121};
+
+    minpack_hybrd(trial_hybrd_fcn, n, x, fvec, tol, 2000, 1, 1, 0.0, diag, 2, 100.0, 0,
+            &info, &nfev, fjac, 9, r, 45, qtf, wa1, wa2, wa3, wa4, NULL);
+
+    if (!check(info, 1, "Unexpected info value")) return 1;
+    if (!check(nfev, 14, "Unexpected number of function calls")) return 1;
+
+    if (!check(enorm(n, fvec), 0.0, tol, "Unexpected residual")) return 1;
+
+    for(int i = 0; i < 9; i++) {
+        if (!check(x[i], reference[i], 10*tol, "Unexpected solution")) return 1;
+    }
+
+    return 0;
+}
+
+void
+trial_hybrj_fcn(int n, const double* x, double* fvec, double* fjac, int ldfjac,
+        int* iflag, void* udata) {
+    if (*iflag == 1) {
+        trial_hybrd_fcn(n, x, fvec, iflag, udata);
+    } else {
+        for(int k = 0; k < n; k++) {
+            for(int j = 0; j < n; j++) {
+                fjac[k*ldfjac + j] = 0.0;
+            }
+            fjac[k*ldfjac + k] = 3.0 - 4.0*x[k];
+            if (k != 0) fjac[k*ldfjac + k - 1] = -1.0;
+            if (k != n-1) fjac[k*ldfjac + k + 1] = -2.0;
+        }
+        
+    }
+}
+
+int
+test_hybrj1 (void)
+{
+    int n = 9;
+    int info = 0;
+    int lwa = 180;
+    double x[9] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
+    double fvec[9], fjac[9*9];
+    double wa[180];
+    double tol = sqrt(minpack_dpmpar(1));
+    double reference[9] = {
+        -0.5706545, -0.6816283, -0.7017325,
+        -0.7042129, -0.7013690, -0.6918656,
+        -0.6657920, -0.5960342, -0.4164121};
+
+    minpack_hybrj1(trial_hybrj_fcn, n, x, fvec, fjac, n, tol, &info, wa, lwa, NULL);
+
+    if (!check(info, 1, "Unexpected info value")) return 1;
+
+    if (!check(enorm(n, fvec), 0.0, tol, "Unexpected residual")) return 1;
+
+    for(int i = 0; i < 9; i++) {
+        if (!check(x[i], reference[i], 10*tol, "Unexpected solution")) return 1;
+    }
+
+    return 0;
+}
+
+int
+test_hybrj (void)
+{
+    int n = 9;
+    int info = 0, nfev = 0, njev = 0;
+    double x[9] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0};
+    double diag[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
+    double fvec[9], fjac[9*9];
+    double r[45], qtf[9], wa1[9], wa2[9], wa3[9], wa4[9];
+    double tol = sqrt(minpack_dpmpar(1));
+    double reference[9] = {
+        -0.5706545, -0.6816283, -0.7017325,
+        -0.7042129, -0.7013690, -0.6918656,
+        -0.6657920, -0.5960342, -0.4164121};
+
+    minpack_hybrj(trial_hybrj_fcn, n, x, fvec, fjac, n, tol, 2000, diag, 2, 100.0, 0,
+            &info, &nfev, &njev, r, 45, qtf, wa1, wa2, wa3, wa4, NULL);
+
+    if (!check(info, 1, "Unexpected info value")) return 1;
+    if (!check(nfev, 15, "Unexpected number of function evaluations")) return 1;
+    if (!check(njev, 1, "Unexpected number of jacobian evaluations")) return 1;
+
+    if (!check(enorm(n, fvec), 0.0, tol, "Unexpected residual")) return 1;
+
+    for(int i = 0; i < 9; i++) {
+        if (!check(x[i], reference[i], 10*tol, "Unexpected solution")) return 1;
+    }
+
+    return 0;
+}
+
+void
+trial_lmder_fcn(int m, int n, const double* x, double* fvec, double* fjac,
+                int ldfjac, int* iflag, void* data) {
+    assert(!!data);
+    double* y = (double*)data;
+    assert(m == 15);
+    assert(n == 3);
+
+    if (*iflag == 1) {
+        for (int i = 0; i < m; i++) {
+            double tmp1 = i + 1;
+            double tmp2 = 16 - i - 1;
+            double tmp3 = i >= 8 ? tmp2 : tmp1;
+            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
+        }
+    } else if (*iflag == 2) {
+        assert(ldfjac == m);
+        for (int i = 0; i < m; i++) {
+            double tmp1 = i + 1;
+            double tmp2 = 16 - i - 1;
+            double tmp3 = i >= 8 ? tmp2 : tmp1;
+            double tmp4 = pow(x[1]*tmp2 + x[2]*tmp3, 2);
+            fjac[i] = -1.0;
+            fjac[i + ldfjac] = tmp1*tmp2/tmp4;
+            fjac[i + 2*ldfjac] = tmp1*tmp3/tmp4;
+        }
+    }
+}
+
+int
+test_lmder1 (void)
+{
+    const double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1,
+        3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34e0, 2.1e0, 4.39e0};
+    const int m = 15, n = 3;
+    int info = 0;
+    double x[3] = {1.0, 1.0, 1.0}, xp[n];
+    double fvec[m], fvecp[m], err[m];
+    double fjac[m*n];
+    int ipvt[n];
+    double wa[5*n+m];
+    double tol = sqrt(minpack_dpmpar(1));
+
+    minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err);
+    info = 1;
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    info = 2;
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    info = 1;
+    trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, y);
+    minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err);
+
+    for (int i = 0; i < 15; i++) {
+        if (!check(err[i], 1.0, tol, "Unexpected derivatives")) return 1;
+    }
+
+    minpack_lmder1(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, &info, ipvt, wa, 30, y);
+    if (!check(info, 1, "Unexpected info value")) return 1;
+    if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
+    if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
+    if (!check(x[2], 0.2343695e+1, 100*tol, "Unexpected x[2]")) return 1;
+    if (!check(enorm(m, fvec), 0.9063596e-1, tol, "Unexpected residual")) return 1;
+
+    return 0;
+}
+
+int
+test_lmder (void)
+{
+    const double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1,
+        3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34e0, 2.1e0, 4.39e0};
+    const int m = 15, n = 3;
+    int info = 0, nfev = 0, njev = 0;
+    double x[3] = {1.0, 1.0, 1.0}, xp[n];
+    double fvec[m], fvecp[m], err[m];
+    double fjac[m*n];
+    int ipvt[n];
+    double diag[n], qtf[n], wa1[n], wa2[n], wa3[n], wa4[m];
+    double tol = sqrt(minpack_dpmpar(1));
+
+    minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err);
+    info = 1;
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    info = 2;
+    trial_lmder_fcn(m, n, x, fvec, fjac, m, &info, y);
+    info = 1;
+    trial_lmder_fcn(m, n, xp, fvecp, fjac, m, &info, y);
+    minpack_chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err);
+
+    for (int i = 0; i < 15; i++) {
+        if (!check(err[i], 1.0, tol, "Unexpected derivatives")) return 1;
+    }
+
+    minpack_lmder(trial_lmder_fcn, m, n, x, fvec, fjac, m, tol, tol, 0.0, 2000, diag, 1,
+            100.0, 0, &info, &nfev, &njev, ipvt, qtf, wa1, wa2, wa3, wa4, y);
+    if (!check(info, 1, "Unexpected info value")) return 1;
+    if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
+    if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
+    if (!check(x[2], 0.2343695e+1, 100*tol, "Unexpected x[2]")) return 1;
+    if (!check(enorm(m, fvec), 0.9063596e-1, tol, "Unexpected residual")) return 1;
+
+    return 0;
+}
+
+void
+trial_lmdif_fcn(int m, int n, const double* x, double* fvec, int* iflag, void* data) {
+    assert(!!data);
+    double* y = (double*)data;
+    assert(m == 15);
+    assert(n == 3);
+    if (*iflag == 0) return;
+
+    for (int i = 0; i < m; i++) {
+        double tmp1 = i + 1;
+        double tmp2 = 16 - i - 1;
+        double tmp3 = i >= 8 ? tmp2 : tmp1;
+        fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
+    }
+}
+
+int
+test_lmdif1 (void)
+{
+    const int m = 15, n = 3;
+    double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1,
+        3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34e0, 2.1e0, 4.39e0};
+    double x[3] = {1.0, 1.0, 1.0}, fvec[15];
+    int info = 0;
+    double tol = sqrt(minpack_dpmpar(1));
+    int ipvt[n];
+    int lwa = m*n + 5*n + m;
+    double wa[lwa];
+
+    minpack_lmdif1(trial_lmdif_fcn, 15, 3, x, fvec, tol, &info, ipvt, wa, lwa, y);
+    if (!check(info, 1, "Unexpected info value")) return 1;
+    if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
+    if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
+    if (!check(x[2], 0.2343695e+1, 100*tol, "Unexpected x[2]")) return 1;
+    if (!check(enorm(m, fvec), 0.9063596e-1, tol, "Unexpected residual")) return 1;
+
+    return 0;
+}
+
+int
+test_lmdif (void)
+{
+    const int m = 15, n = 3;
+    double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1,
+        3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34e0, 2.1e0, 4.39e0};
+    double x[3] = {1.0, 1.0, 1.0}, fvec[15];
+    int info = 0, nfev = 0;
+    double tol = sqrt(minpack_dpmpar(1));
+    int ipvt[n];
+    double fjac[m*n], diag[n], qtf[n], wa1[n], wa2[n], wa3[n], wa4[m];
+
+    minpack_lmdif(trial_lmdif_fcn, 15, 3, x, fvec, tol, tol, 0.0, 2000, 0.0, diag, 1,
+            100.0, 0, &info, &nfev, fjac, 15, ipvt, qtf, wa1, wa2, wa3, wa4, y);
+    if (!check(info, 1, "Unexpected info value")) return 1;
+    if (!check(x[0], 0.8241058e-1, 100*tol, "Unexpected x[0]")) return 1;
+    if (!check(x[1], 0.1133037e+1, 100*tol, "Unexpected x[1]")) return 1;
+    if (!check(x[2], 0.2343695e+1, 100*tol, "Unexpected x[2]")) return 1;
+    if (!check(enorm(m, fvec), 0.9063596e-1, tol, "Unexpected residual")) return 1;
+
+    return 0;
+}
+
+void
+trial_lmstr_fcn(int m, int n, const double* x, double* fvec, double* fjrow, int* iflag,
+        void* udata) {
+    if (*iflag == 1) {
+        fvec[0] = 10.0 * (x[1] - x[0] * x[0]);
+        fvec[1] = 1.0 - x[0];
+    } else {
+        if (*iflag == 2) {
+            fjrow[0] = -20.0 * x[0];
+            fjrow[1] = 10.0;
+        } else {
+            fjrow[0] = -1.0;
+            fjrow[1] = 0.0;
+        }
+    }
+}
+
+int
+test_lmstr1 (void)
+{
+    const int m = 2, n = 2;
+    double x[2] = {-1.2, 1.0}, fvec[2], fjac[2*2];
+    int info = 0;
+    double tol = sqrt(minpack_dpmpar(1));
+    int ipvt[n];
+    int lwa = m*n + 5*n + m;
+    double wa[lwa];
+
+    minpack_lmstr1(trial_lmstr_fcn, 2, 2, x, fvec, fjac, 2, tol, &info, ipvt, wa, lwa, NULL);
+    if (!check(info, 4, "Unexpected info value")) return 1;
+    if (!check(x[0], 1.0, 100*tol, "Unexpected x[0]")) return 1;
+    if (!check(x[1], 1.0, 100*tol, "Unexpected x[1]")) return 1;
+    if (!check(enorm(m, fvec), 0.0, tol, "Unexpected residual")) return 1;
+
+    return 0;
+}
+
+int
+test_lmstr (void)
+{
+    const int m = 2, n = 2;
+    double x[2] = {-1.2, 1.0}, fvec[2], fjac[2*2], diag[2];
+    int info = 0, nfev = 0, njev = 0;
+    double tol = sqrt(minpack_dpmpar(1));
+    int ipvt[n];
+    double qtf[n], wa1[n], wa2[n], wa3[n], wa4[m];
+
+    minpack_lmstr(trial_lmstr_fcn, 2, 2, x, fvec, fjac, 2, tol, tol, 0.0, 2000, diag, 1,
+            100.0, 0, &info, &nfev, &njev, ipvt, qtf, wa1, wa2, wa3, wa4, NULL);
+    if (!check(info, 4, "Unexpected info value")) return 1;
+    if (!check(nfev, 21, "Unexpected number of function evaluations")) return 1;
+    if (!check(njev, 16, "Unexpected number of jacobian evaluations")) return 1;
+    if (!check(x[0], 1.0, 100*tol, "Unexpected x[0]")) return 1;
+    if (!check(x[1], 1.0, 100*tol, "Unexpected x[1]")) return 1;
+    if (!check(enorm(m, fvec), 0.0, tol, "Unexpected residual")) return 1;
+
+    return 0;
+}
+
+int
+main (void) {
+    int stat = 0;
+
+    stat += run("hybrd1", test_hybrd1);
+    stat += run("hybrd ", test_hybrd);
+    stat += run("hybrj1", test_hybrj1);
+    stat += run("hybrj ", test_hybrj);
+    stat += run("lmder1", test_lmder1);
+    stat += run("lmder ", test_lmder);
+    stat += run("lmdif1", test_lmdif1);
+    stat += run("lmdif ", test_lmdif);
+    stat += run("lmstr1", test_lmstr1);
+    stat += run("lmstr ", test_lmstr);
+
+    if (stat > 0) {
+        fprintf(stderr, "[FAIL] %d test(s) failed\n", stat);
+    } else {
+        fprintf(stderr, "[PASS] all tests passed\n");
+    }
+    return stat == 0 ? EXIT_SUCCESS : EXIT_FAILURE;
+}
diff --git a/test/api/testsuite.h b/test/api/testsuite.h
new file mode 100644
index 0000000..3fb90dc
--- /dev/null
+++ b/test/api/testsuite.h
@@ -0,0 +1,49 @@
+#pragma once
+
+#include <stdbool.h>
+#include <stdio.h>
+#include <math.h>
+
+#define check(x, ...) \
+    _Generic((x), \
+        int: check_int, \
+     double: check_double \
+            )(x, __VA_ARGS__)
+
+static inline bool
+check_int(int actual, int expected, const char *msg)
+{
+    if (expected == actual) {
+        return true;
+    }
+    fprintf(stderr, "FAIL: %s: expected %d, got %d\n", msg, expected, actual);
+    return false;
+}
+
+static inline bool
+check_double(double actual, double expected, double tol, const char *msg)
+{
+    if (fabs(expected - actual) < tol) {
+        return true;
+    }
+    fprintf(stderr, "FAIL: %s: expected %g, got %g\n", msg, expected, actual);
+    return false;
+}
+
+static inline bool
+is_close(double a, double b, double tol)
+{
+    return fabs(a - b) < tol;
+}
+
+static inline int
+run(const char* name, int (*test)(void)) {
+    printf("Testing %s ...", name);
+    int stat = test();
+    if (stat) {
+        printf(" FAILED\n");
+    } else {
+        printf(" OK\n");
+    }
+    return stat;
+}