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

Fix eigenvalue matrix dimensions in Beyn. 

Former-commit-id: 4d3b5269df29897ac09632bc0278a99dce2507d8
This commit is contained in:
Marek Nečada 2020-04-12 00:48:13 +03:00
parent 63aa338891
commit c5c148ca40
3 changed files with 42 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
qpms/beyn.c
View File

@ -276,7 +276,7 @@ BeynSolver *BeynSolver_create(int M, int L)
solver->eigenvalues = gsl_vector_complex_calloc(L);
solver->eigenvalue_errors = gsl_vector_complex_calloc(L);
solver->residuals = gsl_vector_calloc(L);
solver->eigenvectors = gsl_matrix_complex_calloc(M, L);
solver->eigenvectors = gsl_matrix_complex_calloc(L, M);
// storage for singular values, random VHat matrix, etc. used in algorithm
solver->A0 = gsl_matrix_complex_calloc(M,L);

4
tests/CMakeLists.txt
View File

@ -76,6 +76,10 @@ add_executable( tbeyn_gsl tbeyn_gsl.c )
target_link_libraries( tbeyn_gsl qpms gsl lapacke amos m )
target_include_directories( tbeyn_gsl PRIVATE .. )
add_executable( tbeyn_gsl2 tbeyn_gsl2.c )
target_link_libraries( tbeyn_gsl2 qpms gsl lapacke amos m )
target_include_directories( tbeyn_gsl2 PRIVATE .. )
add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array tbeyn )
add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc )

37
tests/tbeyn_gsl2.c Normal file
View File

@ -0,0 +1,37 @@
#include <qpms/beyn.h>
#include <stdio.h>
// Matrix as in Beyn, section 4.11
int M_function(gsl_matrix_complex *target, complex double z, void *no_params) {
int m = target->size1;
gsl_matrix_complex_set_zero(target);
for (int i = 0; i < m; ++i) {
gsl_matrix_complex_set(target, i, i, gsl_complex_rect(i - creal(z), -cimag(z)));
}
return 0;
}
int main() {
complex double z0 = 150+2*I;
double Rx = 5, Ry = 5;
int L = 10, N = 600, dim = 200;
beyn_contour_t *contour = beyn_contour_ellipse(z0, Rx, Ry, N);
beyn_result_gsl_t *result =
beyn_solve_gsl(dim, L, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
contour, 1e-4, 1, 1e-4);
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("%zd: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
}
free(contour);
beyn_result_gsl_free(result);
return 0;
}