Commit 254e11e5 authored by Thomas White's avatar Thomas White
Browse files

Check lots of GSL error codes

parent 1bd591d2
......@@ -231,6 +231,11 @@ static UnitCell *cell_new_from_axes(struct rvec as, struct rvec bs,
*/
m = gsl_matrix_alloc(3, 3);
if ( m == NULL ) {
ERROR("Couldn't allocate memory for matrix\n");
free(cell);
return NULL;
}
gsl_matrix_set(m, 0, 0, as.u);
gsl_matrix_set(m, 0, 1, as.v);
gsl_matrix_set(m, 0, 2, as.w);
......@@ -243,9 +248,34 @@ static UnitCell *cell_new_from_axes(struct rvec as, struct rvec bs,
/* Invert */
perm = gsl_permutation_alloc(m->size1);
if ( perm == NULL ) {
ERROR("Couldn't allocate permutation\n");
free(cell);
gsl_matrix_free(m);
return NULL;
}
inv = gsl_matrix_alloc(m->size1, m->size2);
gsl_linalg_LU_decomp(m, perm, &s);
gsl_linalg_LU_invert(m, perm, inv);
if ( inv == NULL ) {
ERROR("Couldn't allocate inverse\n");
gsl_matrix_free(m);
gsl_permutation_free(perm);
free(cell);
return NULL;
}
if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
ERROR("Couldn't decompose matrix");
gsl_matrix_free(m);
gsl_permutation_free(perm);
free(cell);
return NULL;
}
if ( gsl_linalg_LU_invert(m, perm, inv) ) {
ERROR("Couldn't invert matrix");
gsl_matrix_free(m);
gsl_permutation_free(perm);
free(cell);
return NULL;
}
gsl_permutation_free(perm);
gsl_matrix_free(m);
......@@ -294,7 +324,7 @@ void cell_get_cartesian(UnitCell *cell,
}
void cell_get_reciprocal(UnitCell *cell,
int cell_get_reciprocal(UnitCell *cell,
double *asx, double *asy, double *asz,
double *bsx, double *bsy, double *bsz,
double *csx, double *csy, double *csz)
......@@ -305,6 +335,11 @@ void cell_get_reciprocal(UnitCell *cell,
gsl_permutation *perm;
m = gsl_matrix_alloc(3, 3);
if ( m == NULL ) {
ERROR("Couldn't allocate memory for matrix\n");
free(cell);
return -1;
}
gsl_matrix_set(m, 0, 0, cell->ax);
gsl_matrix_set(m, 0, 1, cell->bx);
gsl_matrix_set(m, 0, 2, cell->cx);
......@@ -315,13 +350,37 @@ void cell_get_reciprocal(UnitCell *cell,
gsl_matrix_set(m, 2, 1, cell->bz);
gsl_matrix_set(m, 2, 2, cell->cz);
/* Invert */
/* Invert */
perm = gsl_permutation_alloc(m->size1);
if ( perm == NULL ) {
ERROR("Couldn't allocate permutation\n");
free(cell);
gsl_matrix_free(m);
return -1;
}
inv = gsl_matrix_alloc(m->size1, m->size2);
gsl_linalg_LU_decomp(m, perm, &s);
gsl_linalg_LU_invert(m, perm, inv);
gsl_permutation_free(perm);
gsl_matrix_free(m);
if ( inv == NULL ) {
ERROR("Couldn't allocate inverse\n");
gsl_matrix_free(m);
gsl_permutation_free(perm);
free(cell);
return -1;
}
if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
ERROR("Couldn't decompose matrix\n");
gsl_matrix_free(m);
gsl_permutation_free(perm);
free(cell);
return -1;
}
if ( gsl_linalg_LU_invert(m, perm, inv) ) {
ERROR("Couldn't invert matrix\n");
gsl_matrix_free(m);
gsl_permutation_free(perm);
free(cell);
return -1;
}
/* Transpose */
gsl_matrix_transpose(inv);
......@@ -337,6 +396,8 @@ void cell_get_reciprocal(UnitCell *cell,
*csz = gsl_matrix_get(inv, 2, 2);
gsl_matrix_free(inv);
return 0;
}
......@@ -430,9 +491,12 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
"----------------------------\n");
}
cell_get_reciprocal(template, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
if ( cell_get_reciprocal(template, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz) ) {
ERROR("Couldn't get reciprocal cell for template.\n");
return NULL;
}
lengths[0] = modulus(asx, asy, asz);
lengths[1] = modulus(bsx, bsy, bsz);
......@@ -446,9 +510,12 @@ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose)
cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz) ) {
ERROR("Couldn't get reciprocal cell.\n");
return NULL;
}
/* Negative values mean 1/n, positive means n, zero means zero */
for ( n1l=-2; n1l<=4; n1l++ ) {
......
......@@ -65,7 +65,7 @@ extern void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az
extern void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz);
extern void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz);
extern void cell_get_reciprocal(UnitCell *cell,
extern int cell_get_reciprocal(UnitCell *cell,
double *asx, double *asy, double *asz,
double *bsx, double *bsy, double *bsz,
double *csx, double *csy, double *csz);
......
......@@ -21,6 +21,7 @@
#include <unistd.h>
#include <getopt.h>
#include <hdf5.h>
#include <gsl/gsl_errno.h>
#include "utils.h"
#include "hdf5-file.h"
......@@ -318,6 +319,8 @@ int main(int argc, char *argv[])
return 1;
}
gsl_set_error_handler_off();
n_images = 0;
n_hits = 0;
do {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment