From 11aa48d2da0b90822adc5d86a4f0e322bba78b45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 10 Sep 2019 14:54:08 +0300 Subject: [PATCH] Beyn minor fixes, cleaner API. Former-commit-id: 44d7147a14194a7b8e5a66552dd3855029b9e370 --- qpms/beyn.c | 30 ++++++++++++++++-------------- qpms/beyn.h | 5 +++-- qpms/qpms_cdefs.pxd | 5 +++-- tests/tbeyn.c | 10 +++++----- tests/tbeyn_gsl.c | 10 +++++----- 5 files changed, 32 insertions(+), 28 deletions(-) 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);