diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c
index 89226d0d5a..8731465cd7 100644
--- a/c/examples/bobyqa/bobyqa_example.c
+++ b/c/examples/bobyqa/bobyqa_example.c
@@ -57,7 +57,7 @@ int main(int argc, char * argv[])
 
     // Call the solver
     prima_result_t result;
-    const int rc = prima_minimize(PRIMA_BOBYQA, problem, options, &result);
+    const int rc = prima_minimize(PRIMA_BOBYQA, &problem, &options, &result);
 
     // Print the result
     printf("x* = {%g, %g}, rc = %d, msg = '%s', evals = %d\n", result.x[0], result.x[1], rc, result.message, result.nf);
diff --git a/c/examples/cobyla/cobyla_example.c b/c/examples/cobyla/cobyla_example.c
index 0e4652d188..f3d6c6daf3 100644
--- a/c/examples/cobyla/cobyla_example.c
+++ b/c/examples/cobyla/cobyla_example.c
@@ -64,7 +64,7 @@ int main(int argc, char * argv[])
 
     // Call the solver
     prima_result_t result;
-    const int rc = prima_minimize(PRIMA_COBYLA, problem, options, &result);
+    const int rc = prima_minimize(PRIMA_COBYLA, &problem, &options, &result);
 
     // Print the result
     printf("x* = {%g, %g}, f* = %g, cstrv = %g, nlconstr = {%g}, rc = %d, msg = '%s', evals = %d\n", result.x[0], result.x[1], result.f, result.cstrv, result.nlconstr[0], rc, result.message, result.nf);
diff --git a/c/examples/lincoa/lincoa_example.c b/c/examples/lincoa/lincoa_example.c
index 1018c71c27..b46f402e09 100644
--- a/c/examples/lincoa/lincoa_example.c
+++ b/c/examples/lincoa/lincoa_example.c
@@ -58,7 +58,7 @@ int main(int argc, char * argv[])
 
     // Call the solver
     prima_result_t result;
-    const int rc = prima_minimize(PRIMA_LINCOA, problem, options, &result);
+    const int rc = prima_minimize(PRIMA_LINCOA, &problem, &options, &result);
 
     // Print the result
     printf("x* = {%g, %g}, f* = %g, cstrv = %g, rc = %d, msg = '%s', evals = %d\n", result.x[0], result.x[1], result.f, result.cstrv, rc, result.message, result.nf);
diff --git a/c/examples/newuoa/newuoa_example.c b/c/examples/newuoa/newuoa_example.c
index 86d4c90bcb..f7544c3d88 100644
--- a/c/examples/newuoa/newuoa_example.c
+++ b/c/examples/newuoa/newuoa_example.c
@@ -52,7 +52,7 @@ int main(int argc, char * argv[])
 
     // Call the solver
     prima_result_t result;
-    const int rc = prima_minimize(PRIMA_NEWUOA, problem, options, &result);
+    const int rc = prima_minimize(PRIMA_NEWUOA, &problem, &options, &result);
 
     // Print the result
     printf("x* = {%g, %g}, rc = %d, msg = '%s', evals = %d\n", result.x[0], result.x[1], rc, result.message, result.nf);
diff --git a/c/examples/uobyqa/uobyqa_example.c b/c/examples/uobyqa/uobyqa_example.c
index 220eab8a79..c5f61e5ae3 100644
--- a/c/examples/uobyqa/uobyqa_example.c
+++ b/c/examples/uobyqa/uobyqa_example.c
@@ -52,7 +52,7 @@ int main(int argc, char * argv[])
 
     // Run the solver
     prima_result_t result;
-    const int rc = prima_minimize(PRIMA_UOBYQA, problem, options, &result);
+    const int rc = prima_minimize(PRIMA_UOBYQA, &problem, &options, &result);
 
     // Print the result
     printf("x* = {%g, %g}, rc = %d, msg = '%s', evals = %d\n", result.x[0], result.x[1], rc, result.message, result.nf);
diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h
index 6a041816d8..b260b958bf 100644
--- a/c/include/prima/prima.h
+++ b/c/include/prima/prima.h
@@ -185,7 +185,6 @@ typedef struct {
 PRIMAC_API
 int prima_init_problem(prima_problem_t *const problem, const int n);
 
-
 // Structure to hold the options
 // In the following, "Default" refers to the value set by `prima_init_options`.
 typedef struct {
@@ -291,8 +290,125 @@ int prima_free_result(prima_result_t *const result);
  * return    : see prima_rc_t enum for return codes
  */
 PRIMAC_API
-int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t problem, const prima_options_t options, prima_result_t *const result);
+int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t *problem, const prima_options_t *options, prima_result_t *const result);
+
+
+typedef struct {
+    prima_algorithm_t algorithm;
+    prima_problem_t problem;
+    prima_options_t options;
+    prima_result_t result;
+} prima_context_t;
+
+PRIMAC_API
+prima_context_t* prima_context_create(const int n);
+
+PRIMAC_API
+void prima_context_destroy(prima_context_t* context);
+
+PRIMAC_API
+int prima_context_solve(prima_algorithm_t algorithm, prima_context_t* context);
+
+/* Accessors. */
+
+PRIMAC_API
+int prima_context_set_callfun(prima_context_t* context, prima_obj_t calfun, void* data);
+
+PRIMAC_API
+int prima_context_set_calcfc(prima_context_t* context, prima_objcon_t calcfc, void* data, int m_nlcon);
+
+PRIMAC_API
+int prima_context_set_x0(prima_context_t* context, double *x0);
+
+PRIMAC_API
+int prima_context_set_bounds(prima_context_t* context, double *xl, double *xu);
+
+PRIMAC_API
+int prima_context_set_ineq(prima_context_t* context, int m_ineq, double Aineq[], double *bineq);
+
+PRIMAC_API
+int prima_context_set_eq(prima_context_t* context, int m_eq, double Aeq[], double *beq);
+
+PRIMAC_API
+int prima_context_set_f0(prima_context_t* context, double f0);
+
+PRIMAC_API
+int prima_context_set_nlconstr0(prima_context_t* context, double* nlconstr0);
+
+PRIMAC_API
+int prima_context_get_n(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_get_m_ineq(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_get_m_eq(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_get_m_nlcon(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_set_rho(prima_context_t* context, double rhobeg, double rhoend);
+
+PRIMAC_API
+double prima_context_get_rhobeg(const prima_context_t* context);
+
+PRIMAC_API
+double prima_context_get_rhoend(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_set_maxfun(prima_context_t* context, int maxfun);
+
+PRIMAC_API
+int prima_context_get_maxfun(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_set_iprint(prima_context_t* context, int iprint);
+
+PRIMAC_API
+int prima_context_get_iprint(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_set_npt(prima_context_t* context, int npt);
+
+PRIMAC_API
+int prima_context_get_npt(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_set_ftarget(prima_context_t* context, double ftarget);
+
+PRIMAC_API
+double prima_context_get_ftarget(const prima_context_t* context);
 
+PRIMAC_API
+int prima_context_set_ctol(prima_context_t* context, double ctol);
+
+PRIMAC_API
+double prima_context_get_ctol(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_set_callback(prima_context_t* context, prima_callback_t callback);
+
+PRIMAC_API
+double* prima_context_get_x(const prima_context_t* context);
+
+PRIMAC_API
+double prima_context_get_f(const prima_context_t* context);
+
+PRIMAC_API
+double prima_context_get_cstrv(const prima_context_t* context);
+
+PRIMAC_API
+double* prima_context_get_nlconstr(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_get_nf(const prima_context_t* context);
+
+PRIMAC_API
+int prima_context_get_status(const prima_context_t* context);
+
+PRIMAC_API
+const char* prima_context_get_message(const prima_context_t* context);
 
 #ifdef __cplusplus
 }
diff --git a/c/prima.c b/c/prima.c
index ce16e276da..45747855a2 100644
--- a/c/prima.c
+++ b/c/prima.c
@@ -63,22 +63,25 @@ int prima_init_options(prima_options_t *const options)
 
 
 // Function to check whether the problem matches the algorithm
-int prima_check_problem(const prima_problem_t problem, const int use_constr, const prima_algorithm_t algorithm)
+int prima_check_problem(const prima_problem_t *problem, const int use_constr, const prima_algorithm_t algorithm)
 {
-    if (algorithm != PRIMA_COBYLA && (problem.calcfc || problem.nlconstr0 || problem.m_nlcon > 0))
+    if (problem == NULL)
+        return PRIMA_NULL_PROBLEM;
+
+    if (algorithm != PRIMA_COBYLA && (problem->calcfc || problem->nlconstr0 || problem->m_nlcon > 0))
         return PRIMA_PROBLEM_SOLVER_MISMATCH_NONLINEAR_CONSTRAINTS;
 
     if ((algorithm != PRIMA_COBYLA && algorithm != PRIMA_LINCOA) &&
-        (problem.m_ineq > 0 || problem.m_eq > 0 || problem.Aineq || problem.bineq || problem.Aeq || problem.beq))
+        (problem->m_ineq > 0 || problem->m_eq > 0 || problem->Aineq || problem->bineq || problem->Aeq || problem->beq))
         return PRIMA_PROBLEM_SOLVER_MISMATCH_LINEAR_CONSTRAINTS;
 
-    if ((algorithm != PRIMA_COBYLA && algorithm != PRIMA_LINCOA && algorithm != PRIMA_BOBYQA) && (problem.xl || problem.xu))
+    if ((algorithm != PRIMA_COBYLA && algorithm != PRIMA_LINCOA && algorithm != PRIMA_BOBYQA) && (problem->xl || problem->xu))
         return PRIMA_PROBLEM_SOLVER_MISMATCH_BOUNDS;
 
-    if (!problem.x0)
+    if (!problem->x0)
         return PRIMA_NULL_X0;
 
-    if ((use_constr && !problem.calcfc) || (!use_constr && !problem.calfun))
+    if ((use_constr && !problem->calcfc) || (!use_constr && !problem->calfun))
         return PRIMA_NULL_FUNCTION;
 
     return 0;
@@ -86,7 +89,7 @@ int prima_check_problem(const prima_problem_t problem, const int use_constr, con
 
 
 // Function to initialize the result
-int prima_init_result(prima_result_t *const result, const prima_problem_t problem)
+int prima_init_result(prima_result_t *const result, const prima_problem_t* problem)
 {
     if (!result)
         return PRIMA_NULL_RESULT;
@@ -109,19 +112,19 @@ int prima_init_result(prima_result_t *const result, const prima_problem_t proble
     result->message = NULL;
 
     // x: returned point
-    result->x = (double*)malloc(problem.n * sizeof(double));
+    result->x = (double*)malloc(problem->n * sizeof(double));
     if (!result->x)
         return PRIMA_MEMORY_ALLOCATION_FAILS;
-    for (int i = 0; i < problem.n; i++)
+    for (int i = 0; i < problem->n; i++)
         result->x[i] = NAN;
 
     // nlconstr: nonlinear constraint values at the returned point, of size m_nlcon (COBYLA only)
-    result->nlconstr = (double*)malloc(problem.m_nlcon * sizeof(double));
+    result->nlconstr = (double*)malloc(problem->m_nlcon * sizeof(double));
     if (!result->nlconstr) {
         free(result->x);
         return PRIMA_MEMORY_ALLOCATION_FAILS;
     }
-    for (int i = 0; i < problem.m_nlcon; i++)
+    for (int i = 0; i < problem->m_nlcon; i++)
         result->nlconstr[i] = NAN;
 
     return 0;
@@ -230,7 +233,7 @@ int lincoa_c(prima_obj_t calfun, const void *data, const int n, double x[], doub
 
 
 // The function that does the minimization using a PRIMA solver
-int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t problem, const prima_options_t options, prima_result_t *const result)
+int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t *problem, const prima_options_t *options, prima_result_t *const result)
 {
     int use_constr = (algorithm == PRIMA_COBYLA);
 
@@ -243,35 +246,35 @@ int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t prob
         // We copy x0 into result->x only after prima_check_problem has succeeded,
         // so that if prima_check_problem failed, result->x will not contained a
         // seemingly valid value.
-        for (int i = 0; i < problem.n; i++) {
-            result->x[i] = problem.x0[i];
+        for (int i = 0; i < problem->n; i++) {
+            result->x[i] = problem->x0[i];
         }
 
         switch (algorithm) {
             case PRIMA_BOBYQA:
-                bobyqa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), problem.xl, problem.xu, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, options.callback, &info);
+                bobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info);
                 result->cstrv = 0.0;
                 break;
 
             case PRIMA_COBYLA:
-                cobyla_c(problem.m_nlcon, problem.calcfc, options.data, problem.n, result->x, &(result->f), &(result->cstrv), result->nlconstr,
-                            problem.m_ineq, problem.Aineq, problem.bineq, problem.m_eq, problem.Aeq, problem.beq,
-                            problem.xl, problem.xu, problem.f0, problem.nlconstr0, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.iprint, options.ctol, options.callback, &info);
+                cobyla_c(problem->m_nlcon, problem->calcfc, options->data, problem->n, result->x, &(result->f), &(result->cstrv), result->nlconstr,
+                            problem->m_ineq, problem->Aineq, problem->bineq, problem->m_eq, problem->Aeq, problem->beq,
+                            problem->xl, problem->xu, problem->f0, problem->nlconstr0, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, options->ctol, options->callback, &info);
                 break;
 
             case PRIMA_LINCOA:
-                lincoa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), &(result->cstrv),
-                            problem.m_ineq, problem.Aineq, problem.bineq, problem.m_eq, problem.Aeq, problem.beq,
-                            problem.xl, problem.xu, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, options.ctol, options.callback, &info);
+                lincoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->cstrv),
+                            problem->m_ineq, problem->Aineq, problem->bineq, problem->m_eq, problem->Aeq, problem->beq,
+                            problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->ctol, options->callback, &info);
                 break;
 
             case PRIMA_NEWUOA:
-                newuoa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, options.callback, &info);
+                newuoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info);
                 result->cstrv = 0.0;
                 break;
 
             case PRIMA_UOBYQA:
-                uobyqa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.iprint, options.callback, &info);
+                uobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, options->callback, &info);
                 result->cstrv = 0.0;
                 break;
 
@@ -285,3 +288,225 @@ int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t prob
 
     return info;
 }
+
+prima_context_t* prima_context_create(const int n)
+{
+    size_t size = sizeof(prima_context_t);
+    prima_context_t* context = malloc(size);
+    if (context == NULL) {
+        return NULL;
+    }
+    memset(context, 0, size);
+    if (prima_init_problem(&context->problem, n) != 0 ||
+        prima_init_options(&context->options) != 0 ||
+        prima_init_result(&context->result, &context->problem) != 0) {
+        /* some error occurred */
+        prima_context_destroy(context);
+        return NULL;
+    }
+    return context;
+}
+
+void prima_context_destroy(prima_context_t* context)
+{
+    if (context != NULL) {
+        prima_free_result(&context->result);
+        free(context);
+    }
+}
+
+int prima_context_solve(prima_algorithm_t algorithm, prima_context_t* context)
+{
+    return prima_minimize(algorithm, &context->problem, &context->options, &context->result);
+}
+
+/* Accessors. */
+
+int prima_context_set_callfun(prima_context_t* context, prima_obj_t calfun, void* data)
+{
+    context->problem.calfun = calfun;
+    context->options.data = data;
+    return 0;
+}
+
+int prima_context_set_calcfc(prima_context_t* context, prima_objcon_t calcfc, void* data, int m_nlcon)
+{
+    context->problem.calcfc = calcfc;
+    context->problem.m_nlcon = m_nlcon;
+    context->options.data = data;
+    return 0;
+}
+
+int prima_context_set_x0(prima_context_t* context, double *x0)
+{
+    context->problem.x0 = x0;
+    return 0;
+}
+
+int prima_context_set_bounds(prima_context_t* context, double *xl, double *xu)
+{
+    context->problem.xl = xl;
+    context->problem.xu = xu;
+    return 0;
+}
+
+int prima_context_set_ineq(prima_context_t* context, int m_ineq, double Aineq[], double *bineq)
+{
+    context->problem.m_ineq = m_ineq;
+    context->problem.Aineq = Aineq;
+    context->problem.bineq = bineq;
+    return 0;
+}
+
+int prima_context_set_eq(prima_context_t* context, int m_eq, double Aeq[], double *beq)
+{
+    context->problem.m_eq = m_eq;
+    context->problem.Aeq = Aeq;
+    context->problem.beq = beq;
+    return 0;
+}
+
+int prima_context_set_f0(prima_context_t* context, double f0)
+{
+    context->problem.f0 = f0;
+    return 0;
+}
+
+int prima_context_set_nlconstr0(prima_context_t* context, double* nlconstr0)
+{
+    context->problem.nlconstr0 = nlconstr0;
+    return 0;
+}
+
+int prima_context_get_n(const prima_context_t* context)
+{
+    return context->problem.n;
+}
+
+int prima_context_get_m_ineq(const prima_context_t* context)
+{
+    return context->problem.m_ineq;
+}
+
+int prima_context_get_m_eq(const prima_context_t* context)
+{
+    return context->problem.m_eq;
+}
+
+int prima_context_get_m_nlcon(const prima_context_t* context)
+{
+    return context->problem.m_nlcon;
+}
+
+int prima_context_set_rho(prima_context_t* context, double rhobeg, double rhoend)
+{
+    context->options.rhobeg = rhobeg;
+    context->options.rhoend = rhoend;
+    return 0;
+}
+
+double prima_context_get_rhobeg(const prima_context_t* context)
+{
+    return context->options.rhobeg;
+}
+
+double prima_context_get_rhoend(const prima_context_t* context)
+{
+    return context->options.rhoend;
+}
+
+int prima_context_set_maxfun(prima_context_t* context, int maxfun)
+{
+    context->options.maxfun = maxfun;
+    return 0;
+}
+
+int prima_context_get_maxfun(const prima_context_t* context)
+{
+    return context->options.maxfun;
+}
+
+int prima_context_set_iprint(prima_context_t* context, int iprint)
+{
+    context->options.iprint = iprint;
+    return 0;
+}
+
+int prima_context_get_iprint(const prima_context_t* context)
+{
+    return context->options.iprint;
+}
+
+int prima_context_set_npt(prima_context_t* context, int npt)
+{
+    context->options.npt = npt;
+    return 0;
+}
+
+int prima_context_get_npt(const prima_context_t* context)
+{
+    return context->options.npt;
+}
+
+int prima_context_set_ftarget(prima_context_t* context, double ftarget)
+{
+    context->options.ftarget = ftarget;
+    return 0;
+}
+
+double prima_context_get_ftarget(const prima_context_t* context)
+{
+    return context->options.ftarget;
+}
+
+int prima_context_set_ctol(prima_context_t* context, double ctol)
+{
+    context->options.ctol = ctol;
+    return 0;
+}
+
+double prima_context_get_ctol(const prima_context_t* context)
+{
+    return context->options.ctol;
+}
+
+int prima_context_set_callback(prima_context_t* context, prima_callback_t callback)
+{
+    context->options.callback = callback;
+    return 0;
+}
+
+double* prima_context_get_x(const prima_context_t* context)
+{
+    return context->result.x;
+}
+
+double prima_context_get_f(const prima_context_t* context)
+{
+    return context->result.f;
+}
+
+double prima_context_get_cstrv(const prima_context_t* context)
+{
+    return context->result.cstrv;
+}
+
+double* prima_context_get_nlconstr(const prima_context_t* context)
+{
+    return context->result.nlconstr;
+}
+
+int prima_context_get_nf(const prima_context_t* context)
+{
+    return context->result.nf;
+}
+
+int prima_context_get_status(const prima_context_t* context)
+{
+    return context->result.status;
+}
+
+const char* prima_context_get_message(const prima_context_t* context)
+{
+    return context->result.message;
+}
diff --git a/c/tests/data.c b/c/tests/data.c
index 82f0b5bb26..6c09e87cde 100644
--- a/c/tests/data.c
+++ b/c/tests/data.c
@@ -144,7 +144,7 @@ int main(int argc, char * argv[])
 
   // Call the solver
   prima_result_t result;
-  int rc = prima_minimize(algorithm, problem, options, &result);
+  int rc = prima_minimize(algorithm, &problem, &options, &result);
 
   // Print the result
   printf("f* = %g, cstrv = %g, nlconstr = {%g}, rc = %d, msg = '%s', evals = %d\n", result.f, result.cstrv, result.nlconstr ? result.nlconstr[0] : 0.0, rc, result.message, result.nf);
diff --git a/c/tests/stress.c b/c/tests/stress.c
index be3323727f..31504bca11 100644
--- a/c/tests/stress.c
+++ b/c/tests/stress.c
@@ -87,7 +87,7 @@ int main(int argc, char * argv[])
   printf("Random seed = %d\n", seed);
   srand(seed);
 
-  // Set up the options 
+  // Set up the options
   prima_options_t options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_RHO;
@@ -168,7 +168,7 @@ int main(int argc, char * argv[])
 
   // Call the solver
   prima_result_t result;
-  rc = prima_minimize(algorithm, problem, options, &result);
+  rc = prima_minimize(algorithm, &problem, &options, &result);
 
   // Print the result
   printf("f* = %g, cstrv = %g, nlconstr = {%g}, rc = %d, msg = '%s', evals = %d\n", result.f, result.cstrv, result.nlconstr ? result.nlconstr[0] : 0.0, rc, result.message, result.nf);