Beginnings of Fox & Holmes method

1 files changed, 31 insertions, 0 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index ea87dec7..03f4c1d5 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -19,6 +19,7 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
+#include <gsl/gsl_eigen.h>
 
 #include "image.h"
 #include "peaks.h"
@@ -45,6 +46,20 @@ static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
 }
 
 
+static void show_eigen(gsl_matrix *e_vec, gsl_vector *e_val, int r)
+{
+	int i, j;
+
+	for ( i=0; i<r; i++ ) {
+		STATUS("[ ");
+		for ( j=0; j<r; j++ ) {
+			STATUS("%+5.2f ", gsl_matrix_get(e_vec, i, j));
+		}
+		STATUS("]   [ %+9.3e ]\n", gsl_vector_get(e_val, i));
+	}
+}
+
+
 static void sum_GI(struct image *images, int n, const char *sym,
                    signed int hat, signed int kat, signed int lat,
                    double *sigma_GI, double *sigma_Gsq)
@@ -205,6 +220,21 @@ static double iterate_scale(struct image *images, int n,
 	}
 	show_matrix_eqn(M, v, n);
 
+	gsl_eigen_symmv_workspace *work;
+	gsl_vector *e_val;
+	gsl_matrix *e_vec;
+	int val;
+
+	work = gsl_eigen_symmv_alloc(n);
+	e_val = gsl_vector_alloc(n);
+	e_vec = gsl_matrix_alloc(n, n);
+	val = gsl_eigen_symmv(M, e_val, e_vec, work);
+	STATUS("gsl_eigen_symmv said %i (%s)\n", val, gsl_strerror(val));
+	gsl_eigen_symmv_free(work);
+
+	show_eigen(e_vec, e_val, n);
+
+#if 0  /* HRS method */
 	shifts = gsl_vector_alloc(n);
 	gsl_linalg_HH_solve(M, v, shifts);
 	max_shift = 0.0;
@@ -220,6 +250,7 @@ static double iterate_scale(struct image *images, int n,
 
 	}
 	gsl_vector_free(shifts);
+#endif
 
 	gsl_matrix_free(M);
 	gsl_vector_free(v);