diff --git a/contrib/kbipack/Makefile~ b/contrib/kbipack/Makefile~ deleted file mode 100755 index 3c361030c38340dcf7d3884a72c814f6c9f1ad52..0000000000000000000000000000000000000000 --- a/contrib/kbipack/Makefile~ +++ /dev/null @@ -1,75 +0,0 @@ -#FLAGS = -O3 -funroll-all-loops -ansi -pedantic -Wall -U DEBUG -# Flags for different diagnostics: -#FLAGS = -g -ansi -pedantic -Wall -fprofile-arcs -ftest-coverage -pg -#FLAGS = -g -ansi -pedantic -Wall -pg - -# If gmp.h is not found de facto (e.g. when it is privately intalled), -# its path should be put here -#INCL = - -#LIBS = -lgmp -#LIBDIR = . -#SRC = gmp_normal_form.c gmp_matrix_io.c gmp_matrix.c gmp_blas.c - -# The compiler options are for gcc. If it is not available, -# suitable flags must be set for the compiler -#CC = gcc -#RM = rm -f -#RANLIB = ranlib - -#libkbi.a: $(SRC) -# $(RM) *.o -# $(CC) -c $(FLAGS) $(SRC) -# ar r libkbi.a *.o -# $(RM) *.o -# $(RANLIB) libkbi.a - -#compute_normal_form: libkbi compute_normal_form.c -# $(CC) $(FLAGS) -L $(LIBDIR) $(INCL) compute_normal_form.c -o compute_normal_form -lkbi $(LIBS) -#clean: -# $(RM) compute_normal_form -# $(RM) *.o - -## - -include ../../variables - -LIB = ../../lib/libGmshKbi${LIBEXT} - -INC = ${DASH}I. - -CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE} -lgmp - -SRC = gmp_normal_form.c gmp_matrix_io.c gmp_matrix.c gmp_blas.c - -OBJ = ${SRC:.c=${OBJEXT}} - -.SUFFIXES: ${OBJEXT} .c - -${LIB}: ${OBJ} - ${AR} ${ARFLAGS}${LIB} ${OBJ} - ${RANLIB} ${LIB} - -cpobj: ${OBJ} - cp -f ${OBJ} ../../lib/ - -.c${OBJEXT}: - ${CC} ${CFLAGS} ${DASH}c $< - -clean: - ${RM} *.o *.obj - -depend: - (sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \ - ${CXX} -MM ${CFLAGS} ${SRC} | sed 's/.o:/$${OBJEXT}:/g' \ - ) > Makefile.new - cp Makefile Makefile.bak - cp Makefile.new Makefile - rm -f Makefile.new - -# DO NOT DELETE THIS LINE -gmp_normal_form${OBJEXT}: gmp_normal_form.c gmp_blas.h gmp_matrix.h \ - gmp_normal_form.h -gmp_matrix_io${OBJEXT}: gmp_matrix_io.c gmp_matrix_io.h gmp_matrix.h gmp_blas.h -gmp_matrix${OBJEXT}: gmp_matrix.c gmp_matrix.h gmp_blas.h -gmp_blas${OBJEXT}: gmp_blas.c gmp_blas.h diff --git a/contrib/kbipack/gmp_blas.o b/contrib/kbipack/gmp_blas.o deleted file mode 100755 index 02064581725915c72fb204793a3f9436e14c5cfc..0000000000000000000000000000000000000000 Binary files a/contrib/kbipack/gmp_blas.o and /dev/null differ diff --git a/contrib/kbipack/gmp_matrix.c~ b/contrib/kbipack/gmp_matrix.c~ deleted file mode 100755 index 2afe8d1dd2d6b34312a4859e01d672a1480bf947..0000000000000000000000000000000000000000 --- a/contrib/kbipack/gmp_matrix.c~ +++ /dev/null @@ -1,678 +0,0 @@ -/* - Source file for integer-oriented matrix, relying on the arbitrary - precision integers from the GNU gmp package. - - Copyright (C) 28.10.2003 Saku Suuriniemi TUT/CEM - - This program 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 2 of the License, or - any later version. - - This program 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 this program; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - Saku Suuriniemi, TUT/Electromagetics - P.O.Box 692, FIN-33101 Tampere, Finland - saku.suuriniemi@tut.fi - - $Id: gmp_matrix.c~,v 1.1 2009-03-30 14:10:57 matti Exp $ -*/ - - -#include<stdlib.h> -#include"gmp_matrix.h" - -int gmp_kaka(int laalaa){ - printf("kissa %d", laalaa); - return laalaa; -} - -gmp_matrix * -create_gmp_matrix(size_t r, size_t c, - const mpz_t * e) -{ - gmp_matrix * new_matrix; - size_t ind; - - new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix)); - if(new_matrix == NULL) - { - return NULL; - } - - new_matrix -> storage = (mpz_t *) calloc(r*c, sizeof(mpz_t)); - if(new_matrix -> storage == NULL) - { - free(new_matrix); - return NULL; - } - - new_matrix -> rows = r; - new_matrix -> cols = c; - - for(ind = 0; ind < r*c; ind ++) - { - mpz_init (new_matrix -> storage[ind]); - mpz_set (new_matrix -> storage[ind], e[ind]); - } - - return new_matrix; -} - - -gmp_matrix * -create_gmp_matrix_identity(size_t dim) -{ - gmp_matrix * new_matrix; - size_t ind; - - new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix)); - if(new_matrix == NULL) - { - return NULL; - } - - new_matrix -> storage = (mpz_t *) calloc(dim*dim, sizeof(mpz_t)); - if(new_matrix -> storage == NULL) - { - free(new_matrix); - return NULL; - } - - new_matrix -> rows = dim; - new_matrix -> cols = dim; - - for(ind = 0; ind < dim*dim; ind ++) - { - mpz_init_set_si (new_matrix -> storage[ind], 0); - } - - for(ind = 0; ind < dim; ind ++) - { - mpz_set_ui(new_matrix -> storage[ind*(dim+1)], 1); - } - - return new_matrix; -} - - -gmp_matrix * -create_gmp_matrix_zero(size_t rows, size_t cols) -{ - gmp_matrix * new_matrix; - size_t ind; - - new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix)); - if(new_matrix == NULL) - { - return NULL; - } - - new_matrix -> storage = (mpz_t *) calloc(rows*cols, sizeof(mpz_t)); - if(new_matrix -> storage == NULL) - { - free(new_matrix); - return NULL; - } - - new_matrix -> rows = rows; - new_matrix -> cols = cols; - - for(ind = 0; ind < rows*cols; ind ++) - { - mpz_init_set_si (new_matrix -> storage[ind], 0); - } - - return new_matrix; -} - -int -destroy_gmp_matrix(gmp_matrix * m) -{ - size_t ind, nmb_storage;; - - if(m == NULL) - { - return EXIT_FAILURE; - } - - if(m -> storage == NULL) - { - free(m); - return EXIT_FAILURE; - } - - nmb_storage = (m -> rows)*(m -> cols); - - for(ind = 0; ind < nmb_storage; ind++) - { - mpz_clear(m -> storage[ind]); - } - - free(m -> storage); - free(m); - return EXIT_SUCCESS; -} - - -size_t -gmp_matrix_rows(const gmp_matrix * m) -{ - if(m == NULL) - { - return 0; - } - - return m -> rows; -} - - -size_t -gmp_matrix_cols(const gmp_matrix * m) -{ - if(m == NULL) - { - return 0; - } - - return m -> cols; -} - -/* elem <- (matrix(row, col)) */ -int -gmp_matrix_get_elem(mpz_t elem, size_t row, size_t col, - const gmp_matrix * m) -{ -#ifdef DEBUG - - if(m == NULL) - { - return EXIT_FAILURE; - } - -#endif - - mpz_set(elem, m -> storage[(col-1)*(m -> rows)+row-1]); - return EXIT_SUCCESS; -} - -int -gmp_matrix_swap_rows(size_t row1, size_t row2, gmp_matrix * m) -{ - if(m == NULL) - { - return EXIT_FAILURE; - } - if((row1 < 1) || (row1 > m->rows) || (row2 < 1) || (row2 > m->rows)) - { - return EXIT_FAILURE; - } - - /* printf("Swapping rows %i %i\n", row1, row2); */ - gmp_blas_swap(m -> cols, - &(m -> storage[row1-1]), m -> rows, - &(m -> storage[row2-1]), m -> rows); - - return EXIT_SUCCESS; -} - -int -gmp_matrix_swap_cols(size_t col1, size_t col2, gmp_matrix * m) -{ - if(m == NULL) - { - return EXIT_FAILURE; - } - if((col1 < 1) || (col1 > m->cols) || (col2 < 1) || (col2 > m->cols)) - { - return EXIT_FAILURE; - } - - /* printf("Swapping cols %i %i\n", col1, col2); */ - gmp_blas_swap(m -> rows, - &(m -> storage[(m->rows)*(col1-1)]), 1, - &(m -> storage[(m->rows)*(col2-1)]), 1); - - return EXIT_SUCCESS; -} - -int -gmp_matrix_negate_row(size_t row, gmp_matrix * m) -{ - mpz_t minus_one; - - if(m == NULL) - { - return EXIT_FAILURE; - } - if((row < 1) || (row > m->rows)) - { - return EXIT_FAILURE; - } - - mpz_init(minus_one); - mpz_set_si(minus_one, -1); - gmp_blas_scal(m -> cols, minus_one, (&m -> storage[row-1]), m -> rows); - mpz_clear(minus_one); - return EXIT_SUCCESS; -} - -int -gmp_matrix_negate_col(size_t col, gmp_matrix * m) -{ - mpz_t minus_one; - if(m == NULL) - { - return EXIT_FAILURE; - } - if((col < 1) || (col > m->cols)) - { - return EXIT_FAILURE; - } - - mpz_init(minus_one); - mpz_set_si(minus_one, -1); - gmp_blas_scal(m -> rows, minus_one, - &(m -> storage[(m->rows)*(col-1)]), 1); - mpz_clear(minus_one); - return EXIT_SUCCESS; -} - -/* row2 <- a*row1 + row2*/ -int -gmp_matrix_add_row(mpz_t a, size_t row1, size_t row2, - gmp_matrix * m) -{ - if(m == NULL) - { - return EXIT_FAILURE; - } - if((row1 < 1) || (row1 > m->rows) || (row2 < 1) || (row2 > m->rows)) - { - return EXIT_FAILURE; - } - - gmp_blas_axpy(m->cols, a, - (const mpz_t *) &(m->storage[row1-1]), m->rows, - &(m->storage[row2-1]), m->rows); - - return EXIT_SUCCESS; -} - -int -gmp_matrix_add_col(mpz_t a, size_t col1, size_t col2, - gmp_matrix * m) -{ - if(m == NULL) - { - return EXIT_FAILURE; - } - if((col1 < 1) || (col1 > m->cols) || (col2 < 1) || (col2 > m->cols)) - { - return EXIT_FAILURE; - } - - gmp_blas_axpy(m->rows, a, - (const mpz_t *) &(m -> storage[(m->rows)*(col1-1)]), 1, - &(m -> storage[(m->rows)*(col2-1)]), 1); - - return EXIT_SUCCESS; -} - -/* row1 <- a*row1 + b*row2 - row2 <- c*row1 + d*row2 */ -int -gmp_matrix_row_rot(mpz_t a, mpz_t b, size_t row1, - mpz_t c, mpz_t d, size_t row2, - gmp_matrix * m) -{ - if(m == NULL) - { - return EXIT_FAILURE; - } - if((row1 < 1) || (row1 > m->rows) || (row2 < 1) || (row2 > m->rows)) - { - return EXIT_FAILURE; - } - - gmp_blas_rot(m->cols, - a, b, &(m->storage[row1-1]), m->rows, - c, d, &(m->storage[row2-1]), m->rows); - - return EXIT_SUCCESS; -} -int -gmp_matrix_col_rot(mpz_t a, mpz_t b, size_t col1, - mpz_t c, mpz_t d, size_t col2, - gmp_matrix * m) -{ - if(m == NULL) - { - return EXIT_FAILURE; - } - if((col1 < 1) || (col1 > m->cols) || (col2 < 1) || (col2 > m->cols)) - { - return EXIT_FAILURE; - } - -/* printf("a: %i b: %i c: %i d:%i col1: %i col2: %i\n", */ -/* mpz_get_si(a), */ -/* mpz_get_si(b), */ -/* mpz_get_si(c), */ -/* mpz_get_si(d), */ -/* col1, col2); */ - - gmp_blas_rot(m->rows, - a, b, &(m -> storage[(m->rows)*(col1-1)]), 1, - c, d, &(m -> storage[(m->rows)*(col2-1)]), 1); - - return EXIT_SUCCESS; -} - -size_t -gmp_matrix_col_inz (size_t r1, size_t r2, size_t c, - gmp_matrix * m) -{ - size_t result; - - if(m == NULL) - { - return 0; - } - if((r1 < 1) || (r1 > m->rows) || - (r2 < 1) || (r2 > m->rows) || - (r2 < r1) || (c < 1) || (c > m->cols)) - { - return 0; - } - - result = gmp_blas_inz(r2-r1+1, - (const mpz_t *) &(m->storage[(c-1)*(m->rows)+r1-1]), - 1); - - if(result > r2-r1+1) - { - return 0; - } - - return result; -} - -size_t -gmp_matrix_row_inz (size_t r, size_t c1, size_t c2, - gmp_matrix * m) -{ - size_t result; - - if(m == NULL) - { - return 0; - } - if((r < 1) || (r > m->rows) || - (c1 < 1) || (c1 > m->cols) || - (c2 < c1) || (c2 < 1) || (c2 > m->cols)) - { - return 0; - } - - result = gmp_blas_inz(c2-c1+1, - (const mpz_t *) &(m->storage[(c1-1)*(m->rows)+r-1]), - m->rows); - - if(result > c2-c1+1) - { - return 0; - } - - return result; -} - - -int -gmp_matrix_is_diagonal(const gmp_matrix * M) -{ - size_t i,j; - size_t rows, cols; - - if(M == NULL) - { - return 0; - } - - rows = M->rows; - cols = M->cols; - - for(j = 1; j <= cols; j ++) - { - for(i = 1; i <= rows; i ++) - { - if((mpz_cmp_si(M->storage[(i-1)+(j-1)*rows], 0) != 0) && - (i != j)) - { - return 0; - } - } - } - - return 1; -} - - -int -gmp_matrix_transp(gmp_matrix * M) -{ - mpz_t * new_storage; - size_t i,j; - size_t rows, cols; - - if(M == NULL) - { - return EXIT_FAILURE; - } - - rows = M->rows; - cols = M->cols; - - new_storage = (mpz_t *) calloc(rows*cols, sizeof(mpz_t)); - if(new_storage == NULL) - { - return EXIT_FAILURE; - } - - for(i = 1; i <= rows; i++) - { - for(j = 1; j <= cols; j++) - { - mpz_init_set(new_storage[(j-1)+(i-1)*cols], - M-> storage[(i-1)+(j-1)*rows]); - mpz_clear(M-> storage[(i-1)+(j-1)*rows]); - } - } - free(M->storage); - - M -> storage = new_storage; - M -> rows = cols; - M -> cols = rows; - - return EXIT_SUCCESS; -} - - -int -gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B) -{ - mpz_t * new_storage; - size_t i,j; - size_t rows_A, cols_A, rows_B, cols_B; - - if((A == NULL) || (B == NULL)) - { - return EXIT_FAILURE; - } - - rows_A = A->rows; - cols_A = A->cols; - rows_B = B->rows; - cols_B = B->cols; - - if(cols_A != rows_B) - { - return EXIT_FAILURE; - } - - /* Create new storage for the product */ - new_storage = (mpz_t *) calloc(rows_A*cols_B, sizeof(mpz_t)); - if(new_storage == NULL) - { - return EXIT_FAILURE; - } - - /* Compute the product to the storage */ - for(j = 1; j <= cols_B; j++) - { - for(i = 1; i <= rows_A; i++) - { - mpz_init (new_storage[(i-1)+(j-1)*rows_A]); - gmp_blas_dot(&(new_storage[(i-1)+(j-1)*rows_A]), - cols_A, - (const mpz_t *) &(A->storage[i-1]), rows_A, - (const mpz_t *) &(B->storage[(j-1)*rows_B]), 1); - } - } - - /* Get rid of the old storage */ - for(i = 1; i <= rows_A*cols_A; i++) - { - mpz_clear (A->storage[i-1]); - } - free(A->storage); - - /* Update A */ - A -> storage = new_storage; - A -> cols = cols_B; - - return EXIT_SUCCESS; -} - -int -gmp_matrix_left_mult(const gmp_matrix * A, gmp_matrix * B) -{ - mpz_t * new_storage; - size_t i,j; - size_t rows_A, cols_A, rows_B, cols_B; - - if((A == NULL) || (B == NULL)) - { - return EXIT_FAILURE; - } - - rows_A = A->rows; - cols_A = A->cols; - rows_B = B->rows; - cols_B = B->cols; - - if(cols_A != rows_B) - { - return EXIT_FAILURE; - } - - /* Create new storage for the product */ - new_storage = (mpz_t *) calloc(rows_A*cols_B, sizeof(mpz_t)); - if(new_storage == NULL) - { - return EXIT_FAILURE; - } - - /* Compute the product to the storage */ - for(j = 1; j <= cols_B; j++) - { - for(i = 1; i <= rows_A; i++) - { - mpz_init (new_storage[(i-1)+(j-1)*rows_A]); - gmp_blas_dot(&(new_storage[(i-1)+(j-1)*rows_A]), - cols_A, - (const mpz_t *) &(A->storage[i-1]), rows_A, - (const mpz_t *) &(B->storage[(j-1)*rows_B]), 1); - } - } - - /* Get rid of the old storage */ - for(i = 1; i <= rows_B*cols_B; i++) - { - mpz_clear (B->storage[i-1]); - } - free(B->storage); - - /* Update A */ - B -> storage = new_storage; - B -> rows = rows_A; - - return EXIT_SUCCESS; -} - - -int gmp_matrix_printf(const gmp_matrix * m) -{ - - mpz_t outz; - size_t rows, cols, i, j; - - if(m == NULL) - { - return EXIT_FAILURE; - } - - rows = m->rows; - cols = m->cols; - mpz_init(outz); - for(i = 1; i <= rows ; i++) - { - for(j = 1; j <= cols ; j++) - { - gmp_matrix_get_elem(outz, i, j, m); - mpz_out_str (stdout, 10, outz); - printf(" "); - } - printf("\n"); - } - - mpz_clear(outz); - - return EXIT_SUCCESS; -} - -int gmp_matrix_fprintf(FILE* stream, const gmp_matrix * m) -{ - mpz_t outz; - size_t rows, cols, i, j; - - if(m == NULL) - { - return EXIT_FAILURE; - } - - rows = m->rows; - cols = m->cols; - mpz_init(outz); - for(i = 1; i <= rows ; i++) - { - for(j = 1; j <= cols ; j++) - { - gmp_matrix_get_elem(outz, i, j, m); - mpz_out_str (stream, 10, outz); - printf(" "); - } - printf("\n"); - } - - mpz_clear(outz); - - return EXIT_SUCCESS; -} diff --git a/contrib/kbipack/gmp_matrix.h~ b/contrib/kbipack/gmp_matrix.h~ deleted file mode 100755 index 217b5fe7e1d1279cf617579cdc2abefaedfc9c31..0000000000000000000000000000000000000000 --- a/contrib/kbipack/gmp_matrix.h~ +++ /dev/null @@ -1,144 +0,0 @@ -/* - Header file for integer-oriented matrix, relying on the arbitrary - precision integers from the GNU gmp package. - - Copyright (C) 28.10.2003 Saku Suuriniemi TUT/CEM - - This program 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 2 of the License, or - any later version. - - This program 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 this program; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - Saku Suuriniemi, TUT/Electromagetics - P.O.Box 692, FIN-33101 Tampere, Finland - saku.suuriniemi@tut.fi - - $Id: gmp_matrix.h~,v 1.1 2009-03-30 14:10:57 matti Exp $ -*/ - -#ifndef __GMP_MATRIX_H__ -#define __GMP_MATRIX_H__ - -#include"gmp_blas.h" - -typedef struct -{ - size_t rows; - size_t cols; - mpz_t * storage; -} gmp_matrix; - -int gmp_kaka(int laalaa); - -/* Sets the values of "elems" column by column. The user is - responsible for sufficient supply. */ -gmp_matrix * -create_gmp_matrix(size_t rows, size_t cols, - const mpz_t * elems); - -gmp_matrix * -create_gmp_matrix_identity(size_t dim); - -gmp_matrix * -create_gmp_matrix_zero(size_t rows, size_t cols); - -int -destroy_gmp_matrix(gmp_matrix *); - -size_t -gmp_matrix_rows(const gmp_matrix *); - -size_t -gmp_matrix_cols(const gmp_matrix *); - -/* elem <- (matrix(row, col)) */ -int -gmp_matrix_get_elem(mpz_t elem, size_t row, size_t col, - const gmp_matrix *); - -int -gmp_matrix_swap_rows(size_t row1, size_t row2, gmp_matrix *); - -int -gmp_matrix_swap_cols(size_t col1, size_t col2, gmp_matrix *); - -int -gmp_matrix_negate_row(size_t row, gmp_matrix *); - -int -gmp_matrix_negate_col(size_t col, gmp_matrix *); - -/* row2 <- a*row1 + row2*/ -int -gmp_matrix_add_row(mpz_t a, size_t row1, size_t row2, - gmp_matrix *); -int -gmp_matrix_add_col(mpz_t a, size_t col1, size_t col2, - gmp_matrix *); - -/* row1 <- a*row1 + b*row2 - row2 <- c*row1 + d*row2 */ -int -gmp_matrix_row_rot(mpz_t a, mpz_t b, size_t row1, - mpz_t c, mpz_t d, size_t row2, - gmp_matrix *); -int -gmp_matrix_col_rot(mpz_t a, mpz_t b, size_t col1, - mpz_t c, mpz_t d, size_t col2, - gmp_matrix *); - -/* 0 for no, 1 for yes */ -int -gmp_matrix_is_diagonal(const gmp_matrix * M); - -/* Finds a nonzero in a subcolumn M(r1:r2,c). */ -/* Returns zero if no nonzeros found. */ -size_t -gmp_matrix_col_inz (size_t r1, size_t r2, size_t c, - gmp_matrix * M); - -/* Finds a nonzero in a subrow M(r,c1:c2). */ -/* Returns zero if no nonzeros found. */ -size_t -gmp_matrix_row_inz (size_t r, size_t c1, size_t c2, - gmp_matrix * M); - -int -gmp_matrix_transp(gmp_matrix * M); - -/* A <- A*B */ -int -gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B); - -/* B <- A*B */ -int -gmp_matrix_left_mult(const gmp_matrix * A, gmp_matrix * B); - -/* (TBD ?)Implement the BLAS style GEMM? Place it here at all? */ -/* (T) (T) */ -/* A <- B * C */ -/* Give 1 if the transpose of the matrix is to be used, else 0 */ -/* int */ -/* gmp_matrix_mult(gmp_matrix * A, */ -/* const gmp_matrix * B, int transpose_B, */ -/* const gmp_matrix * C, int transpose_C); */ - -/* - Mainly for diagnostics ... - -------------------------- -*/ - -int gmp_matrix_printf(const gmp_matrix *); -int gmp_matrix_fprintf(FILE*, const gmp_matrix *); - - -#endif diff --git a/contrib/kbipack/gmp_matrix.o b/contrib/kbipack/gmp_matrix.o deleted file mode 100755 index b151210c86fa837e700afce35afff9326bbec789..0000000000000000000000000000000000000000 Binary files a/contrib/kbipack/gmp_matrix.o and /dev/null differ diff --git a/contrib/kbipack/gmp_matrix_io.o b/contrib/kbipack/gmp_matrix_io.o deleted file mode 100755 index c7231192ab4007c70d68c083fb81cab8a3fddd0f..0000000000000000000000000000000000000000 Binary files a/contrib/kbipack/gmp_matrix_io.o and /dev/null differ diff --git a/contrib/kbipack/gmp_normal_form.c~ b/contrib/kbipack/gmp_normal_form.c~ deleted file mode 100755 index 5f695a53d00542bd82baea3889dcebefac08aea4..0000000000000000000000000000000000000000 --- a/contrib/kbipack/gmp_normal_form.c~ +++ /dev/null @@ -1,784 +0,0 @@ -/* - Implementation for integer computation of Hermite and Smith normal - forms of matrices of modest size. - - Implementation: Dense matrix with GMP-library's mpz_t elements to - hold huge integers. - - Algorithm: Kannan - Bachem algorithm with improvement by - Chou and Collins. Expects a large number of unit invariant - factors and takes advantage of them as they appear. - - References: - [1] Ravindran Kannan, Achim Bachem: - "Polynomial algorithms for computing the Smith and Hermite normal - forms of an integer matrix", - SIAM J. Comput., vol. 8, no. 5, pp. 499-507, 1979. - [2] Tsu-Wu J.Chou, George E. Collins: - "Algorithms for the solution of systems of linear Diophantine - equations", - SIAM J. Comput., vol. 11, no. 4, pp. 687-708, 1982. - [3] GMP homepage http://www.swox.com/gmp/ - [4] GNU gmp page http://www.gnu.org/software/gmp/ - - Copyright (C) 30.10.2003 Saku Suuriniemi TUT/CEM - - This program 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 2 of the License, or - any later version. - - This program 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 this program; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - Saku Suuriniemi, TUT/Electromagetics - P.O.Box 692, FIN-33101 Tampere, Finland - saku.suuriniemi@tut.fi - - $Id: gmp_normal_form.c~,v 1.1 2009-03-30 14:10:57 matti Exp $ -*/ - -#include<stdlib.h> -#include"gmp_blas.h" -#include"gmp_matrix.h" -#include"gmp_normal_form.h" - - -/* The unaltered matrix and identity left and right factors */ -static gmp_normal_form * -create_gmp_trivial_normal_form(gmp_matrix * A, - inverted_flag left_inverted, - inverted_flag right_inverted) -{ - size_t rows, cols; - - gmp_normal_form * new_nf; - if(A == NULL) - { - return NULL; - } - - new_nf = (gmp_normal_form *) malloc(sizeof(gmp_normal_form)); - if(new_nf == NULL) - { - destroy_gmp_matrix(A); - return NULL; - } - - rows = A -> rows; - cols = A -> cols; - - if((rows == 0) || (cols == 0)) - { - destroy_gmp_matrix(A); - return NULL; - } - - new_nf->left = create_gmp_matrix_identity(rows); - if(new_nf->left == NULL) - { - destroy_gmp_matrix(A); - free(new_nf); - return NULL; - } - - new_nf->right = create_gmp_matrix_identity(cols); - if(new_nf->right == NULL) - { - destroy_gmp_matrix(A); - destroy_gmp_matrix(new_nf->left); - free(new_nf); - return NULL; - } - - new_nf -> canonical = A; - new_nf -> left_inverted = left_inverted; - new_nf -> right_inverted = right_inverted; - - return new_nf; -} - -static int -gmp_Hermite_eliminate_step(gmp_matrix * L, gmp_matrix * U, - size_t col, inverted_flag right_inverted) -{ - size_t ind, row_limit; - mpz_t pivot, elem; - mpz_t bez1, bez2, gc_div; - mpz_t cff1, cff2; - - mpz_init(pivot); - mpz_init(elem); - mpz_init(bez1); - mpz_init(bez2); - mpz_init(cff1); - mpz_init(cff2); - mpz_init(gc_div); - - row_limit = (L->rows >= col) ? - col-1 : - L->rows; - - for(ind = 1; ind <= row_limit; ind++) - { - gmp_matrix_get_elem(elem, ind, col, L); - - /* Eliminate only if nonzero */ - if(mpz_sgn (elem) != 0) - { - gmp_matrix_get_elem(pivot, ind, ind, L); - - /* Extended Euclid's: - bez1*pivot+bez2*elem = gc_div */ - gmp_blas_rotg(gc_div, bez1, bez2, pivot, elem); - - /* Make cff1 = -elem/gc_div */ - mpz_divexact(cff1, elem, gc_div); - mpz_neg (cff1, cff1); - /* Make cff2 = pivot/gc_div */ - mpz_divexact(cff2, pivot, gc_div); - - /* Update the HNF canonical matrix */ - gmp_matrix_col_rot(bez1, bez2, ind, - cff1, cff2, col, - L); - - /* Update the unimodular factor matrix */ - if(right_inverted == INVERTED) - { - gmp_matrix_col_rot(bez1, bez2, ind, - cff1, cff2, col, - U); - } - else - { - - /* [bez1, -elem/gc_div] [pivot/gc_div, elem/gc_div] - [bez2, pivot/gc_div] [-bez2, bez1 ] = I_2 */ - mpz_neg(cff1, cff1); - mpz_neg(bez2, bez2); - gmp_matrix_row_rot(cff2, cff1, ind, - bez2, bez1, col, - U); - } - } - } - mpz_clear(pivot); - mpz_clear(elem); - mpz_clear(bez1); - mpz_clear(bez2); - mpz_clear(cff1); - mpz_clear(cff2); - mpz_clear(gc_div); - - return EXIT_SUCCESS; -} - - - -static int -gmp_Hermite_reduce_step(gmp_matrix * L, gmp_matrix * U, - size_t col, inverted_flag right_inverted) -{ - - size_t i, j; - mpz_t pivot, elem; - mpz_t cff; - - mpz_init(pivot); - mpz_init(elem); - mpz_init(cff); - - if(col > (L->rows)) - { - return EXIT_FAILURE; - } - - /* printf("Col = %i\n", col); - printf("L before\n"); - gmp_matrix_printf(L);*/ - - for(j = col-1; j >= 1; j--) - { - for(i = j+1; i <= col; i++) - { - gmp_matrix_get_elem(pivot, i, i, L); - gmp_matrix_get_elem(elem, i, j, L); - /* printf(" i %i j %i\n",i,j); */ - - if((mpz_cmpabs(elem, pivot) >= 0) || (mpz_sgn(elem) < 0)) - { - /* The objective: - 0 <= elem + k*pivot < pivot - -elem <= k*pivot < pivot - elem - -elem/pivot <= k < - elem/pivot + 1 - - Solution: - k = ceil(-elem/pivot) - */ - - /* Make cff = -elem */ - mpz_neg (cff, elem); - mpz_cdiv_q (cff, cff, pivot); - - /* printf("col %i j %i\n", i, j); - printf("elem %i k %i pivot %i\n", - mpz_get_si(elem), - mpz_get_si(cff), - mpz_get_si(pivot));*/ - - - gmp_matrix_add_col(cff, i, j, L); - - /* Update the unimodular factor matrix */ - if(right_inverted == INVERTED) - { - gmp_matrix_add_col(cff, i, j, U); - } - - /* [1 0] [ 1 0] = [1 0] - [cff 1] [-cff 1] [0 1] */ - else - { - mpz_neg (cff, cff); - gmp_matrix_add_row(cff, j, i, U); - } - - /* printf("\n"); - gmp_matrix_printf(L);*/ - - } - } - } - -/* printf("L after\n"); */ -/* gmp_matrix_printf(L); */ - - mpz_clear(pivot); - mpz_clear(elem); - mpz_clear(cff); - return EXIT_SUCCESS; -} - - -static int -gmp_normal_form_make_Hermite(gmp_normal_form * nf) -{ - size_t rows, cols; - size_t pivot_ind; - size_t schur, colind; - mpz_t pivot; - gmp_matrix * canonical; - gmp_matrix * left; - gmp_matrix * right; - - if(nf == NULL) - { - return EXIT_FAILURE; - } - - /* OK, it's safe to get to business */ - canonical = nf->canonical; - left = nf->left; - right = nf->right; - rows = canonical -> rows; - cols = canonical -> cols; - - mpz_init(pivot); - - /* "schur" denotes where the present Schur complement starts */ - schur = 1; - - while((schur <= rows) && (schur <= cols)) - { - /* Eliminate a column */ - if (schur > 1) - { - gmp_Hermite_eliminate_step(canonical, right, - schur, nf->right_inverted); - } - - /* Find a pivot */ - pivot_ind = gmp_matrix_col_inz(schur, rows, schur, canonical); - - - /* If no nonzeros was found, the column is all zero, hence - settled with. Permute it to the end and decrement cols. */ - if(pivot_ind == 0) - { - gmp_matrix_swap_cols(schur, cols, canonical); - - if(nf -> right_inverted == INVERTED) - { - gmp_matrix_swap_cols(schur, cols, right); - } - else - { - gmp_matrix_swap_rows(schur, cols, right); - } - - cols--; - - /* When the whole column was zeroed, the diagonal - elements may have got reduced. Reduce the sub- - diagonals as well*/ - - if(schur > 1) - { - gmp_Hermite_reduce_step (canonical, right, schur-1, - nf -> right_inverted); - } - } - - /* A nonzero pivot was found. Permute it to the diagonal position, - make it positive, and reduce the off-diagonals. - The schur complement now starts from the next diagonal. */ - else - { - pivot_ind = schur+pivot_ind-1; - gmp_matrix_swap_rows(schur, pivot_ind, canonical); - - if(nf->left_inverted == INVERTED) - { - gmp_matrix_swap_rows(schur, pivot_ind, left); - } - else - { - gmp_matrix_swap_cols(schur, pivot_ind, left); - } - - /* Make the pivot positive */ - gmp_matrix_get_elem(pivot, schur, schur, canonical); - - if(mpz_cmp_si(pivot, 0) < 0) - { - gmp_matrix_negate_col(schur, canonical); - - if(nf->right_inverted == INVERTED) - { - gmp_matrix_negate_col(schur, right); - } - else - { - gmp_matrix_negate_row(schur, right); - } - } - - /* Reduce off-diagonals */ - gmp_Hermite_reduce_step (canonical, right, schur, nf -> right_inverted); - - schur++; - } - } - - /* The Schur complement is now empty. There may still be uneliminated - columns left (in case of a wide matrix) */ - - colind = schur; - while (colind <= cols) - { - gmp_Hermite_eliminate_step (canonical, right, colind, nf->right_inverted); - gmp_Hermite_reduce_step (canonical, right, schur-1, nf->right_inverted); - colind++; - } - - mpz_clear(pivot); - - return EXIT_SUCCESS; -} - - - -gmp_normal_form * -create_gmp_Hermite_normal_form(gmp_matrix * A, - inverted_flag left_inverted, - inverted_flag right_inverted) -{ - gmp_normal_form * new_nf; - - if(A == NULL) - { - return NULL; - } - - new_nf = - create_gmp_trivial_normal_form(A, left_inverted, right_inverted); - - if(new_nf == NULL) - { - /* A has been destroyed already */ - return NULL; - } - - if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - return NULL; - } - - return new_nf; -} - - -gmp_normal_form * -create_gmp_Smith_normal_form(gmp_matrix * A, - inverted_flag left_inverted, - inverted_flag right_inverted) -{ - - gmp_normal_form * new_nf; - gmp_matrix * canonical; - gmp_matrix * elbow; - inverted_flag elbow_flag; - size_t rows, cols; - size_t i, j; - size_t subd_ind, row_undivisible; - size_t last_ready_row, first_ready_col, ind; - mpz_t pivot; - mpz_t remainder; - - if(A == NULL) - { - return NULL; - } - - new_nf = - create_gmp_trivial_normal_form(A, left_inverted, right_inverted); - - if(new_nf == NULL) - { - /* A has been destroyed already */ - return NULL; - } - - canonical = new_nf -> canonical; - mpz_init(pivot); - mpz_init(remainder); - rows = canonical->rows; - cols = canonical->cols; - last_ready_row = 0; - first_ready_col = (cols < rows) ? (cols+1) : (rows+1); - - while(first_ready_col > last_ready_row + 1) - { - /*******/ - /* HNF */ - /*******/ - - if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - -#ifdef DEBUG - printf("HNF\n"); - gmp_matrix_printf (new_nf -> left); - gmp_matrix_printf (new_nf -> canonical); - gmp_matrix_printf (new_nf -> right); -#endif - - /**********************/ - /* Find ready columns */ - /**********************/ - - /* If a diagonal entry is zero, so is the corresponding - column. The zero diagonals always reside in the end. - Seek until zero diagonal encountered, but stay within the matrix! */ - ind = 1; - while ( (mpz_cmp_si(canonical -> storage[(ind-1)+(ind-1)*rows], 0) - != 0) - && - (ind < first_ready_col) ) - { - ind ++; - } - first_ready_col = ind; - - /* Note: The number of ready cols is settled after the first HNF, - but the check is cheap. */ - - /**********************************************/ - /* Permute unit diagonals such that they lead */ - /**********************************************/ - - /* If the recently computed HNF has ones on the diagonal, their - corresponding rows are all zero (except the diagonal). - They are then settled, because the next LHNF kills the elements - on their columns. */ - - ind = last_ready_row+1; - - /* Stay within the nonzero cols of the matrix */ - while (ind < first_ready_col) - { - /* Unit diagonal encountered */ - if(mpz_cmp_si ( canonical->storage[(ind-1) + (ind-1)*rows], - 1) == 0) - { - /* If not in the beginning, permute to extend the leading minor - with unit diagonals */ - if(ind != last_ready_row+1) - { - gmp_matrix_swap_rows(last_ready_row+1, ind, - new_nf -> canonical); - - if(left_inverted == INVERTED) - { - gmp_matrix_swap_rows(last_ready_row+1, ind, - new_nf -> left); - } - else - { - gmp_matrix_swap_cols(last_ready_row+1, ind, - new_nf -> left); - } - - gmp_matrix_swap_cols(last_ready_row+1, ind, - new_nf -> canonical); - - if(right_inverted == INVERTED) - { - gmp_matrix_swap_cols(last_ready_row+1, ind, - new_nf -> right); - } - else - { - gmp_matrix_swap_rows(last_ready_row+1, ind, - new_nf -> right); - } - } - last_ready_row ++; - } - ind++; - } - -#ifdef DEBUG - printf("Leading units\n"); - gmp_matrix_printf (new_nf -> left); - gmp_matrix_printf (new_nf -> canonical); - gmp_matrix_printf (new_nf -> right); -#endif - - /********/ - /* LHNF */ - /********/ - - /* Extravagant strategy: Compute HNF of an explicit tranpose. */ - - /* 1. Transpose everything in the normal form */ - if(gmp_matrix_transp(canonical) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - if(gmp_matrix_transp(new_nf->left) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - if(gmp_matrix_transp(new_nf->right) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - - elbow = new_nf ->left; - new_nf ->left = new_nf ->right; - new_nf ->right = elbow; - elbow_flag = new_nf ->left_inverted; - new_nf ->left_inverted = new_nf ->right_inverted; - new_nf ->right_inverted = elbow_flag; - - /* 2. Compute HNF of the transpose */ - if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - - /* 3. Transpose everything back */ - if(gmp_matrix_transp(canonical) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - if(gmp_matrix_transp(new_nf->left) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - if(gmp_matrix_transp(new_nf->right) != EXIT_SUCCESS) - { - destroy_gmp_normal_form(new_nf); - mpz_clear(pivot); - mpz_clear(remainder); - return NULL; - } - elbow = new_nf ->left; - new_nf ->left = new_nf ->right; - new_nf ->right = elbow; - elbow_flag = new_nf ->left_inverted; - new_nf ->left_inverted = new_nf ->right_inverted; - new_nf ->right_inverted = elbow_flag; - - -#ifdef DEBUG - printf("LHNF\n"); - gmp_matrix_printf (new_nf -> left); - gmp_matrix_printf (new_nf -> canonical); - gmp_matrix_printf (new_nf -> right); -#endif - - - /*****************************************************/ - /* Check if more of the leading normal form is ready */ - /*****************************************************/ - - /* The matrix is in LHNF, i.e. it is upper triangular. - If the row trailing the top left diagonal element is - zero, the diagonal element may be an invariant factor - on its final position, and the stage may be ready. - - The stage may not still be ready: The leading diagonal element - of D may not divide the rest of the Schur complement. */ - - subd_ind = 0; - row_undivisible = 0; - - /* Explanation of loop conditions: - 1.) No relative primes found from Schur complement - 2.) Stay within the Schur complement - 3.) If a nonzero is found from the trailing row, the stage is - definitely not ready */ - while (row_undivisible == 0 && - last_ready_row + 1 < first_ready_col && - subd_ind == 0) - { - subd_ind = - gmp_matrix_row_inz(last_ready_row+1, - last_ready_row+2, cols, - canonical); - - /* printf("subd_ind %i\n", subd_ind); - printf("last_ready_row %i\n", last_ready_row); */ - - /* No nonzeros found, the stage may be ready */ - if (subd_ind == 0) - { - mpz_set (pivot, - canonical->storage[(last_ready_row)+ - (last_ready_row)*rows]); - - /* Check whether the pivot divides all elements in the Schur - complement */ - row_undivisible = 0; - for(j = last_ready_row+2; j < first_ready_col; j++) - { - for(i = last_ready_row+2; i <= j; i++) - { - mpz_tdiv_r (remainder, - canonical->storage[(i-1)+ - (j-1)*rows], - pivot); - - if(mpz_cmp_si(remainder, 0) !=0) - { - row_undivisible = i; - /* col_undivisible = j; */ - } - } - } - - /* printf("Row undivisible %i\n", row_undivisible); */ - - /* If a relative prime was found from the Schur complement, - add that row to the first row of the Schur complement */ - if(row_undivisible != 0) - { - mpz_set_si (remainder, 1); - gmp_matrix_add_row(remainder, - row_undivisible, last_ready_row+1, - canonical); - - /* [ 1 0] [1 0] = [1 0] - [-1 1] [1 1] [0 1] */ - - if(left_inverted == INVERTED) - { - gmp_matrix_add_row(remainder, - row_undivisible, last_ready_row+1, - new_nf->left); - } - else - { - mpz_neg (remainder, remainder); - gmp_matrix_add_col(remainder, - last_ready_row+1, row_undivisible, - new_nf->left); - } - } - /* The Schur complement is smaller by one again */ - else - { - last_ready_row++; - } - } - } - } /* The main loop ends here */ - - mpz_clear(pivot); - mpz_clear(remainder); - return new_nf; -} - -int destroy_gmp_normal_form(gmp_normal_form * nf) -{ - int status; - - if(nf == NULL) - { - return EXIT_FAILURE; - } - - status = EXIT_SUCCESS; - if(destroy_gmp_matrix(nf -> left) == EXIT_FAILURE) - { - status = EXIT_FAILURE; - } - - if(destroy_gmp_matrix(nf -> canonical) == EXIT_FAILURE) - { - status = EXIT_FAILURE; - } - - if(destroy_gmp_matrix(nf -> right) == EXIT_FAILURE) - { - status = EXIT_FAILURE; - } - free(nf); - - return status; -} - - diff --git a/contrib/kbipack/gmp_normal_form.o b/contrib/kbipack/gmp_normal_form.o deleted file mode 100755 index f165a904233fc099325f1bb789f948c5207669de..0000000000000000000000000000000000000000 Binary files a/contrib/kbipack/gmp_normal_form.o and /dev/null differ