From c5c148ca4099689368dfacf004814b01874d744a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Sun, 12 Apr 2020 00:48:13 +0300
Subject: [PATCH] Fix eigenvalue matrix dimensions in Beyn.

Former-commit-id: 4d3b5269df29897ac09632bc0278a99dce2507d8
---
 qpms/beyn.c          |  2 +-
 tests/CMakeLists.txt |  4 ++++
 tests/tbeyn_gsl2.c   | 37 +++++++++++++++++++++++++++++++++++++
 3 files changed, 42 insertions(+), 1 deletion(-)
 create mode 100644 tests/tbeyn_gsl2.c

diff --git a/qpms/beyn.c b/qpms/beyn.c
index b04649f..36c5605 100644
--- a/qpms/beyn.c
+++ b/qpms/beyn.c
@@ -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);
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 17ded6f..ccd039a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -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 )
diff --git a/tests/tbeyn_gsl2.c b/tests/tbeyn_gsl2.c
new file mode 100644
index 0000000..1da782d
--- /dev/null
+++ b/tests/tbeyn_gsl2.c
@@ -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;
+}
+
+