diff --git a/qpms/beyn.c b/qpms/beyn.c
index 3c3b1ed..d2495af 100644
--- a/qpms/beyn.c
+++ b/qpms/beyn.c
@@ -284,6 +284,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
 
   gsl_matrix_complex_free(V0);
 
+  // FIXME!!! make clear relation between KRetained and K in the results!
   int KRetained = 0;
   gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m);
   gsl_vector_complex *MVk = gsl_vector_complex_alloc(m);
@@ -319,7 +320,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
 }
 
 
-int beyn_solve_gsl(beyn_result_gsl_t **result, const size_t m, const size_t l,
+beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l,
     beyn_function_M_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function,
     void *params, const beyn_contour_t *contour,
     double rank_tol, double res_tol)
@@ -399,15 +400,17 @@ int beyn_solve_gsl(beyn_result_gsl_t **result, const size_t m, const size_t l,
   gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors);
   // TODO Original did also fabs on the complex and real parts ^^^.
 
-  QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_gsl_t));
-  (*result)->eigval = solver->eigenvalues;
-  (*result)->eigval_err = solver->eigenvalue_errors;
-  (*result)->residuals = solver->residuals;
-  (*result)->eigvec = solver->eigenvectors;
+  beyn_result_gsl_t *result;
+  QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t));
+  result->eigval = solver->eigenvalues;
+  result->eigval_err = solver->eigenvalue_errors;
+  result->residuals = solver->residuals;
+  result->eigvec = solver->eigenvectors;
+  result->neig = K;
 
   BeynSolver_free_partial(solver);
 
-  return K;
+  return result;
 }
 
 // Wrapper of pure C array M-matrix function to GSL.
@@ -440,16 +443,14 @@ static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target,
       (complex double *) Vhat->data, z, p->param);
 }
 
-int beyn_solve(beyn_result_t **result, size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat,
+beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat,
     void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) {
   struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params};
-  QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_t));
-  struct beyn_result_gsl_t *result_gsl;
-  int retval = beyn_solve_gsl(&result_gsl, m, l, beyn_function_M_carr2gsl,
+  return beyn_result_from_beyn_result_gsl(
+    beyn_solve_gsl(m, l, beyn_function_M_carr2gsl,
       (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, 
-      (void *) &p, contour, rank_tol, res_tol);
-  *result = beyn_result_from_beyn_result_gsl(result_gsl);
-  return retval;
+      (void *) &p, contour, rank_tol, res_tol)
+    );
 }
 
 beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) {
@@ -457,6 +458,7 @@ beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) {
   QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_t));
   result->gsl = g;
   result->neig = result->gsl->neig;
+  result->vlen = result->gsl->eigvec->size2;
   result->eigval = (complex double *) result->gsl->eigval->data;
   result->eigval_err = (complex double *) result->gsl->eigval_err->data;
   result->residuals = result->gsl->residuals->data;
diff --git a/qpms/beyn.h b/qpms/beyn.h
index b32025e..cbcb64d 100644
--- a/qpms/beyn.h
+++ b/qpms/beyn.h
@@ -51,6 +51,7 @@ void beyn_result_gsl_free(beyn_result_gsl_t *result);
 /// Beyn algorithm result structure (pure C array version).
 typedef struct beyn_result_t {
 	size_t neig; ///< Number of eigenvalues found.
+	size_t vlen; ///< Vector space dimension (also the leading dimension of eigvec).
 	complex double *eigval;
 	complex double *eigval_err;
 	double *residuals;
@@ -69,7 +70,7 @@ beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g);
 
 void beyn_result_free(beyn_result_t *result);
 
-int beyn_solve_gsl(beyn_result_gsl_t **result,
+beyn_result_gsl_t *beyn_solve_gsl(
 		size_t m, ///< Dimension of the matrix \a M.
 		size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
 		beyn_function_M_gsl_t M, ///< Function providing the matrix \f$ M(z) \f$.
@@ -80,7 +81,7 @@ int beyn_solve_gsl(beyn_result_gsl_t **result,
 		double res_tol ///< (default: `0.0`) TODO DOC.
 	      );
 
-int beyn_solve(beyn_result_t **result,
+beyn_result_t *beyn_solve(
 		size_t m, ///< Dimension of the matrix \a M.
 		size_t l, ///< Number of columns of the random matrix  \f$ \hat V \f$ (larger than the expected number of solutions).
 		beyn_function_M_t M, ///< Function providing the matrix \f$ M(z) \f$.
diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd
index 941e51c..046c343 100644
--- a/qpms/qpms_cdefs.pxd
+++ b/qpms/qpms_cdefs.pxd
@@ -554,6 +554,7 @@ cdef extern from "beyn.h":
         pass
     ctypedef struct beyn_result_t:
         size_t neig
+        size_t vlen
         cdouble *eigval
         cdouble *eigval_err
         double *residuals
@@ -568,11 +569,11 @@ cdef extern from "beyn.h":
     void beyn_result_gsl_free(beyn_result_gsl_t *result)
     void beyn_result_free(beyn_result_t *result)
 
-    int beyn_solve_gsl(beyn_result_gsl_t **result, size_t m, size_t l, beyn_function_M_gsl_t M,
+    beyn_result_gsl_t *beyn_solve_gsl(size_t m, size_t l, beyn_function_M_gsl_t M,
             beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, void *params, const beyn_contour_t *contour,
             double rank_tol, double res_tol)
 
-    int beyn_solve(beyn_result_t **result, size_t m, size_t l, beyn_function_M_t M,
+    beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M,
             beyn_function_M_inv_Vhat_t M_inv_Vhat, void *params, const beyn_contour_t *contour,
             double rank_tol, double res_tol)
 
diff --git a/tests/tbeyn.c b/tests/tbeyn.c
index 1e0322c..8cdb4ea 100644
--- a/tests/tbeyn.c
+++ b/tests/tbeyn.c
@@ -24,13 +24,13 @@ int main() {
   int L = 10, N = 50, dim = 400;
   beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
 
-  beyn_result_t *result;
-  int K = beyn_solve(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
+  beyn_result_t *result =
+    beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
       contour, 1e-4, 0);
-  printf("Found %d eigenvalues:\n", K);
-  for (int i = 0; i < K; ++i) {
+  printf("Found %zd eigenvalues:\n", result->neig);
+  for (size_t i = 0; i < result->neig; ++i) {
     complex double eig = result->eigval[i];
-    printf("%d: %g%+gj\n", i, creal(eig), cimag(eig));
+    printf("%zd: %g%+gj\n", i, creal(eig), cimag(eig));
   }
   free(contour);
   beyn_result_free(result);
diff --git a/tests/tbeyn_gsl.c b/tests/tbeyn_gsl.c
index 8d492d3..36c46a7 100644
--- a/tests/tbeyn_gsl.c
+++ b/tests/tbeyn_gsl.c
@@ -26,13 +26,13 @@ int main() {
   int L = 10, N = 50, dim = 400;
   beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
 
-  beyn_result_gsl_t *result;
-  int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
+  beyn_result_gsl_t *result =
+    beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
       contour, 1e-4, 0);
-  printf("Found %d eigenvalues:\n", K);
-  for (int i = 0; i < K; ++i) {
+  printf("Found %zd eigenvalues:\n", result->neig);
+  for (size_t i = 0; i < result->neig; ++i) {
     gsl_complex eig = gsl_vector_complex_get(result->eigval, i);
-    printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
+    printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
   }
   free(contour);
   beyn_result_gsl_free(result);