Change matrix notation to match ITA chapter 5.1

6 files changed, 88 insertions, 81 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index cc17f108..605aa449 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -405,6 +405,10 @@ static IntegerMatrix *centering_transformation(UnitCell *in,
 	ua = cell_get_unique_axis(in);
 	cen = cell_get_centering(in);
 
+	/* Write the matrices exactly as they appear in ITA Table 5.1.3.1.
+	 * C is "P", and Ci is "Q=P^-1".  Vice-versa if the transformation
+	 * should go the opposite way to what's written in the first column. */
+
 	if ( (cen=='P') || (cen=='R') ) {
 		*new_centering = cen;
 		*new_latt = lt;
@@ -451,12 +455,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in,
 
 	if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) {
 		/* Obverse setting */
-		C = intmat_create_3x3(1, -1,  0,
-		                      0,  1, -1,
-		                      1,  1,  1);
-		Ci = create_rtnl_mtx( 2,3,  1,3,  1,3,
-		                     -1,3,  1,3,  1,3,
-		                     -1,3, -2,3,  1,3);
+		C = intmat_create_3x3( 1,  0,  1,
+		                      -1,  1,  1,
+		                       0, -1,  1);
+		Ci = create_rtnl_mtx( 2,3, -1,3, -1,3,
+		                      1,3,  1,3, -2,3,
+		                      1,3,  1,3,  1,3);
 		assert(lt == L_HEXAGONAL);
 		assert(ua == 'c');
 		*new_latt = L_RHOMBOHEDRAL;
@@ -465,12 +469,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in,
 	}
 
 	if ( cen == 'A' ) {
-		C = intmat_create_3x3(0,  0, -1,
-		                      1,  1,  0,
-		                      1, -1,  0);
-		Ci = create_rtnl_mtx( 0,1,  1,2,  1,2,
+		C = intmat_create_3x3( 1,  0,  0,
+		                       0,  1,  1,
+		                       0, -1,  1);
+		Ci = create_rtnl_mtx( 1,1,  0,1,  0,1,
 		                      0,1,  1,2, -1,2,
-		                     -1,1,  0,1,  0,1);
+		                      0,1,  1,2,  1,2);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
 			*new_centering = 'P';
@@ -483,12 +487,12 @@ static IntegerMatrix *centering_transformation(UnitCell *in,
 	}
 
 	if ( cen == 'B' ) {
-		C = intmat_create_3x3(1,  1,  0,
-		                      0,  0,  1,
-		                      1, -1,  0);
-		Ci = create_rtnl_mtx( 1,2,  0,1,  1,2,
-		                      1,2,  0,1, -1,2,
-		                      0,1,  1,1,  0,1);
+		C = intmat_create_3x3( 1,  0,  1,
+		                       0,  1,  0,
+		                      -1,  0,  1);
+		Ci = create_rtnl_mtx( 1,2,  0,1, -1,2,
+		                      0,1,  1,1,  0,1,
+		                      1,2,  0,1,  1,2);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
 			*new_centering = 'P';
@@ -501,11 +505,11 @@ static IntegerMatrix *centering_transformation(UnitCell *in,
 	}
 
 	if ( cen == 'C' ) {
-		C = intmat_create_3x3(1, -1,  0,
-		                      1,  1,  0,
-		                      0,  0,  1);
-		Ci = create_rtnl_mtx( 1,2,  1,2,  0,1,
-		                     -1,2,  1,2,  0,1,
+		C = intmat_create_3x3( 1,  1,  0,
+		                      -1,  1,  0,
+		                       0,  0,  1);
+		Ci = create_rtnl_mtx( 1,2, -1,2,  0,1,
+		                      1,2,  1,2,  0,1,
 		                      0,1,  0,1,  1,1);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
@@ -2047,13 +2051,13 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 
 			/* Form the matrix using the first candidate for c */
 			rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
-			rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]);
-			rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]);
-			rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]);
+			rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]);
+			rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]);
+			rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]);
 			rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
-			rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]);
-			rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]);
-			rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]);
+			rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]);
+			rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]);
+			rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]);
 			rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
 
 			/* Check angle between a and b */
@@ -2066,13 +2070,13 @@ int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
 			for ( ic=0; ic<ncand_c; ic++ ) {
 
 				rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
-				rtnl_mtx_set(m, 0, 1, cand_a[3*ia+1]);
-				rtnl_mtx_set(m, 0, 2, cand_a[3*ia+2]);
-				rtnl_mtx_set(m, 1, 0, cand_b[3*ib+0]);
+				rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]);
+				rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]);
+				rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]);
 				rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
-				rtnl_mtx_set(m, 1, 2, cand_b[3*ib+2]);
-				rtnl_mtx_set(m, 2, 0, cand_c[3*ic+0]);
-				rtnl_mtx_set(m, 2, 1, cand_c[3*ic+1]);
+				rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]);
+				rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]);
+				rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]);
 				rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
 
 				if ( rtnl_cmp(rtnl_mtx_det(m),rtnl_zero()) == 0 ) continue;
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 130a24b1..5e00ac9c 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -368,13 +368,13 @@ static int cell_invert(double ax, double ay, double az,
 		return 1;
 	}
 	gsl_matrix_set(m, 0, 0, ax);
-	gsl_matrix_set(m, 0, 1, bx);
-	gsl_matrix_set(m, 0, 2, cx);
 	gsl_matrix_set(m, 1, 0, ay);
-	gsl_matrix_set(m, 1, 1, by);
-	gsl_matrix_set(m, 1, 2, cy);
 	gsl_matrix_set(m, 2, 0, az);
+	gsl_matrix_set(m, 0, 1, bx);
+	gsl_matrix_set(m, 1, 1, by);
 	gsl_matrix_set(m, 2, 1, bz);
+	gsl_matrix_set(m, 0, 2, cx);
+	gsl_matrix_set(m, 1, 2, cy);
 	gsl_matrix_set(m, 2, 2, cz);
 
 	/* Invert */
@@ -410,13 +410,13 @@ static int cell_invert(double ax, double ay, double az,
 	gsl_matrix_transpose(inv);
 
 	*asx = gsl_matrix_get(inv, 0, 0);
-	*bsx = gsl_matrix_get(inv, 0, 1);
-	*csx = gsl_matrix_get(inv, 0, 2);
 	*asy = gsl_matrix_get(inv, 1, 0);
-	*bsy = gsl_matrix_get(inv, 1, 1);
-	*csy = gsl_matrix_get(inv, 1, 2);
 	*asz = gsl_matrix_get(inv, 2, 0);
+	*bsx = gsl_matrix_get(inv, 0, 1);
+	*bsy = gsl_matrix_get(inv, 1, 1);
 	*bsz = gsl_matrix_get(inv, 2, 1);
+	*csx = gsl_matrix_get(inv, 0, 2);
+	*csy = gsl_matrix_get(inv, 1, 2);
 	*csz = gsl_matrix_get(inv, 2, 2);
 
 	gsl_matrix_free(inv);
@@ -630,27 +630,27 @@ UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m)
 
 	c = gsl_matrix_alloc(3, 3);
 	gsl_matrix_set(c, 0, 0, asx);
-	gsl_matrix_set(c, 0, 1, asy);
-	gsl_matrix_set(c, 0, 2, asz);
-	gsl_matrix_set(c, 1, 0, bsx);
+	gsl_matrix_set(c, 1, 0, asy);
+	gsl_matrix_set(c, 2, 0, asz);
+	gsl_matrix_set(c, 0, 1, bsx);
 	gsl_matrix_set(c, 1, 1, bsy);
-	gsl_matrix_set(c, 1, 2, bsz);
-	gsl_matrix_set(c, 2, 0, csx);
-	gsl_matrix_set(c, 2, 1, csy);
+	gsl_matrix_set(c, 2, 1, bsz);
+	gsl_matrix_set(c, 0, 2, csx);
+	gsl_matrix_set(c, 1, 2, csy);
 	gsl_matrix_set(c, 2, 2, csz);
 
 	res = gsl_matrix_calloc(3, 3);
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c, m, 0.0, res);
 
 	out = cell_new_from_cell(in);
 	cell_set_cartesian(out, gsl_matrix_get(res, 0, 0),
-	                        gsl_matrix_get(res, 0, 1),
-	                        gsl_matrix_get(res, 0, 2),
 	                        gsl_matrix_get(res, 1, 0),
-	                        gsl_matrix_get(res, 1, 1),
-	                        gsl_matrix_get(res, 1, 2),
 	                        gsl_matrix_get(res, 2, 0),
+	                        gsl_matrix_get(res, 0, 1),
+	                        gsl_matrix_get(res, 1, 1),
 	                        gsl_matrix_get(res, 2, 1),
+	                        gsl_matrix_get(res, 0, 2),
+	                        gsl_matrix_get(res, 1, 2),
 	                        gsl_matrix_get(res, 2, 2));
 
 	gsl_matrix_free(res);
@@ -673,27 +673,27 @@ UnitCell *cell_transform_gsl_reciprocal(UnitCell *in, gsl_matrix *m)
 
 	c = gsl_matrix_alloc(3, 3);
 	gsl_matrix_set(c, 0, 0, asx);
-	gsl_matrix_set(c, 0, 1, asy);
-	gsl_matrix_set(c, 0, 2, asz);
-	gsl_matrix_set(c, 1, 0, bsx);
+	gsl_matrix_set(c, 1, 0, asy);
+	gsl_matrix_set(c, 2, 0, asz);
+	gsl_matrix_set(c, 0, 1, bsx);
 	gsl_matrix_set(c, 1, 1, bsy);
