diff options
| author | Thomas White <taw@physics.org> | 2011-02-03 10:40:18 +0100 |
|---|---|---|
| committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:13 +0100 |
| commit | 20e6c9e8f6055d4c98cc32f8e4dd59e6180b416c (patch) | |
| tree | 84378acdacf1e8031ddebcf6f971f26dc55a0312 /src/hrs-scaling.c | |
| parent | 6854d52ef4e7268d50e90ea112af73e42a34a829 (diff) | |
Solve the diagonalised equation the easy way
Diffstat (limited to 'src/hrs-scaling.c')
| -rw-r--r-- | src/hrs-scaling.c | 18 |
1 files changed, 8 insertions, 10 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 30af5a61..73321afb 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -242,21 +242,20 @@ static double iterate_scale(struct image *images, int n, gsl_eigen_symmv_free(work); val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC); - /* Set up the diagonalised normal equations */ - gsl_matrix *D; + /* Rotate the "solution vector" */ gsl_vector *rprime; - D = gsl_matrix_calloc(n, n); - for ( a=0; a<n; a++ ) { - gsl_matrix_set(D, a, a, gsl_vector_get(e_val, a)); - } rprime = gsl_vector_alloc(n); val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime); - show_matrix_eqn(D, rprime, n); - /* Solve */ + /* Solve (now easy) */ gsl_vector *sprime; sprime = gsl_vector_alloc(n); - gsl_linalg_HH_solve(D, rprime, sprime); + for ( a=0; a<n-1; a++ ) { + double num, den; + num = gsl_vector_get(rprime, a); + den = gsl_vector_get(e_val, a); + gsl_vector_set(sprime, a, num/den); + } gsl_vector_set(sprime, n-1, 0.0); /* Set last shift to zero */ /* Rotate back */ @@ -281,7 +280,6 @@ static double iterate_scale(struct image *images, int n, gsl_vector_free(sprime); gsl_matrix_free(e_vec); gsl_vector_free(e_val); - gsl_matrix_free(D); gsl_matrix_free(M); gsl_vector_free(v); |