diff options
| author | Thomas White <taw@physics.org> | 2023-07-07 15:02:45 +0200 |
|---|---|---|
| committer | Thomas White <taw@physics.org> | 2023-07-28 13:22:05 +0200 |
| commit | 4f7c412d650e8493e0a356f4159be650cffb4d0b (patch) | |
| tree | 23aec2e098197b121f27bbbcf5af2a4c419c1027 /libcrystfel/src/predict-refine.c | |
| parent | 1857f41d235359d9ed6c018b63ccda63edba9d2f (diff) | |
Factorise matrix operations
This makes the code much clearer. Note that two opposing sign errors
have been fixed in the gradient calculation.
Diffstat (limited to 'libcrystfel/src/predict-refine.c')
| -rw-r--r-- | libcrystfel/src/predict-refine.c | 17 |
1 files changed, 6 insertions, 11 deletions
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index fda3eaa4..00de8bb8 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -156,8 +156,7 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, gsl_matrix *M; double mu; gsl_matrix *dMdp; - gsl_matrix *minusMinvdMdp; - gsl_matrix *minusMinvdMdpMinv; + gsl_matrix *gM; /* M^-1 * dM/dx * M^-1 */ double fs, ss; get_indices(refl, &h, &k, &l); @@ -287,19 +286,15 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, } - minusMinvdMdp = gsl_matrix_calloc(3, 3); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, Minv, dMdp, 0.0, minusMinvdMdp); - minusMinvdMdpMinv = gsl_matrix_calloc(3, 3); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, minusMinvdMdp, Minv, 0.0, minusMinvdMdpMinv); - gsl_blas_dgemv(CblasNoTrans, 1.0, minusMinvdMdpMinv, t, 0.0, v); + gM = matrix_mult3(Minv, dMdp, Minv); + gsl_blas_dgemv(CblasNoTrans, -1.0, gM, t, 0.0, v); - *fsg = -(mu*gsl_vector_get(v, 1) - mu*fs*gsl_vector_get(v, 0)); - *ssg = -(mu*gsl_vector_get(v, 2) - mu*ss*gsl_vector_get(v, 0)); + *fsg = mu*(gsl_vector_get(v, 1) - fs*gsl_vector_get(v, 0)); + *ssg = mu*(gsl_vector_get(v, 2) - ss*gsl_vector_get(v, 0)); gsl_vector_free(v); - gsl_matrix_free(minusMinvdMdpMinv); + gsl_matrix_free(gM); gsl_matrix_free(dMdp); - gsl_matrix_free(minusMinvdMdp); gsl_vector_free(t); return 0; |