Revert "Remove row-crossing stuff since it should not be needed"

This reverts commit a043b3cb7cd57da3f5f7feb9d470d1e771f5b664.

Conflicts:

	src/hrs-scaling.c

1 files changed, 24 insertions, 8 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 86a5ffab..eee459a2 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -137,19 +137,23 @@ static double iterate_scale(struct image *images, int n,
 	gsl_vector *shifts;
 	int l;
 	double max_shift;
+	int crossed = 0;
 	int n_ref;
 
-	M = gsl_matrix_calloc(n, n);
-	v = gsl_vector_calloc(n);
+	M = gsl_matrix_calloc(n-1, n-1);
+	v = gsl_vector_calloc(n-1);
 	n_ref = num_items(obs);
 
 	for ( l=0; l<n; l++ ) {  /* "Equation number": one equation per frame */
 
 		int m;  /* Frame (scale factor) number */
 		int h;
+		int lcomp;
 		double vc_tot = 0.0;
 		struct image *imagel = &images[l];
 
+		if ( l == crossed ) continue;
+
 		/* Determine the "solution" vector component */
 		for ( h=0; h<n_ref; h++ ) {
 
@@ -172,15 +176,19 @@ static double iterate_scale(struct image *images, int n,
 
 		}
 
-		gsl_vector_set(v, l, vc_tot);
+		lcomp = l;
+		if ( l > crossed ) lcomp--;
+		gsl_vector_set(v, lcomp, vc_tot);
 
 		/* Now fill in the matrix components */
 		for ( m=0; m<n; m++ ) {
 
 			double mc_tot = 0.0;
+			int mcomp;
 			struct image *imagem = &images[m];
 
 			if ( m > l ) continue;  /* Matrix is symmetric */
+			if ( m == crossed ) continue;
 
 			for ( h=0; h<n_ref; h++ ) {
 
@@ -211,14 +219,17 @@ static double iterate_scale(struct image *images, int n,
 
 			}
 
-			gsl_matrix_set(M, l, m, mc_tot);
-			gsl_matrix_set(M, m, l, mc_tot);
+			mcomp = m;
+			if ( m > crossed ) mcomp--;
+			gsl_matrix_set(M, lcomp, mcomp, mc_tot);
+			gsl_matrix_set(M, mcomp, lcomp, mc_tot);
 
 		}
 
 
 	}
-	show_matrix_eqn(M, v, n);
+	show_matrix_eqn(M, v, n-1);
+
 
 	gsl_eigen_symmv_workspace *work;
 	gsl_vector *e_val;
@@ -235,14 +246,19 @@ static double iterate_scale(struct image *images, int n,
 	show_eigen(e_vec, e_val, n);
 
 #if 0  /* HRS method */
-	shifts = gsl_vector_alloc(n);
+	shifts = gsl_vector_alloc(n-1);
+
 	gsl_linalg_HH_solve(M, v, shifts);
 	max_shift = 0.0;
 	for ( l=0; l<n-1; l++ ) {
 
 		double shift = gsl_vector_get(shifts, l);
+		int limg;
+
+		limg = l;
+		if ( limg >= crossed ) limg++;
 
-		images[l].osf += shift;
+		images[limg].osf += shift;
 
 		if ( fabs(shift) > fabs(max_shift) ) {
 			max_shift = fabs(shift);