Avoid converting IntegerMatrix to RationalMatrix

3 files changed, 46 insertions, 12 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 6e7cebb5..dbc845bc 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -1844,6 +1844,23 @@ static void rtnl_mult_in_place(RationalMatrix *T, RationalMatrix *M)
 }
 
 
+/* Replace the elements of T with those of T*M */
+static void rtnl_int_mult_in_place(RationalMatrix *T, IntegerMatrix *M)
+{
+	int i, j;
+	RationalMatrix *tmp;
+
+	tmp = rtnl_mtx_new(3, 3);
+	rtnl_mtx_intmatmult(T, M, tmp);
+	for ( i=0; i<3; i++ ) {
+		for ( j=0; j<3; j++ ) {
+			rtnl_mtx_set(T, i, j, rtnl_mtx_get(tmp, i, j));
+		}
+	}
+	rtnl_mtx_free(tmp);
+}
+
+
 /* Cell volume from G6 components */
 static double g6_volume(struct g6 g)
 {
@@ -2243,24 +2260,16 @@ UnitCell *compare_reindexed_cell_parameters(UnitCell *cell_in,
 		//intmat_print(P);
 
 		/* Calculate combined matrix: CB.RiB.RA.CiA */
-		RationalMatrix *rRA = rtnl_mtx_from_intmat(RA);
-		RationalMatrix *rP = rtnl_mtx_from_intmat(P);
 		IntegerMatrix *RiB = intmat_inverse(RB);
-		RationalMatrix *rRiB = rtnl_mtx_from_intmat(RiB);
-		RationalMatrix *rCB = rtnl_mtx_from_intmat(CB);
 		RationalMatrix *comb = rtnl_mtx_identity(3);
 
 		rtnl_mult_in_place(comb, CiA);
-		rtnl_mult_in_place(comb, rRA);
-		rtnl_mult_in_place(comb, rP);
-		rtnl_mult_in_place(comb, rRiB);
-		rtnl_mult_in_place(comb, rCB);
+		rtnl_int_mult_in_place(comb, RA);
+		rtnl_int_mult_in_place(comb, P);
+		rtnl_int_mult_in_place(comb, RiB);
+		rtnl_int_mult_in_place(comb, CB);
 
-		rtnl_mtx_free(rRA);
-		rtnl_mtx_free(rP);
 		intmat_free(RiB);
-		rtnl_mtx_free(rRiB);
-		rtnl_mtx_free(rCB);
 
 		match = cell_transform_rational(cell_in, comb);
 		//STATUS("Original cell transformed to look like reference:\n");
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index ce286d7a..f8d3cb68 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -553,6 +553,29 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
 }
 
 
+void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B,
+                         RationalMatrix *ans)
+{
+	int i, j;
+
+	assert(ans->rows == A->rows);
+
+	for ( i=0; i<ans->rows; i++ ) {
+		for ( j=0; j<ans->cols; j++ ) {
+			int k;
+			Rational sum = rtnl_zero();
+			for ( k=0; k<A->rows; k++ ) {
+				Rational add;
+				add = rtnl_mul(rtnl_mtx_get(A, i, k),
+				               rtnl(intmat_get(B, k, j), 1));
+				sum = rtnl_add(sum, add);
+			}
+			rtnl_mtx_set(ans, i, j, sum);
+		}
+	}
+}
+
+
 /* Given a "P-matrix" (see ITA chapter 5.1), calculate the fractional
  * coordinates of point "vec" in the original axes, given its fractional
  * coordinates in the transformed axes.
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
index 8681bb06..012b929e 100644
--- a/libcrystfel/src/rational.h
+++ b/libcrystfel/src/rational.h
@@ -92,6 +92,8 @@ extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m);
 extern void rtnl_mtx_free(RationalMatrix *mtx);
 extern void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
                              RationalMatrix *ans);
+extern void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B,
+                                RationalMatrix *ans);
 extern int transform_fractional_coords_rtnl(const RationalMatrix *P,
                                             const Rational *ivec,
                                             Rational *ans);