aboutsummaryrefslogtreecommitdiff
path: root/src/hrs-scaling.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-01-16 17:03:08 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:11 +0100
commit87e90281aeb032410ede9c28f68c3a5d03e402fe (patch)
tree567b2793f23695723cf9b4ff306d5364e9861a7b /src/hrs-scaling.c
parent34901d9de3fc22f57ff318015fa62f35955f4d4f (diff)
Beginnings of Fox & Holmes method
Diffstat (limited to 'src/hrs-scaling.c')
-rw-r--r--src/hrs-scaling.c31
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);