From 24e2744b8b53dde715f59d98f969d44a10530675 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Wed, 28 Aug 2019 15:33:46 +0200
Subject: Rationalise matrix mutliplication

compare_reindexed_cell_parameters still needs to be updated
---
 libcrystfel/src/cell-utils.c     | 46 ++++--------------------------
 libcrystfel/src/integer_matrix.c |  4 +--
 libcrystfel/src/integer_matrix.h |  4 +--
 libcrystfel/src/rational.c       | 61 +++++++++++++++++++++++++++++++++-------
 libcrystfel/src/rational.h       | 11 +++++---
 libcrystfel/src/symmetry.c       | 10 +++----
 6 files changed, 72 insertions(+), 64 deletions(-)

(limited to 'libcrystfel')

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index dbc845bc..de8fe7ae 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -1612,7 +1612,6 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in
 	Rational *cand_c;
 	int ncand_a, ncand_b, ncand_c;
 	int ia, ib;
-	RationalMatrix *MCB;
 	RationalMatrix *CiAMCB;
 	double min_dist = +INFINITY;
 
@@ -1647,8 +1646,6 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in
 	}
 
 	M = rtnl_mtx_new(3, 3);
-	MCB = rtnl_mtx_new(3, 3);
-	CiAMCB = rtnl_mtx_new(3, 3);
 	for ( ia=0; ia<ncand_a; ia++ ) {
 		for ( ib=0; ib<ncand_b; ib++ ) {
 
@@ -1656,6 +1653,7 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in
 			double at, bt, ct, alt, bet, gat;
 			double dist;
 			int ic = 0;
+			RationalMatrix *MCB;
 
 			/* Form the matrix using the first candidate for c */
 			rtnl_mtx_set(M, 0, 0, cand_a[3*ia+0]);
@@ -1710,8 +1708,9 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in
 				dist = g6_distance(test, reference);
 				if ( dist < min_dist ) {
 					min_dist = dist;
-					rtnl_mtx_mtxmult(M, CB, MCB);
-					rtnl_mtx_mtxmult(CiA, MCB, CiAMCB);
+					MCB = rtnlmtx_times_rtnlmtx(M, CB);
+					CiAMCB = rtnlmtx_times_rtnlmtx(CiA, MCB);
+					rtnl_mtx_free(MCB);
 				}
 
 				cell_free(test);
@@ -1721,7 +1720,6 @@ int compare_derivative_cell_parameters(UnitCell *cell_in, UnitCell *reference_in
 	}
 
 	rtnl_mtx_free(M);
-	rtnl_mtx_free(MCB);
 	free(cand_a);
 	free(cand_b);
 	free(cand_c);
@@ -1817,7 +1815,7 @@ static void debug_lattice(struct g6 g, double eps, int step)
 static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M)
 {
 	int i, j;
-	IntegerMatrix *tmp = intmat_intmat_mult(T, M);
+	IntegerMatrix *tmp = intmat_times_intmat(T, M);
 	for ( i=0; i<3; i++ ) {
 		for ( j=0; j<3; j++ ) {
 			intmat_set(T, i, j, intmat_get(tmp, i, j));
@@ -1827,40 +1825,6 @@ static void mult_in_place(IntegerMatrix *T, IntegerMatrix *M)
 }
 
 
-/* Replace the elements of T with those of T*M */
-static void rtnl_mult_in_place(RationalMatrix *T, RationalMatrix *M)
-{
-	int i, j;
-	RationalMatrix *tmp;
-
-	tmp = rtnl_mtx_new(3, 3);
-	rtnl_mtx_mtxmult(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);
-}
-
-
-/* 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)
 {
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index c633a620..16aff5bc 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -218,8 +218,8 @@ signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl)
  * \returns A newly allocated \ref IntegerMatrix containing the answer, or NULL on
  * error.
  **/
-IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
-                                  const IntegerMatrix *b)
+IntegerMatrix *intmat_times_intmat(const IntegerMatrix *a,
+                                   const IntegerMatrix *b)
 {
 	unsigned int i, j;
 	IntegerMatrix *ans;
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
index 5a365f0a..25c0fb99 100644
--- a/libcrystfel/src/integer_matrix.h
+++ b/libcrystfel/src/integer_matrix.h
@@ -74,8 +74,8 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i
 extern signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl);
 
 /* Matrix-matrix multiplication */
-extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
-                                         const IntegerMatrix *b);
+extern IntegerMatrix *intmat_times_intmat(const IntegerMatrix *a,
+                                          const IntegerMatrix *b);
 
 /* Inverse */
 extern IntegerMatrix *intmat_inverse(const IntegerMatrix *m);
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
index f8d3cb68..d59da59c 100644
--- a/libcrystfel/src/rational.c
+++ b/libcrystfel/src/rational.c
@@ -528,14 +528,18 @@ void rtnl_mtx_print(const RationalMatrix *m)
 }
 
 
-void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
-                      RationalMatrix *ans)
+RationalMatrix *rtnlmtx_times_intmat(const RationalMatrix *A,
+                                     const IntegerMatrix *B)
 {
 	int i, j;
+	RationalMatrix *ans;
+	unsigned int B_rows, B_cols;
 
-	assert(ans->cols == B->cols);
-	assert(ans->rows == A->rows);
-	assert(A->cols == B->rows);
+	intmat_size(B, &B_rows, &B_cols);
+	assert(A->cols == B_rows);
+
+	ans = rtnl_mtx_new(A->rows, B_cols);
+	if ( ans == NULL ) return NULL;
 
 	for ( i=0; i<ans->rows; i++ ) {
 		for ( j=0; j<ans->cols; j++ ) {
@@ -544,21 +548,27 @@ void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
 			for ( k=0; k<A->rows; k++ ) {
 				Rational add;
 				add = rtnl_mul(rtnl_mtx_get(A, i, k),
-				               rtnl_mtx_get(B, k, j));
+				               rtnl(intmat_get(B, k, j), 1));
 				sum = rtnl_add(sum, add);
 			}
 			rtnl_mtx_set(ans, i, j, sum);
 		}
 	}
+
+	return ans;
 }
 
 
-void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B,
-                         RationalMatrix *ans)
+RationalMatrix *rtnlmtx_times_rtnlmtx(const RationalMatrix *A,
+                                      const RationalMatrix *B)
 {
 	int i, j;
+	RationalMatrix *ans;
+
+	assert(A->cols == B->rows);
 
-	assert(ans->rows == A->rows);
+	ans = rtnl_mtx_new(A->rows, B->cols);
+	if ( ans == NULL ) return NULL;
 
 	for ( i=0; i<ans->rows; i++ ) {
 		for ( j=0; j<ans->cols; j++ ) {
@@ -567,12 +577,43 @@ void rtnl_mtx_intmatmult(const RationalMatrix *A, const IntegerMatrix *B,
 			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));
+				               rtnl_mtx_get(B, k, j));
+				sum = rtnl_add(sum, add);
+			}
+			rtnl_mtx_set(ans, i, j, sum);
+		}
+	}
+
+	return ans;
+}
+
+
+RationalMatrix *intmat_times_rtnlmtx(const IntegerMatrix *A, const RationalMatrix *B)
+{
+	int i, j;
+	RationalMatrix *ans;
+	unsigned int A_rows, A_cols;
+
+	intmat_size(A, &A_rows, &A_cols);
+
+	ans = rtnl_mtx_new(A_rows, B->cols);
+	if ( ans == NULL ) return NULL;
+
+	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(intmat_get(A, i, k), 1),
+				               rtnl_mtx_get(B, k, j));
 				sum = rtnl_add(sum, add);
 			}
 			rtnl_mtx_set(ans, i, j, sum);
 		}
 	}
+
+	return ans;
 }
 
 
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
index 012b929e..880bdbeb 100644
--- a/libcrystfel/src/rational.h
+++ b/libcrystfel/src/rational.h
@@ -90,10 +90,13 @@ extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m);
 extern RationalMatrix *rtnl_mtx_identity(int rows);
 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 RationalMatrix *rtnlmtx_times_rtnlmtx(const RationalMatrix *A,
+                                             const RationalMatrix *B);
+extern RationalMatrix *rtnlmtx_times_intmat(const RationalMatrix *A,
+                                            const IntegerMatrix *B);
+extern RationalMatrix *intmat_times_rtnlmtx(const IntegerMatrix *a,
+                                            const RationalMatrix *b);
+
 extern int transform_fractional_coords_rtnl(const RationalMatrix *P,
                                             const Rational *ivec,
                                             Rational *ans);
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 6739daaf..e834e357 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -326,7 +326,7 @@ static void expand_ops(SymOpList *s)
 				int k, nk;
 				int found;
 
-				m = intmat_intmat_mult(opi, opj);
+				m = intmat_times_intmat(opi, opj);
 				assert(m != NULL);
 
 				nk = num_ops(s);
@@ -380,13 +380,13 @@ static void transform_ops(SymOpList *s, IntegerMatrix *P)
 
 		IntegerMatrix *r, *f;
 
-		r = intmat_intmat_mult(P, s->ops[i]);
+		r = intmat_times_intmat(P, s->ops[i]);
 		if ( r == NULL ) {
 			ERROR("Matrix multiplication failed.\n");
 			return;
 		}
 
-		f = intmat_intmat_mult(r, Pi);
+		f = intmat_times_intmat(r, Pi);
 		if ( f == NULL ) {
 			ERROR("Matrix multiplication failed.\n");
 			return;
@@ -1340,7 +1340,7 @@ static int check_mult(const IntegerMatrix *ans,
 	int val;
 	IntegerMatrix *m;
 
-	m = intmat_intmat_mult(a, b);
+	m = intmat_times_intmat(a, b);
 	assert(m != NULL);
 
 	val = intmat_equals(ans, m);
@@ -1403,7 +1403,7 @@ static int order(const IntegerMatrix *m)
 
 		IntegerMatrix *anew;
 
-		anew = intmat_intmat_mult(m, a);
+		anew = intmat_times_intmat(m, a);
 		assert(anew != NULL);
 		intmat_free(a);
 		a = anew;
-- 
cgit v1.2.3