QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Beyn minor fixes, cleaner API. 

Former-commit-id: 44d7147a14194a7b8e5a66552dd3855029b9e370
This commit is contained in:
Marek Nečada 2019-09-10 14:54:08 +03:00
parent 3bc59096bc
commit 11aa48d2da
5 changed files with 32 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
qpms/beyn.c
View File

@ -284,6 +284,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun
gsl_matrix_complex_free(V0); gsl_matrix_complex_free(V0);
// FIXME!!! make clear relation between KRetained and K in the results!
int KRetained = 0; int KRetained = 0;
gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m); gsl_matrix_complex *Mmat = gsl_matrix_complex_alloc(m, m);
gsl_vector_complex *MVk = gsl_vector_complex_alloc(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, 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, void *params, const beyn_contour_t *contour,
double rank_tol, double res_tol) 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); gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors);
// TODO Original did also fabs on the complex and real parts ^^^. // TODO Original did also fabs on the complex and real parts ^^^.
QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_gsl_t)); beyn_result_gsl_t *result;
(*result)->eigval = solver->eigenvalues; QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t));
(*result)->eigval_err = solver->eigenvalue_errors; result->eigval = solver->eigenvalues;
(*result)->residuals = solver->residuals; result->eigval_err = solver->eigenvalue_errors;
(*result)->eigvec = solver->eigenvectors; result->residuals = solver->residuals;
result->eigvec = solver->eigenvectors;
result->neig = K;
BeynSolver_free_partial(solver); BeynSolver_free_partial(solver);
return K; return result;
} }
// Wrapper of pure C array M-matrix function to GSL. // 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); (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) { 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}; struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params};
QPMS_CRASHING_MALLOC(*result, sizeof(beyn_result_t)); return beyn_result_from_beyn_result_gsl(
struct beyn_result_gsl_t *result_gsl; beyn_solve_gsl(m, l, beyn_function_M_carr2gsl,
int retval = beyn_solve_gsl(&result_gsl, m, l, beyn_function_M_carr2gsl,
(p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL,
(void *) &p, contour, rank_tol, res_tol); (void *) &p, contour, rank_tol, res_tol)
*result = beyn_result_from_beyn_result_gsl(result_gsl); );
return retval;
} }
beyn_result_t *beyn_result_from_beyn_result_gsl(beyn_result_gsl_t *g) { 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)); QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_t));
result->gsl = g; result->gsl = g;
result->neig = result->gsl->neig; result->neig = result->gsl->neig;
result->vlen = result->gsl->eigvec->size2;
result->eigval = (complex double *) result->gsl->eigval->data; result->eigval = (complex double *) result->gsl->eigval->data;
result->eigval_err = (complex double *) result->gsl->eigval_err->data; result->eigval_err = (complex double *) result->gsl->eigval_err->data;
result->residuals = result->gsl->residuals->data; result->residuals = result->gsl->residuals->data;

5
qpms/beyn.h
View File

@ -51,6 +51,7 @@ void beyn_result_gsl_free(beyn_result_gsl_t *result);
/// Beyn algorithm result structure (pure C array version). /// Beyn algorithm result structure (pure C array version).
typedef struct beyn_result_t { typedef struct beyn_result_t {
size_t neig; ///< Number of eigenvalues found. 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;
complex double *eigval_err; complex double *eigval_err;
double *residuals; 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); 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 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). 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$. 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. 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 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). 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$. beyn_function_M_t M, ///< Function providing the matrix \f$ M(z) \f$.

5
qpms/qpms_cdefs.pxd
View File

@ -554,6 +554,7 @@ cdef extern from "beyn.h":
pass pass
ctypedef struct beyn_result_t: ctypedef struct beyn_result_t:
size_t neig size_t neig
size_t vlen
cdouble *eigval cdouble *eigval
cdouble *eigval_err cdouble *eigval_err
double *residuals double *residuals
@ -568,11 +569,11 @@ cdef extern from "beyn.h":
void beyn_result_gsl_free(beyn_result_gsl_t *result) void beyn_result_gsl_free(beyn_result_gsl_t *result)
void beyn_result_free(beyn_result_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, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, void *params, const beyn_contour_t *contour,
double rank_tol, double res_tol) 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, beyn_function_M_inv_Vhat_t M_inv_Vhat, void *params, const beyn_contour_t *contour,
double rank_tol, double res_tol) double rank_tol, double res_tol)

10
tests/tbeyn.c
View File

@ -24,13 +24,13 @@ int main() {
int L = 10, N = 50, dim = 400; int L = 10, N = 50, dim = 400;
beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
beyn_result_t *result; beyn_result_t *result =
int K = beyn_solve(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, beyn_solve(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
contour, 1e-4, 0); contour, 1e-4, 0);
printf("Found %d eigenvalues:\n", K); printf("Found %zd eigenvalues:\n", result->neig);
for (int i = 0; i < K; ++i) { for (size_t i = 0; i < result->neig; ++i) {
complex double eig = result->eigval[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); free(contour);
beyn_result_free(result); beyn_result_free(result);

10
tests/tbeyn_gsl.c
View File

@ -26,13 +26,13 @@ int main() {
int L = 10, N = 50, dim = 400; int L = 10, N = 50, dim = 400;
beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N); beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
beyn_result_gsl_t *result; beyn_result_gsl_t *result =
int K = beyn_solve_gsl(&result, dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/, beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
contour, 1e-4, 0); contour, 1e-4, 0);
printf("Found %d eigenvalues:\n", K); printf("Found %zd eigenvalues:\n", result->neig);
for (int i = 0; i < K; ++i) { for (size_t i = 0; i < result->neig; ++i) {
gsl_complex eig = gsl_vector_complex_get(result->eigval, 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); free(contour);
beyn_result_gsl_free(result); beyn_result_gsl_free(result);