Alternative approach in recursive "Delta" to avoid overflows.
[qpms.git] / tests / tbeyn_gsl2.c
blob1da782d3dfc34853831c8755177c0e38bd27a271
1 #include <qpms/beyn.h>
2 #include <stdio.h>
5 // Matrix as in Beyn, section 4.11
6 int M_function(gsl_matrix_complex *target, complex double z, void *no_params) {
7 int m = target->size1;
9 gsl_matrix_complex_set_zero(target);
10 for (int i = 0; i < m; ++i) {
11 gsl_matrix_complex_set(target, i, i, gsl_complex_rect(i - creal(z), -cimag(z)));
14 return 0;
17 int main() {
18 complex double z0 = 150+2*I;
19 double Rx = 5, Ry = 5;
20 int L = 10, N = 600, dim = 200;
21 beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
23 beyn_result_gsl_t *result =
24 beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
25 contour, 1e-4, 1, 1e-4);
26 printf("Found %zd eigenvalues:\n", result->neig);
27 for (size_t i = 0; i < result->neig; ++i) {
28 gsl_complex eig = gsl_vector_complex_get(result->eigval, i);
29 printf("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
32 free(contour);
33 beyn_result_gsl_free(result);
34 return 0;