-	gsl_matrix_set(c, 1, 2, bsz);
-	gsl_matrix_set(c, 2, 0, csx);
-	gsl_matrix_set(c, 2, 1, csy);
+	gsl_matrix_set(c, 2, 1, bsz);
+	gsl_matrix_set(c, 0, 2, csx);
+	gsl_matrix_set(c, 1, 2, csy);
 	gsl_matrix_set(c, 2, 2, csz);
 
 	res = gsl_matrix_calloc(3, 3);
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c, m, 0.0, res);
 
 	out = cell_new_from_cell(in);
 	cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0),
-	                         gsl_matrix_get(res, 0, 1),
-	                         gsl_matrix_get(res, 0, 2),
 	                         gsl_matrix_get(res, 1, 0),
-	                         gsl_matrix_get(res, 1, 1),
-	                         gsl_matrix_get(res, 1, 2),
 	                         gsl_matrix_get(res, 2, 0),
+	                         gsl_matrix_get(res, 0, 1),
+	                         gsl_matrix_get(res, 1, 1),
 	                         gsl_matrix_get(res, 2, 1),
+	                         gsl_matrix_get(res, 0, 2),
+	                         gsl_matrix_get(res, 1, 2),
 	                         gsl_matrix_get(res, 2, 2));
 
 	gsl_matrix_free(res);
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index 6e2ade6a..1b448771 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -276,23 +276,26 @@ int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec,
  * number of columns in @m, and the size of the result equals the number of rows
  * in @m.
  *
+ * The multiplication looks like this:
+ *    (a1, a2, a3) = (vec1, vec2, vec3) m
+ * Therefore matching the notation in ITA chapter 5.1.
+ *
  * Returns: a newly allocated array of signed integers containing the answer,
  * or NULL on error.
  **/
 signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec)
 {
 	signed int *ans;
-	unsigned int i;
+	unsigned int j;
 
 	ans = malloc(m->rows * sizeof(signed int));
 	if ( ans == NULL ) return NULL;
 
-	for ( i=0; i<m->rows; i++ ) {
-
-		unsigned int j;
+	for ( j=0; j<m->cols; j++ ) {
 
-		ans[i] = 0;
-		for ( j=0; j<m->cols; j++ ) {
+		unsigned int i;
+		ans[j] = 0;
+		for ( i=0; i<m->rows; i++ ) {
 			ans[i] += intmat_get(m, i, j) * vec[j];
 		}
 
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index 20a73bc8..ab1eff91 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -551,8 +551,8 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
 {
 	int i, j;
 
-	assert(ans->cols == A->cols);
-	assert(ans->rows == B->rows);
+	assert(ans->cols == B->cols);
+	assert(ans->rows == A->rows);
 	assert(A->cols == B->rows);
 
 	for ( i=0; i<ans->rows; i++ ) {
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 9ee80a16..05ee7c6b 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -197,9 +197,9 @@ static void add_symop_v(SymOpList *ops,
 	m = intmat_new(3, 3);
 	assert(m != NULL);
 
-	for ( i=0; i<3; i++ ) intmat_set(m, 0, i, h[i]);
-	for ( i=0; i<3; i++ ) intmat_set(m, 1, i, k[i]);
-	for ( i=0; i<3; i++ ) intmat_set(m, 2, i, l[i]);
+	for ( i=0; i<3; i++ ) intmat_set(m, i, 0, h[i]);
+	for ( i=0; i<3; i++ ) intmat_set(m, i, 1, k[i]);
+	for ( i=0; i<3; i++ ) intmat_set(m, i, 2, l[i]);
 
 	free(h);
 	free(k);
diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y
index bebf823f..be31c21d 100644
--- a/libcrystfel/src/symop.y
+++ b/libcrystfel/src/symop.y
@@ -101,13 +101,13 @@ symoplist:
 
 symop:
   axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[0]);
-                                     rtnl_mtx_set(m, 0, 1, $1[1]);
-                                     rtnl_mtx_set(m, 0, 2, $1[2]);
-                                     rtnl_mtx_set(m, 1, 0, $3[0]);
+                                     rtnl_mtx_set(m, 1, 0, $1[1]);
+                                     rtnl_mtx_set(m, 2, 0, $1[2]);
+                                     rtnl_mtx_set(m, 0, 1, $3[0]);
                                      rtnl_mtx_set(m, 1, 1, $3[1]);
-                                     rtnl_mtx_set(m, 1, 2, $3[2]);
-                                     rtnl_mtx_set(m, 2, 0, $5[0]);
-                                     rtnl_mtx_set(m, 2, 1, $5[1]);
+                                     rtnl_mtx_set(m, 2, 1, $3[2]);
+                                     rtnl_mtx_set(m, 0, 2, $5[0]);
+                                     rtnl_mtx_set(m, 1, 2, $5[1]);
                                      rtnl_mtx_set(m, 2, 2, $5[2]);
                                    }
 ;