diff options
-rw-r--r-- | libcrystfel/CMakeLists.txt | 2 | ||||
-rw-r--r-- | libcrystfel/src/integer_matrix.c | 107 | ||||
-rw-r--r-- | libcrystfel/src/integer_matrix.h | 16 | ||||
-rw-r--r-- | libcrystfel/src/rational.c | 472 | ||||
-rw-r--r-- | libcrystfel/src/rational.h | 94 | ||||
-rw-r--r-- | tests/CMakeLists.txt | 4 | ||||
-rw-r--r-- | tests/rational_check.c | 194 |
7 files changed, 886 insertions, 3 deletions
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt index 247b925e..bd4f44ca 100644 --- a/libcrystfel/CMakeLists.txt +++ b/libcrystfel/CMakeLists.txt @@ -62,6 +62,7 @@ set(LIBCRYSTFEL_SOURCES src/peakfinder8.c src/taketwo.c src/xgandalf.c + src/rational.c ) if (HAVE_FFTW) @@ -101,6 +102,7 @@ set(LIBCRYSTFEL_HEADERS src/peakfinder8.h src/taketwo.h src/xgandalf.h + src/rational.h ) add_library(${PROJECT_NAME} SHARED diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c index 1e616e43..6e2ade6a 100644 --- a/libcrystfel/src/integer_matrix.c +++ b/libcrystfel/src/integer_matrix.c @@ -35,7 +35,9 @@ #include <string.h> #include <assert.h> +#include "rational.h" #include "integer_matrix.h" +#include "utils.h" /** @@ -95,7 +97,7 @@ IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols) * * Returns: a newly allocated copy of @m, or NULL on error/ **/ -IntegerMatrix *intmat_copy(IntegerMatrix *m) +IntegerMatrix *intmat_copy(const IntegerMatrix *m) { IntegerMatrix *p; int i, j; @@ -127,6 +129,19 @@ void intmat_free(IntegerMatrix *m) } +void intmat_size(const IntegerMatrix *m, unsigned int *rows, unsigned int *cols) +{ + if ( m == NULL ) { + *rows = 0; + *cols = 0; + return; + } + + *rows = m->rows; + *cols = m->cols; +} + + /** * intmat_set: * @m: An %IntegerMatrix @@ -162,6 +177,96 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j) } +/* intmat_intvec_mult: + * @m: An %IntegerMatrix + * @vec: An array of floats + * @ans: An array of floats in which to store the result + * + * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the + * number of columns in @m, and the size of the result equals the number of rows + * in @m. + * + * Returns: non-zero on error + **/ +int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, float *ans) +{ + unsigned int i; + + for ( i=0; i<m->rows; i++ ) { + + unsigned int j; + + ans[i] = 0.0; + for ( j=0; j<m->cols; j++ ) { + ans[i] += intmat_get(m, i, j) * vec[j]; + } + + } + + return 0; +} + + +/* intmat_rationalvec_mult: + * @m: An %IntegerMatrix + * @vec: An array of %Rational + * @ans: An array of %Rational in which to store the result + * + * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the + * number of columns in @m, and the size of the result equals the number of rows + * in @m. + * + * Returns: non-zero on error + **/ +int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec, + Rational *ans) +{ + unsigned int i; + + + for ( i=0; i<m->rows; i++ ) { + + unsigned int j; + + ans[i] = rtnl_zero(); + for ( j=0; j<m->cols; j++ ) { + Rational t; + t = rtnl_mul(vec[j], rtnl(intmat_get(m, i, j), 1)); + ans[i] = rtnl_add(ans[i], t); + } + + } + + return 0; +} + + +/* intmat_solve_rational: + * @m: An %IntegerMatrix + * @vec: An array of %Rational + * @ans: An array of %Rational in which to store the result + * + * Solves the matrix equation m*ans = vec, where @ans and @vec are + * vectors of %Rational. + * + * This is just a convenience function which creates a %RationalMatrix out of + * @m and then calls rtnl_mtx_solve(). + * + * Returns: non-zero on error + **/ +int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, + Rational *ans) +{ + RationalMatrix *rm; + int r; + + rm = rtnl_mtx_from_intmat(m); + r = rtnl_mtx_solve(rm, vec, ans); + rtnl_mtx_free(rm); + return r; +} + + /** * intmat_intvec_mult: * @m: An %IntegerMatrix diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h index c4983aef..aeea38dd 100644 --- a/libcrystfel/src/integer_matrix.h +++ b/libcrystfel/src/integer_matrix.h @@ -40,17 +40,22 @@ **/ typedef struct _integermatrix IntegerMatrix; +#include "rational.h" + #ifdef __cplusplus extern "C" { #endif /* Alloc/dealloc */ extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols); -extern IntegerMatrix *intmat_copy(IntegerMatrix *m); +extern IntegerMatrix *intmat_copy(const IntegerMatrix *m); extern IntegerMatrix *intmat_identity(int size); extern void intmat_free(IntegerMatrix *m); /* Get/set */ +extern void intmat_size(const IntegerMatrix *m, unsigned int *rows, + unsigned int *cols); + extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j, signed int v); extern signed int intmat_get(const IntegerMatrix *m, @@ -65,14 +70,21 @@ extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed i signed int m21, signed int m22, signed int m23, signed int m31, signed int m32, signed int m33); -/* Matrix-(int)vector multiplication */ +/* Matrix-vector multiplication */ +extern int intmat_floatvec_mult(const IntegerMatrix *m, const float *vec, + float *ans); extern signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec); +extern int intmat_rationalvec_mult(const IntegerMatrix *m, const Rational *vec, + Rational *ans); /* Matrix-matrix multiplication */ extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a, const IntegerMatrix *b); +extern int intmat_solve_rational(const IntegerMatrix *m, const Rational *vec, + Rational *ans); + /* Inverse */ extern IntegerMatrix *intmat_inverse(const IntegerMatrix *m); diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c new file mode 100644 index 00000000..bc783711 --- /dev/null +++ b/libcrystfel/src/rational.c @@ -0,0 +1,472 @@ +/* + * rational.c + * + * A small rational number library + * + * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2019 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <assert.h> +#include <locale.h> + +#include "rational.h" +#include "integer_matrix.h" +#include "utils.h" + + +/** + * SECTION:rational + * @short_description: Rational numbers + * @title: Rational numbers + * @section_id: + * @see_also: + * @include: "rational.h" + * @Image: + * + * A rational number library + */ + + +/* Eucliden algorithm for finding greatest common divisor */ +static signed int gcd(signed int a, signed int b) +{ + while ( b != 0 ) { + signed int t = b; + b = a % b; + a = t; + } + return a; +} + + +static void squish(Rational *rt) +{ + signed int g; + + if ( rt->num == 0 ) { + rt->den = 1; + return; + } + + g = gcd(rt->num, rt->den); + assert(g != 0); + rt->num /= g; + rt->den /= g; + + if ( rt->den < 0 ) { + rt->num = -rt->num; + rt->den = -rt->den; + } +} + + +Rational rtnl_zero() +{ + Rational r; + r.num = 0; + r.den = 1; + return r; +} + + +Rational rtnl(signed long long int num, signed long long int den) +{ + Rational r; + r.num = num; + r.den = den; + squish(&r); + return r; +} + + +double rtnl_as_double(Rational r) +{ + return (double)r.num/r.den; +} + + +static void overflow(long long int c, long long int a, long long int b) +{ + setlocale(LC_ALL, ""); + ERROR("Overflow detected in rational number library.\n"); + ERROR("%'lli < %'lli * %'lli\n", c, a, b); + abort(); +} + + +static void check_overflow(long long int c, long long int a, long long int b) +{ + if ( (a==0) || (b==0) ) { + if ( c != 0 ) overflow(c,a,b); + } else if ( (llabs(c) < llabs(a)) || (llabs(c) < llabs(b)) ) { + overflow(c,a,b); + } +} + + +Rational rtnl_mul(Rational a, Rational b) +{ + Rational r; + r.num = a.num * b.num; + r.den = a.den * b.den; + check_overflow(r.num, a.num, b.num); + check_overflow(r.den, a.den, b.den); + squish(&r); + return r; +} + + +Rational rtnl_div(Rational a, Rational b) +{ + signed int t = b.num; + b.num = b.den; + b.den = t; + return rtnl_mul(a, b); +} + + +Rational rtnl_add(Rational a, Rational b) +{ + Rational r, trt1, trt2; + + trt1.num = a.num * b.den; + trt2.num = b.num * a.den; + check_overflow(trt1.num, a.num, b.den); + check_overflow(trt2.num, b.num, a.den); + + trt1.den = a.den * b.den; + trt2.den = trt1.den; + check_overflow(trt1.den, a.den, b.den); + + r.num = trt1.num + trt2.num; + r.den = trt1.den; + squish(&r); + return r; +} + + +Rational rtnl_sub(Rational a, Rational b) +{ + b.num = -b.num; + return rtnl_add(a, b); +} + + +/* -1, 0 +1 respectively for a<b, a==b, a>b */ +signed int rtnl_cmp(Rational a, Rational b) +{ + Rational trt1, trt2; + + trt1.num = a.num * b.den; + trt2.num = b.num * a.den; + + trt1.den = a.den * b.den; + trt2.den = a.den * b.den; + + if ( trt1.num > trt2.num ) return +1; + if ( trt1.num < trt2.num ) return -1; + return 0; +} + + +Rational rtnl_abs(Rational a) +{ + Rational r = a; + squish(&r); + if ( r.num < 0 ) r.num = -r.num; + return r; +} + + +/** + * rtnl_format + * @rt: A %Rational + * + * Formats @rt as a string + * + * Returns: a string which should be freed by the caller + */ +char *rtnl_format(Rational rt) +{ + char *v = malloc(32); + if ( v == NULL ) return NULL; + if ( rt.den == 1 ) { + snprintf(v, 31, "%lli", rt.num); + } else { + snprintf(v, 31, "%lli/%lli", rt.num, rt.den); + } + return v; +} + + +/** + * SECTION:rational_matrix + * @short_description: Rational matrices + * @title: Rational matrices + * @section_id: + * @see_also: + * @include: "rational.h" + * @Image: + * + * A rational matrix library + */ + + +struct _rationalmatrix +{ + unsigned int rows; + unsigned int cols; + Rational *v; +}; + + +/** + * rtnl_mtx_new: + * @rows: Number of rows that the new matrix is to have + * @cols: Number of columns that the new matrix is to have + * + * Allocates a new %RationalMatrix with all elements set to zero. + * + * Returns: a new %RationalMatrix, or NULL on error. + **/ +RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols) +{ + RationalMatrix *m; + + m = malloc(sizeof(RationalMatrix)); + if ( m == NULL ) return NULL; + + m->v = calloc(rows*cols, sizeof(Rational)); + if ( m->v == NULL ) { + free(m); + return NULL; + } + + m->rows = rows; + m->cols = cols; + + return m; +} + + +RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m) +{ + RationalMatrix *n; + int i; + + n = rtnl_mtx_new(m->rows, m->cols); + if ( n == NULL ) return NULL; + + for ( i=0; i<m->rows*m->cols; i++ ) { + n->v[i] = m->v[i]; + } + + return n; +} + + +Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j) +{ + assert(m != NULL); + return m->v[j+m->cols*i]; +} + + +void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v) +{ + assert(m != NULL); + m->v[j+m->cols*i] = v; +} + + +RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m) +{ + RationalMatrix *n; + unsigned int rows, cols; + int i, j; + + intmat_size(m, &rows, &cols); + n = rtnl_mtx_new(rows, cols); + if ( n == NULL ) return NULL; + + for ( i=0; i<rows; i++ ) { + for ( j=0; j<cols; j++ ) { + n->v[j+cols*i] = rtnl(intmat_get(m, i, j), 1); + } + } + + return n; +} + + +void rtnl_mtx_free(RationalMatrix *mtx) +{ + if ( mtx == NULL ) return; + free(mtx->v); + free(mtx); +} + + +/* rtnl_mtx_solve: + * @m: A %RationalMatrix + * @vec: An array of %Rational + * @ans: An array of %Rational in which to store the result + * + * Solves the matrix equation m*ans = vec, where @ans and @vec are + * vectors of %Rational. + * + * In this version, @m must be square. + * + * The number of columns in @m must equal the length of @ans, and the number + * of rows in @m must equal the length of @vec, but note that there is no way + * for this function to check that this is the case. + * + * Returns: non-zero on error + **/ +int rtnl_mtx_solve(const RationalMatrix *m, const Rational *ivec, Rational *ans) +{ + RationalMatrix *cm; + Rational *vec; + int h, k; + int i; + + if ( m->rows != m->cols ) return 1; + + /* Copy the matrix and vector because the calculation will + * be done in-place */ + cm = rtnl_mtx_copy(m); + if ( cm == NULL ) return 1; + + vec = malloc(m->rows*sizeof(Rational)); + if ( vec == NULL ) return 1; + for ( h=0; h<m->rows; h++ ) vec[h] = ivec[h]; + + /* Gaussian elimination with partial pivoting */ + h = 0; + k = 0; + while ( h<=m->rows && k<=m->cols ) { + + int prow = 0; + Rational pval = rtnl_zero(); + Rational t; + + /* Find the row with the largest value in column k */ + for ( i=h; i<m->rows; i++ ) { + if ( rtnl_cmp(rtnl_abs(rtnl_mtx_get(cm, i, k)), pval) > 0 ) { + pval = rtnl_abs(rtnl_mtx_get(cm, i, k)); + prow = i; + } + } + + if ( rtnl_cmp(pval, rtnl_zero()) == 0 ) { + k++; + continue; + } + + /* Swap 'prow' with row h */ + for ( i=0; i<m->cols; i++ ) { + t = rtnl_mtx_get(cm, h, i); + rtnl_mtx_set(cm, h, i, rtnl_mtx_get(cm, prow, i)); + rtnl_mtx_set(cm, prow, i, t); + } + t = vec[prow]; + vec[prow] = vec[h]; + vec[h] = t; + + /* Divide and subtract rows below */ + for ( i=h+1; i<m->rows; i++ ) { + + int j; + Rational dval; + + dval = rtnl_div(rtnl_mtx_get(cm, i, k), + rtnl_mtx_get(cm, h, k)); + + for ( j=0; j<m->cols; j++ ) { + Rational t = rtnl_mtx_get(cm, i, j); + Rational p = rtnl_mul(dval, rtnl_mtx_get(cm, h, j)); + t = rtnl_sub(t, p); + rtnl_mtx_set(cm, i, j, t); + } + + /* Divide the right hand side as well */ + Rational t = vec[i]; + Rational p = rtnl_mul(dval, vec[h]); + vec[i] = rtnl_sub(t, p); + } + + h++; + k++; + + + } + + /* Back-substitution */ + for ( i=m->rows-1; i>=0; i-- ) { + int j; + Rational sum = rtnl_zero(); + for ( j=i+1; j<m->cols; j++ ) { + Rational av; + av = rtnl_mul(rtnl_mtx_get(cm, i, j), ans[j]); + sum = rtnl_add(sum, av); + } + sum = rtnl_sub(vec[i], sum); + ans[i] = rtnl_div(sum, rtnl_mtx_get(cm, i, i)); + } + + free(vec); + rtnl_mtx_free(cm); + + return 0; +} + + +/** + * rtnl_mtx_print + * @m: A %RationalMatrix + * + * Prints @m to stderr. + * + */ +void rtnl_mtx_print(const RationalMatrix *m) +{ + unsigned int i, j; + + for ( i=0; i<m->rows; i++ ) { + + fprintf(stderr, "[ "); + for ( j=0; j<m->cols; j++ ) { + char *v = rtnl_format(m->v[j+m->cols*i]); + fprintf(stderr, "%4s ", v); + free(v); + } + fprintf(stderr, "]\n"); + } +} diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h new file mode 100644 index 00000000..8b182163 --- /dev/null +++ b/libcrystfel/src/rational.h @@ -0,0 +1,94 @@ +/* + * rational.h + * + * A small rational number library + * + * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2019 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef RATIONAL_H +#define RATIONAL_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +/** + * Rational + * + * The Rational is an opaque-ish data structure representing a rational number. + * + * "Opaque-ish" means that the structure isn't technically opaque, allowing you + * to assign and allocate them easily. But you shouldn't look at or set its + * contents, except by using the accessor functions. + **/ +typedef struct { + /* Private, don't modify */ + signed long long int num; + signed long long int den; +} Rational; + + +/** + * RationalMatrix + * + * The RationalMatrix is an opaque data structure representing a matrix of + * rational numbers. + **/ +typedef struct _rationalmatrix RationalMatrix; + +#include "integer_matrix.h" + +#ifdef __cplusplus +extern "C" { +#endif + +extern Rational rtnl_zero(void); +extern Rational rtnl(signed long long int num, signed long long int den); +extern double rtnl_as_double(Rational r); + +extern Rational rtnl_mul(Rational a, Rational b); +extern Rational rtnl_div(Rational a, Rational b); +extern Rational rtnl_add(Rational a, Rational b); +extern Rational rtnl_sub(Rational a, Rational b); + +extern signed int rtnl_cmp(Rational a, Rational b); +extern Rational rtnl_abs(Rational a); + +extern char *rtnl_format(Rational rt); + +extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols); +extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m); +extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j); +extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v); +extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m); +extern void rtnl_mtx_free(RationalMatrix *mtx); +extern int rtnl_mtx_solve(const RationalMatrix *m, const Rational *vec, + Rational *ans); +extern void rtnl_mtx_print(const RationalMatrix *m); + +#ifdef __cplusplus +} +#endif + +#endif /* RATIONAL_H */ diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 02f3ca02..beb141c3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -83,3 +83,7 @@ if (HAVE_OPENCL) add_test(gpu_sim_check gpu_sim_check) endif (HAVE_OPENCL) +add_executable(rational_check rational_check.c) +target_include_directories(rational_check PRIVATE ${COMMON_INCLUDES}) +target_link_libraries(rational_check ${COMMON_LIBRARIES}) +add_test(rational_check rational_check) diff --git a/tests/rational_check.c b/tests/rational_check.c new file mode 100644 index 00000000..8957cfdd --- /dev/null +++ b/tests/rational_check.c @@ -0,0 +1,194 @@ +/* + * rational_check.c + * + * Check that rational numbers work + * + * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2019 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdlib.h> +#include <stdio.h> +#include <stdarg.h> +#include <gsl/gsl_rng.h> + +#include <rational.h> +#include <utils.h> + + +static Rational gen_rtnl(gsl_rng *rng, double *dptr, int sz) +{ + Rational r; + signed int num, den; + do { + num = gsl_rng_uniform_int(rng, 2*sz) - sz; + den = gsl_rng_uniform_int(rng, 2*sz) - sz; + } while ( den == 0 ); + r = rtnl(num, den); + *dptr = (double)num/den; + return r; +} + + +static int test_rational(gsl_rng *rng) +{ + Rational r1, r2, r3; + double rd1, rd2, rd3; + char op; + + r1 = gen_rtnl(rng, &rd1, 100); + r2 = gen_rtnl(rng, &rd2, 100); + + switch ( gsl_rng_uniform_int(rng, 3) ) { + + case 0: + r3 = rtnl_mul(r1, r2); + rd3 = rd1*rd2; + op = '*'; + break; + + case 1: + r3 = rtnl_div(r1, r2); + rd3 = rd1/rd2; + op = '/'; + break; + + case 2: + r3 = rtnl_add(r1, r2); + rd3 = rd1+rd2; + op = '+'; + break; + + case 3: + r3 = rtnl_sub(r1, r2); + rd3 = rd1-rd2; + op = '-'; + break; + + default: + abort(); + + } + + if ( isinf(rd3) && isinf(rtnl_as_double(r3)) ) return 0; + + if ( fabs(rtnl_as_double(r3) - rd3) > 0.001 ) { + ERROR("%10s %c %10s = %10s", rtnl_format(r1), op, + rtnl_format(r2), rtnl_format(r3)); + ERROR(" ---> %10f %c %10f = %10f\n", rd1, op, rd2, rd3); + return 1; + } + + return 0; +} + + +static int test_rational_matrix(gsl_rng *rng) +{ + const int size = 3; + RationalMatrix *rm; + Rational *rvec; + Rational *rans; + gsl_matrix *gm; + gsl_vector *gvec; + gsl_vector *gans; + int i, j; + int filt; + int fail = 0; + + rm = rtnl_mtx_new(size, size); + gm = gsl_matrix_alloc(size, size); + rvec = malloc(size*sizeof(Rational)); + gvec = gsl_vector_alloc(size); + rans = malloc(size*sizeof(Rational)); + if ( (rm==NULL) || (gm==NULL) || (rvec==NULL) + || (gvec==NULL) || (rans==NULL) ) return 1; + + /* Find a matrix equation which actually works */ + do { + for ( i=0; i<size; i++ ) { + double d; + Rational r; + for ( j=0; j<size; j++ ) { + r = gen_rtnl(rng, &d, 5); + rtnl_mtx_set(rm, i, j, r); + gsl_matrix_set(gm ,i, j, d); + } + r = gen_rtnl(rng, &d, 5); + rvec[i] = r; + gsl_vector_set(gvec ,i, d); + } + + gans = solve_svd(gvec, gm, &filt, 0); + } while ( (gans==NULL) || (filt!=0) || (isnan(gsl_vector_get(gans,0))) ); + + rtnl_mtx_solve(rm, rvec, rans); + + for ( i=0; i<size; i++ ) { + if ( fabs(rtnl_as_double(rans[i]) - gsl_vector_get(gans, i) > 0.001) ) + { + STATUS("%25s %10f %10f\n", rtnl_format(rans[i]), + rtnl_as_double(rans[i]), + gsl_vector_get(gans, i)); + fail = 1; + } + } + + gsl_matrix_free(gm); + rtnl_mtx_free(rm); + free(rans); + free(rvec); + gsl_vector_free(gvec); + gsl_vector_free(gans); + return fail; +} + + +int main(int argc, char *argv[]) +{ + int fail = 0; + gsl_rng *rng; + int i; + + rng = gsl_rng_alloc(gsl_rng_mt19937); + gsl_set_error_handler_off(); + + STATUS("Arithmetic test...\n"); + for ( i=0; i<1000; i++ ) { + fail += test_rational(rng); + } + + STATUS("Matrix test...\n"); + for ( i=0; i<1000; i++ ) { + fail += test_rational_matrix(rng); + } + + gsl_rng_free(rng); + + if ( fail != 0 ) return 1; + return 0; +} |