Add solve_inv() for simpler matrix equations

2 files changed, 20 insertions, 0 deletions

diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index 5204e5f1..0ff09b58 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -167,6 +167,25 @@ gsl_matrix *matrix_invert(gsl_matrix *m)
 }
 
 
+gsl_vector *solve_inv(gsl_vector *v, gsl_matrix *M)
+{
+	gsl_matrix *Minv;
+	gsl_vector *s;
+	int n;
+
+	n = v->size;
+	if ( v->size != M->size1 ) return NULL;
+	if ( v->size != M->size2 ) return NULL;
+
+	Minv =  matrix_invert(M);
+	if ( Minv == NULL ) return NULL;
+
+	s = gsl_vector_calloc(n);
+	gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, s);
+	return s;
+}
+
+
 static int check_eigen(gsl_vector *e_val, int verbose)
 {
 	int i;
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index 30a7c4ef..08e059a7 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -75,6 +75,7 @@ extern "C" {
 extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v);
 extern void show_matrix(gsl_matrix *M);
 extern void show_vector(gsl_vector *M);
+extern gsl_vector *solve_inv(gsl_vector *v, gsl_matrix *M);
 extern gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
                             int verbose);
 extern gsl_matrix *matrix_mult2(gsl_matrix *A, gsl_matrix *B);