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

Beyn's method test (fails) 

Former-commit-id: 535a94df45d1ddd6b23b6c94a8167ca969389099
This commit is contained in:
Marek Nečada 2019-08-23 11:16:46 +03:00
parent cf9892279d
commit 142d97d484
4 changed files with 50 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
qpms/beyn.c
View File

@ -417,7 +417,7 @@ int BeynSolve(BeynSolver *Solver, beyn_function_M_t M_function,
} }
*/ */
return 0; return K;
} }
/***************************************************************/ /***************************************************************/

5
qpms/beyn.h
View File

@ -33,6 +33,7 @@
#include <complex.h> #include <complex.h>
#include <gsl/gsl_matrix.h> #include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
/// User-supplied function that provides the operator M(z) whose "roots" are to be found. /// User-supplied function that provides the operator M(z) whose "roots" are to be found.
typedef int (*beyn_function_M_t)(gsl_matrix_complex *target_M, complex double z, void *params); typedef int (*beyn_function_M_t)(gsl_matrix_complex *target_M, complex double z, void *params);
@ -83,4 +84,8 @@ int BeynSolve(BeynSolver *Solver,
beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params, beyn_function_M_t M_function, beyn_function_M_inv_Vhat_t M_inv_Vhat_function, void *params,
double complex z0, double Rx, double Ry, int N); double complex z0, double Rx, double Ry, int N);
static inline complex double gsl_complex_tostd(gsl_complex z) { return GSL_REAL(z) + I*GSL_IMAG(z); }
static inline gsl_complex gsl_complex_fromstd(complex double z) { return gsl_complex_rect(creal(z), cimag(z)); }
#endif // BEYN_H #endif // BEYN_H

6
tests/CMakeLists.txt
View File

@ -15,7 +15,11 @@ add_executable( test_scalar_ewald32 test_scalar_ewald32.c )
target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m ) target_link_libraries( test_scalar_ewald32 qpms gsl lapacke amos m )
target_include_directories( test_scalar_ewald32 PRIVATE .. ) target_include_directories( test_scalar_ewald32 PRIVATE .. )
add_custom_target( mytests DEPENDS test_single_translations_vs_calc test_vswf_translations test_vswf_translations_array ) add_executable( tbeyn tbeyn.c )
target_link_libraries( tbeyn qpms gsl lapacke amos m )
target_include_directories( tbeyn 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 ) add_test( NAME single_vs_array_translation_coeffs COMMAND test_single_translations_vs_calc )

39
tests/tbeyn.c Normal file
View File

@ -0,0 +1,39 @@
#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_complex d = gsl_complex_fromstd( 2*m - 4*z / (6*m) );
gsl_complex od = gsl_complex_fromstd( -m - z / (6*m) );
gsl_matrix_complex_set_zero(target);
for (int i = 0; i < m; ++i) {
gsl_matrix_complex_set(target, i, i, (gsl_complex) d);
if(i > 0) gsl_matrix_complex_set(target, i, i-1, (gsl_complex) od);
if(i < m - 1) gsl_matrix_complex_set(target, i, i+1, (gsl_complex) od);
}
return 0;
}
int main() {
complex double z0 = 150+2*I;
double Rx = 148, Ry = 148;
int L = 10, N = 50, dim = 400;
BeynSolver * solver = CreateBeynSolver(dim, L);
int K = BeynSolve(solver, M_function, NULL /*M_inv_Vhat_function*/, NULL /*params*/,
z0, Rx, Ry, N);
printf("Found %d eigenvalues:\n", K);
for (int i = 0; i < K; ++i) {
gsl_complex eig = gsl_vector_complex_get(solver->Eigenvalues, i);
printf("%d: %g%+gj\n", i, GSL_REAL(eig), GSL_IMAG(eig));
}
DestroyBeynSolver(solver);
return 0;
}