From ac0b79e376e42f7e88175a3c9164356035164fe7 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Sat, 13 Feb 2010 10:09:35 +0000
Subject: [PATCH] Patch by Richard.
---
Geo/ChainComplex.h | 6 +-
contrib/kbipack/CMakeLists.txt | 8 +-
contrib/kbipack/compute_normal_form.cpp | 355 +++++++++++
contrib/kbipack/gmp_blas.cpp | 227 +++++++
contrib/kbipack/gmp_blas.h | 14 +-
contrib/kbipack/gmp_matrix.cpp | 793 ++++++++++++++++++++++++
contrib/kbipack/gmp_matrix.h | 2 +-
contrib/kbipack/gmp_matrix_io.cpp | 135 ++++
contrib/kbipack/gmp_matrix_io.h | 2 +-
contrib/kbipack/gmp_normal_form.cpp | 783 +++++++++++++++++++++++
contrib/kbipack/mpz.cpp | 2 +-
11 files changed, 2308 insertions(+), 19 deletions(-)
create mode 100644 contrib/kbipack/compute_normal_form.cpp
create mode 100644 contrib/kbipack/gmp_blas.cpp
create mode 100644 contrib/kbipack/gmp_matrix.cpp
create mode 100644 contrib/kbipack/gmp_matrix_io.cpp
create mode 100644 contrib/kbipack/gmp_normal_form.cpp
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index eeaec9b163..6ecb04e7b3 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -27,15 +27,11 @@
#if defined(HAVE_GMP)
-#include "gmp.h"
-extern "C" {
+ #include "gmp.h"
#include "gmp_normal_form.h"
-}
#else
-extern "C" {
#include "mpz.h"
#include "gmp_normal_form.h"
-}
#endif
// A class representing a chain complex of a cell complex.
diff --git a/contrib/kbipack/CMakeLists.txt b/contrib/kbipack/CMakeLists.txt
index ae9cf789de..253ced8236 100644
--- a/contrib/kbipack/CMakeLists.txt
+++ b/contrib/kbipack/CMakeLists.txt
@@ -4,10 +4,10 @@
# bugs and problems to <gmsh@geuz.org>.
set(SRC
- gmp_normal_form.c
- gmp_matrix_io.c
- gmp_matrix.c
- gmp_blas.c
+ gmp_normal_form.cpp
+ gmp_matrix_io.cpp
+ gmp_matrix.cpp
+ gmp_blas.cpp
mpz.cpp
)
diff --git a/contrib/kbipack/compute_normal_form.cpp b/contrib/kbipack/compute_normal_form.cpp
new file mode 100644
index 0000000000..f86e043cc2
--- /dev/null
+++ b/contrib/kbipack/compute_normal_form.cpp
@@ -0,0 +1,355 @@
+/*
+ COMPUTE_NORMAL_FORM
+
+ Hermite and Smith normal forms of integer matrices.
+
+ 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) 3.11.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: compute_normal_form.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+
+#include<stdlib.h>
+#include<stdio.h>
+#include<string.h>
+#include"gmp_matrix_io.h"
+#include"gmp_normal_form.h"
+
+typedef enum {HERMITE, SMITH} nf_flag;
+typedef enum {NOT_WANTED, WANTED} w_flag;
+
+static int nf_parse_options(int arg_counter, char ** args,
+ char ** p_infile_name, char ** p_outfile_name,
+ nf_flag * p_normal_form_requested,
+ inverted_flag * p_left_inverted,
+ inverted_flag * p_right_inverted,
+ w_flag * p_help_wanted,
+ w_flag * p_verbose_wanted)
+{
+ size_t ind1, ind2;
+ size_t outfilename_length;
+
+ /* If no infile is specified, the error is indisputable.
+ => Print help and exit.
+ Infile should be the last argument. */
+
+ if(strncmp(args[arg_counter-1], "-", 1) == 0)
+ {
+ *(p_help_wanted) = WANTED;
+ return EXIT_SUCCESS;
+ }
+
+ /* Should I find arguments? */
+ if(arg_counter > 2)
+ {
+ for(ind1 = 1; ind1 < arg_counter-1; ind1++)
+ {
+ /* No multiple infile names ! */
+ if(strncmp(args[ind1], "-", 1) != 0)
+ {
+ *(p_help_wanted) = WANTED;
+ return EXIT_SUCCESS;
+ }
+
+ /* Parse the option letters after '-' one by one */
+ for(ind2 = 1; ind2 < strlen(args[ind1]); ind2++)
+ {
+ switch(args[ind1][ind2])
+ {
+ case 'h':
+ *(p_help_wanted) = WANTED;
+ return EXIT_SUCCESS;
+ break;
+ case 'H':
+ *(p_normal_form_requested) = HERMITE;
+ break;
+ case 'S':
+ *(p_normal_form_requested) = SMITH;
+ break;
+ case 'l':
+ *(p_left_inverted) = INVERTED;
+ break;
+ case 'r':
+ *(p_right_inverted) = INVERTED;
+ break;
+ case 'v':
+ *(p_verbose_wanted) = WANTED;
+ break;
+ }
+ }
+ }
+ }
+
+ /* Get the infile name */
+ * p_infile_name = (char *) calloc(strlen(args[arg_counter-1]) + 1,
+ sizeof(char));
+ if(* p_infile_name == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ strcpy(* p_infile_name, args[arg_counter-1]);
+
+ /* Create the outfile name */
+ outfilename_length = strcspn( * p_infile_name, ".");
+ * p_outfile_name = (char *) calloc(outfilename_length + 1,
+ sizeof(char));
+ if( * p_outfile_name == NULL)
+ {
+ free( * p_infile_name);
+ return EXIT_FAILURE;
+ }
+
+ strncpy(* p_outfile_name, * p_infile_name, outfilename_length);
+ (* p_outfile_name)[outfilename_length] = '\0';
+
+ return EXIT_SUCCESS;
+}
+
+static int nf_print_help(void)
+{
+ fprintf(stdout, "Usage: compute_normal_form [-hHSlrv] infile\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Description: Compute the Hermite/Smith normal form of an integer matrix.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Input: \"infile\" is an ascii file, which contains an integer matrix\n");
+ fprintf(stdout, " in the sparse coordinate format:\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, " # All comments before the data\n");
+ fprintf(stdout, " rows columns nonzeros\n");
+ fprintf(stdout, " row1 col1 value1\n");
+ fprintf(stdout, " row2 col2 value2\n");
+ fprintf(stdout, " ...\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, " No multiple infiles, please.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Options: -h This help\n");
+ fprintf(stdout, " -H Compute the Hermite normal form (default)\n");
+ fprintf(stdout, " -S Compute the Smith normal form\n");
+ fprintf(stdout, " -l Left factor inverted (A = P^(T)LU instead of PA = LU for Hermite, \n");
+ fprintf(stdout, " U^(-1)A = SV instead of A = USV for Smith)\n");
+ fprintf(stdout, " -r Right factor inverted (PAU^(-1) = L instead of PA = LU for Hermite, \n");
+ fprintf(stdout, " AV^(-1) = US instead of A = USV for Smith)\n");
+ fprintf(stdout, " -v Verbose. Print the matrix and the decomposition to sdtout, \n");
+ fprintf(stdout, " so that the user may check that the problem is read and written\n");
+ fprintf(stdout, " correctly.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Output: For inpur matrix in e.g. \"infile.coord\", writes ascii files\n");
+ fprintf(stdout, " \"infile.left\", \"infile.can\", and \"infile.right\", which \n");
+ fprintf(stdout, " contain the left factor, canonical form, and the right factor, \n");
+ fprintf(stdout, " respectively. The matrices are written in the sparse coordinate\n");
+ fprintf(stdout, " format.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Note: The return values may be big enough to not fit the standard integer types.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Copyright (C) 3.11.2003 Saku Suuriniemi, Tampere University of Technology/Electromagnetics\n");
+ fprintf(stdout, " P.O.Box 692, FIN-33101 Tampere, Finland. saku.suuriniemi@tut.fi\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Licenced under GNU General Public Licence (GPL) version 2 or later, ABSOLUTELY NO WARRANTY.\n");
+
+ return EXIT_SUCCESS;
+}
+
+
+static char* append_suffix(const char * file_base_name, const char * suffix)
+{
+ char * filename;
+
+ if((file_base_name == NULL) || (suffix == NULL))
+ {
+ return NULL;
+ }
+
+ filename =
+ (char*) calloc(strlen(file_base_name)+strlen(suffix)+2,
+ sizeof(char));
+
+ if(filename == NULL)
+ {
+ return NULL;
+ }
+
+ strcpy(filename, file_base_name);
+ strcat(filename, ".");
+ strcat(filename, suffix);
+
+ return filename;
+}
+
+
+
+int main(int argc, char * argv[])
+{
+ char * infile_name;
+ char * file_base_name;
+ char * outfile_name;
+ gmp_matrix * A;
+ gmp_normal_form * normal_form;
+ w_flag help_wanted;
+ w_flag verbose_wanted;
+ nf_flag normal_form_requested;
+ inverted_flag left_inverted;
+ inverted_flag right_inverted;
+
+ if(argc == 1)
+ {
+ fprintf(stdout, "Please try \"compute_normal_form -h\" for help.\n");
+ return EXIT_SUCCESS;
+ }
+
+ /* Default settings before option parsing */
+ normal_form_requested = HERMITE;
+ left_inverted = NOT_INVERTED;
+ right_inverted = NOT_INVERTED;
+ help_wanted = NOT_WANTED;
+ verbose_wanted = NOT_WANTED;
+
+ if( nf_parse_options(argc, argv,
+ &infile_name, &file_base_name,
+ &normal_form_requested,
+ &left_inverted, &right_inverted,
+ &help_wanted,
+ &verbose_wanted)
+ == EXIT_FAILURE )
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* In this case, no dynamical memory has been allocated by
+ nf_parse_options ! */
+ if(help_wanted == WANTED)
+ {
+ return nf_print_help();
+ }
+
+ if(verbose_wanted == WANTED)
+ {
+ fprintf(stdout, "Computing ");
+ if(normal_form_requested == HERMITE)
+ {
+ fprintf(stdout, "Hermite");
+ }
+ else
+ {
+ fprintf(stdout, "Smith");
+ }
+ fprintf(stdout, " normal form\n");
+ if(left_inverted == INVERTED)
+ {
+ fprintf(stdout, "The left factor will be inverted.\n");
+ }
+ if(right_inverted == INVERTED)
+ {
+ fprintf(stdout, "The right factor will be inverted.\n");
+ }
+ fprintf(stdout, "Matrix read from file %s\n", infile_name);
+ }
+
+ /*************************************************************/
+ /* We have all we need to get to the computational business: */
+ /*************************************************************/
+
+ /* Read the matrix */
+ A = gmp_matrix_read_coord(infile_name);
+ free(infile_name);
+
+ if(verbose_wanted == WANTED)
+ {
+ fprintf(stdout, "Input matrix:\n");
+ gmp_matrix_fprintf(stdout, A);
+ }
+
+ /* Compute the requested normal form */
+ switch(normal_form_requested)
+ {
+ case SMITH:
+ normal_form = create_gmp_Smith_normal_form(A, left_inverted, right_inverted);
+ break;
+
+ default:
+ normal_form = create_gmp_Hermite_normal_form(A, left_inverted, right_inverted);
+ break;
+ }
+
+
+ if(verbose_wanted == WANTED)
+ {
+ fprintf(stdout, "Left factor:\n");
+ gmp_matrix_fprintf(stdout, normal_form->left);
+ fprintf(stdout, "Canonical form:\n");
+ gmp_matrix_fprintf(stdout, normal_form->canonical);
+ fprintf(stdout, "Right factor:\n");
+ gmp_matrix_fprintf(stdout, normal_form->right);
+ }
+
+ /* Save the result to the output files */
+ outfile_name = append_suffix(file_base_name, "left");
+ if( gmp_matrix_write_coord(outfile_name, normal_form -> left)
+ == EXIT_SUCCESS)
+ {
+ if(verbose_wanted == WANTED)
+ fprintf(stdout, "Left factor successfully writen to file %s\n", outfile_name);
+ }
+
+ free(outfile_name);
+
+ outfile_name = append_suffix(file_base_name, "can");
+ if( gmp_matrix_write_coord(outfile_name, normal_form -> canonical)
+ == EXIT_SUCCESS)
+ {
+ if(verbose_wanted == WANTED)
+ fprintf(stdout, "Canonical form successfully writen to file %s\n", outfile_name);
+ }
+ free(outfile_name);
+
+ outfile_name = append_suffix(file_base_name, "right");
+ if( gmp_matrix_write_coord(outfile_name, normal_form -> right)
+ == EXIT_SUCCESS)
+ {
+ if(verbose_wanted == WANTED)
+ fprintf(stdout, "Right factor successfully writen to file %s\n", outfile_name);
+ }
+ free(outfile_name);
+
+ free(file_base_name);
+ destroy_gmp_normal_form(normal_form);
+
+ return EXIT_SUCCESS;
+}
diff --git a/contrib/kbipack/gmp_blas.cpp b/contrib/kbipack/gmp_blas.cpp
new file mode 100644
index 0000000000..4d0dd74eab
--- /dev/null
+++ b/contrib/kbipack/gmp_blas.cpp
@@ -0,0 +1,227 @@
+/*
+ Source file for integer-oriented elementary versions of some
+ subroutines of BLAS. There routines are chosen to facilitate the
+ computation of Hermite and Smith normal forms of matrices of
+ modest size.
+
+ 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_blas.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+/* #include<stdio.h> */
+#include"gmp_blas.h"
+
+void gmp_blas_swap(size_t n, mpz_t* x, size_t incx, mpz_t* y, size_t incy)
+{
+ mpz_t elbow;
+ size_t ind;
+
+ mpz_init(elbow);
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_set(elbow , x[ind*incx]);
+ mpz_set(x[ind*incx], y[ind*incy]);
+ mpz_set(y[ind*incy], elbow);
+ }
+ mpz_clear(elbow);
+
+ return;
+}
+
+void gmp_blas_scal(size_t n, mpz_t a, mpz_t* x, size_t incx)
+{
+ size_t ind;
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_mul (x[ind*incx], x[ind*incx], a);
+ }
+
+ return;
+}
+
+void gmp_blas_copy(size_t n, mpz_t* x, size_t incx, mpz_t* y, size_t incy)
+{
+ size_t ind;
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_set(y[ind*incy], x[ind*incx]);
+ }
+
+ return;
+}
+
+void gmp_blas_axpy(size_t n, mpz_t a, mpz_t* x, size_t incx,
+ mpz_t* y, size_t incy)
+{
+ size_t ind;
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_addmul (y[ind*incy], a, x[ind*incx]);
+ }
+
+ return;
+}
+
+void
+gmp_blas_dot(mpz_t * d, size_t n,
+ mpz_t* x, size_t incx,
+ mpz_t* y, size_t incy)
+{
+ size_t ind;
+
+ mpz_set_si(*d,0);
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_addmul (*d, x[ind*incx], y[ind*incy]);
+ }
+
+ return;
+}
+
+/* x <- ax + by
+ y <- cx + dy */
+void gmp_blas_rot(size_t n,
+ mpz_t a, mpz_t b, mpz_t* x, size_t incx,
+ mpz_t c, mpz_t d, mpz_t* y, size_t incy)
+{
+ mpz_t ax;
+ mpz_t by;
+ mpz_t cx;
+ mpz_t dy;
+
+ size_t ind;
+
+ mpz_init(ax);
+ mpz_init(by);
+ mpz_init(cx);
+ mpz_init(dy);
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_mul (ax, a, x[ind*incx]);
+ mpz_mul (by, b, y[ind*incy]);
+ mpz_mul (cx, c, x[ind*incx]);
+ mpz_mul (dy, d, y[ind*incy]);
+ mpz_add (x[ind*incx], ax, by);
+ mpz_add (y[ind*incy], cx, dy);
+ }
+
+ mpz_clear(ax);
+ mpz_clear(by);
+ mpz_clear(cx);
+ mpz_clear(dy);
+
+ return;
+}
+
+
+/* Returns k such that x[(k-1)*incx] != 0 holds for the first time.
+ If none found, returns n+1. */
+size_t gmp_blas_inz (size_t n, mpz_t* x, size_t incx)
+{
+ size_t ind;
+
+ for(ind = 0; ind < n; ind++)
+ {
+ if(mpz_sgn (x[ind*incx]) != 0)
+ {
+ return ind+1;
+ }
+ }
+
+ /* No nonzeros found */
+ return n+1;
+}
+
+size_t gmp_blas_iamax(size_t n, mpz_t* x, size_t incx)
+{
+ size_t ind;
+ size_t ind_so_far;
+ mpz_t max_so_far;
+
+ mpz_init(max_so_far);
+ mpz_set_si(max_so_far, 0);
+ ind_so_far = 0;
+
+ for(ind = 0; ind < n; ind++)
+ {
+ if(mpz_cmpabs (x[ind*incx], max_so_far) > 0)
+ {
+ mpz_set(max_so_far, x[ind*incx]);
+ ind_so_far = ind;
+ }
+ }
+
+ /* No nonzeros found */
+ if(mpz_sgn (max_so_far) == 0)
+ {
+ mpz_clear(max_so_far);
+ return n + 1;
+ }
+
+ /* Nonzero maximal by modulus element found */
+ mpz_clear(max_so_far);
+ return ind_so_far + 1;
+}
+
+
+size_t gmp_blas_iamin(size_t n, mpz_t* x, size_t incx)
+{
+ size_t ind;
+ size_t ind_so_far;
+ mpz_t min_so_far;
+
+ ind_so_far = gmp_blas_inz (n, x, incx);
+
+ /* No nonzeros found? */
+ if(ind_so_far == n+1)
+ {
+ return n+1;
+ }
+
+ /* OK, There is at leat one nonzero element */
+ mpz_init(min_so_far);
+ mpz_set(min_so_far, x[(ind_so_far-1)*incx]);
+
+ /* Scan through te rest of the elements to see if smaller
+ elements by modulus are found */
+ for(ind = ind_so_far-1; ind < n; ind++)
+ {
+ if(mpz_sgn (x[ind*incx]) != 0)
+ {
+ if(mpz_cmpabs (x[ind*incx], min_so_far) < 0)
+ {
+ mpz_set(min_so_far, x[ind*incx]);
+ ind_so_far = ind + 1;
+ }
+ }
+ }
+
+ /* Nonzero minimal by modulus element found */
+ mpz_clear(min_so_far);
+ return ind_so_far;
+}
diff --git a/contrib/kbipack/gmp_blas.h b/contrib/kbipack/gmp_blas.h
index b2800e5c2f..8a70cd4c6e 100644
--- a/contrib/kbipack/gmp_blas.h
+++ b/contrib/kbipack/gmp_blas.h
@@ -47,19 +47,19 @@ gmp_blas_scal(size_t n, mpz_t a, mpz_t * x, size_t incx);
/* y <- x */
void
-gmp_blas_copy(size_t n, const mpz_t * x, size_t incx, mpz_t * y, size_t incy);
+gmp_blas_copy(size_t n, mpz_t * x, size_t incx, mpz_t * y, size_t incy);
/* y <- ax + y */
void
gmp_blas_axpy(size_t n, mpz_t a,
- const mpz_t * x, size_t incx,
+ mpz_t * x, size_t incx,
mpz_t * y, size_t incy);
/* d <- x^T y The integer * d has to be initialized by e.g. mpz_init() ! */
void
gmp_blas_dot(mpz_t * d, size_t n,
- const mpz_t * x, size_t incx,
- const mpz_t * y, size_t incy);
+ mpz_t * x, size_t incx,
+ mpz_t * y, size_t incy);
/* Givens rotations are obviously impossible. However, the extended Euclid
@@ -92,12 +92,12 @@ void gmp_blas_rot(size_t n,
/* Returns k such that x[(k-1)*incx] != 0 holds for the first time.
In FORTRAN, this is x((k-1)*incx + 1) .NE. 0.
If none found, returns n+1. */
-size_t gmp_blas_inz (size_t n, const mpz_t * x, size_t incx);
+size_t gmp_blas_inz (size_t n, mpz_t * x, size_t incx);
/* Similarly, x[(k-1)*incx] is maximal by modulus */
-size_t gmp_blas_iamax(size_t n, const mpz_t * x, size_t incx);
+size_t gmp_blas_iamax(size_t n, mpz_t * x, size_t incx);
/* Similarly, x[(k-1)*incx] is nonzero and minimal by modulus */
-size_t gmp_blas_iamin(size_t n, const mpz_t * x, size_t incx);
+size_t gmp_blas_iamin(size_t n, mpz_t * x, size_t incx);
#endif
diff --git a/contrib/kbipack/gmp_matrix.cpp b/contrib/kbipack/gmp_matrix.cpp
new file mode 100644
index 0000000000..aae01ec1e7
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix.cpp
@@ -0,0 +1,793 @@
+/*
+ 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.5 2009-05-05 11:35:01 matti Exp $
+*/
+
+
+#include<stdlib.h>
+#include"gmp_matrix.h"
+
+gmp_matrix *
+create_gmp_matrix(size_t r, size_t c,
+ 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_int(size_t r, size_t c,
+ const long int * 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_si (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;
+}
+
+gmp_matrix *
+copy_gmp_matrix(const gmp_matrix * matrix,
+ const size_t start_row, const size_t start_col,
+ const size_t end_row, const size_t end_col)
+{
+ gmp_matrix * new_matrix;
+ size_t ind;
+ size_t r;
+ size_t c;
+ size_t old_rows;
+ size_t old_cols;
+ size_t i;
+ size_t j;
+
+ new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix));
+ if(new_matrix == NULL)
+ {
+ return NULL;
+ }
+
+ r = end_row-start_row+1;
+ c = end_col-start_col+1;
+ if(r < 1 || c < 1) 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;
+
+ old_rows = matrix -> rows;
+ old_cols = matrix -> cols;
+
+ ind = 0;
+ for(j = 1; j <= old_cols; j++){
+ if(j >= start_col && j <= end_col){
+ for(i = 1; i <= old_rows; i++){
+ if(i >= start_row && i <= end_row){
+ mpz_init (new_matrix -> storage[ind]);
+ mpz_set (new_matrix -> storage[ind], matrix -> storage[(j-1)*old_rows+(i-1)]);
+ ind++;
+ }
+ }
+ }
+ }
+
+ 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;
+}
+
+/* (matrix(row, col)) <- elem */
+int
+gmp_matrix_set_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(m -> storage[(col-1)*(m -> rows)+row-1], elem);
+ 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,
+ (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,
+ (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,
+ (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,
+ (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;
+ }
+
+ if(rows == 1){
+ for(i = 1; i <= cols; i++)
+ {
+ mpz_init_set(new_storage[i-1],
+ M-> storage[i-1]);
+ mpz_clear(M-> storage[i-1]);
+ }
+ }
+ else if(cols == 1){
+ for(i = 1; i <= rows; i++)
+ {
+ mpz_init_set(new_storage[i-1],
+ M-> storage[i-1]);
+ mpz_clear(M-> storage[i-1]);
+ }
+ }
+ else{
+ 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,
+ (mpz_t *) &(A->storage[i-1]), rows_A,
+ (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,
+ (mpz_t *) &(A->storage[i-1]), rows_A,
+ (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
index ea0c9af66b..931f382803 100644
--- a/contrib/kbipack/gmp_matrix.h
+++ b/contrib/kbipack/gmp_matrix.h
@@ -41,7 +41,7 @@ typedef struct
responsible for sufficient supply. */
gmp_matrix *
create_gmp_matrix(size_t rows, size_t cols,
- const mpz_t * elems);
+ mpz_t * elems);
gmp_matrix *
create_gmp_matrix_int(size_t rows, size_t cols,
const long int * elems);
diff --git a/contrib/kbipack/gmp_matrix_io.cpp b/contrib/kbipack/gmp_matrix_io.cpp
new file mode 100644
index 0000000000..14b0fe3650
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix_io.cpp
@@ -0,0 +1,135 @@
+/*
+ IO-routines for gmp integer matrices in coordinate format.
+
+ Copyright (C) 19.11.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_io.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+ */
+
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+#include "gmp_matrix_io.h"
+
+gmp_matrix * gmp_matrix_read_coord(char* filename)
+{
+ FILE * p_file;
+ size_t read_values;
+ size_t row, col, nrows, ncols, dummy;
+ int val;
+ char buffer[1000];
+ gmp_matrix * new_matrix;
+
+ p_file = fopen(filename, "r");
+ if(p_file == NULL)
+ return NULL;
+
+ /* Matlab and Octave typically include comments on ascii files. */
+ /* They start with # */
+ fgets(buffer, 999, p_file);
+ while(strncmp(buffer, "#",1) == 0)
+ {
+ fgets(buffer, 999, p_file);
+ }
+
+ /* First read the size of the matrix */
+ read_values = sscanf(buffer, "%u %u %u", &nrows, &ncols, &dummy);
+
+ /* Create the matrix */
+ new_matrix = create_gmp_matrix_zero(nrows, ncols);
+ if(new_matrix == NULL)
+ {
+ fclose(p_file);
+ return NULL;
+ }
+
+ /* Read the values to the matrix */
+ while(read_values != EOF)
+ {
+ read_values = fscanf(p_file, "%u %u %i\n", &row, &col, &val);
+ if( (row <= nrows) && (row > 0) && (col <= ncols) && (col > 0) )
+ {
+ mpz_set_si ( new_matrix->storage[(row-1)+(col-1)*nrows], val );
+ }
+ }
+ fclose(p_file);
+
+ return new_matrix;
+}
+
+int gmp_matrix_write_coord(char* filename, gmp_matrix * M)
+{
+
+ FILE * p_file;
+ size_t written_values;
+ size_t row, col, nrows, ncols, nnz;
+ mpz_t outz;
+
+ if(M == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ nrows = M->rows;
+ ncols = M->cols;
+
+ p_file = fopen(filename, "w");
+ if(p_file == NULL)
+ return 0;
+
+ nnz = 0;
+ /* Compute the number of nonzeros for coordinate format */
+ mpz_init(outz);
+ for(col = 1; col <= ncols ; col++)
+ {
+ for(row = 1; row <= nrows ; row++)
+ {
+ gmp_matrix_get_elem(outz, row, col, M);
+ if( mpz_cmp_si (outz, 0) != 0)
+ {
+ nnz++;
+ }
+ }
+ }
+
+ /* First write the size of the matrix */
+ written_values = fprintf(p_file, "%u %u %u\n", nrows, ncols, nnz);
+
+ /* Write the values to the matrix */
+ for(col = 1; col <= ncols ; col++)
+ {
+ for(row = 1; row <= nrows ; row++)
+ {
+ gmp_matrix_get_elem(outz, row, col, M);
+ if( mpz_cmp_si (outz, 0) != 0)
+ {
+ fprintf(p_file, "%u %u ", row, col);
+ mpz_out_str (p_file, 10, outz);
+ fprintf(p_file, "\n");
+ }
+ }
+ }
+ mpz_clear(outz);
+
+ fclose(p_file);
+
+ return EXIT_SUCCESS;
+}
diff --git a/contrib/kbipack/gmp_matrix_io.h b/contrib/kbipack/gmp_matrix_io.h
index 74df3f2ad1..65c23b9bff 100644
--- a/contrib/kbipack/gmp_matrix_io.h
+++ b/contrib/kbipack/gmp_matrix_io.h
@@ -31,7 +31,7 @@
#include "gmp_matrix.h"
gmp_matrix * gmp_matrix_read_coord(char* filename);
-int gmp_matrix_write_coord(char* filename, const gmp_matrix * M);
+int gmp_matrix_write_coord(char* filename, gmp_matrix * M);
#endif
diff --git a/contrib/kbipack/gmp_normal_form.cpp b/contrib/kbipack/gmp_normal_form.cpp
new file mode 100644
index 0000000000..42785337e7
--- /dev/null
+++ b/contrib/kbipack/gmp_normal_form.cpp
@@ -0,0 +1,783 @@
+/*
+ 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 ( (ind < first_ready_col) &&
+ (mpz_cmp_si(canonical -> storage[(ind-1)+(ind-1)*rows], 0) != 0)
+ )
+ {
+ 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/mpz.cpp b/contrib/kbipack/mpz.cpp
index 3ec1b2a4df..e4658189ec 100755
--- a/contrib/kbipack/mpz.cpp
+++ b/contrib/kbipack/mpz.cpp
@@ -19,7 +19,7 @@
#if ! defined(HAVE_GMP)
-#include "GmshMessage.cpp"
+#include "GmshMessage.h"
#include "limits.h"
void overflow()
--
GitLab