PARI hangs with initialize / freeall / initialize

[Ramón Casero Cañas <8-T@gmx.net>]

For my program the inversion takes very little, and it needs memory, so at
the beginning of the function there's a pari stack initialization, and before
returning, the stack is freed. But two consecutive calls to the function hang
the Linux terminal where the program runs. A CTRL-C says

  ***   user interrupt.

but it doesn't release the terminal, and the CPU load remain 100%.  A kill
command must be executed from another terminal in order to stop the program. 

The test program is compiled and linked with g++ 2.95.4:


---------------------------------------------------------------------

/* File m_modulo2.cpp */

#include <iostream>                      // for std::cout
#include <exception>
#include <typeinfo>
// #include <fstream>                       // for ifstream
typedef enum{FALSE, TRUE} mybool;
#define exor(a, b) ((mybool)(a != b))
extern "C" {
// #include "include/gamatrix.h"
/* mybool inverse_mat_mod2(mybool **_a, unsigned _size); */
mybool inverse_mat_mod2(mybool **_a, unsigned _size, mybool _CHECK);
/* print_c_boolmatrix: print a C boolean matrix */
void print_c_boolmatrix(mybool **_a, unsigned _m, unsigned _n);
}
/* main */
/*
   argv[] (IN): file names array
   ---
   return n:   n files were succesfully read, n > 0
   return 0:   none of the files could be read
   return -1:  unexpected error
   ---
   read steiner problem from files
*/
int main(int argc, char *argv[]) {

  try {

    unsigned m, n;

    std::cout << "Hola mundo" << std::endl;
//     boost::array<BoolMatrix::index, 2> shape = {{ 5, 8 }};
//     BoolMatrix *A;
//     boost::array<BoolMatrix::index, 2> dims = {{20, 15}};
//     A.reshape(dims);
//     print_multi_array(A);

    const unsigned a_size = 32;
    const bool a_const[a_size][a_size] =
      {
        {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
         1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1},
        {1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
         1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0},
        {0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1,
         1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
        {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0,
         0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
        {0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
         1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0},
        {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
         1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0},
        {0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0,
         1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0},
        {0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
         0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1},
        {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
         0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1},
        {0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
         0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
        {0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0,
         0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1},
        {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
         0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0},
        {0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
         0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0},
        {0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0,
         0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0},
        {0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,
         1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0},
        {0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
         1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1},
        {0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1,
         1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1},
        {0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1,
         0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0},
        {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0,
         1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0},
        {0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
         0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1},
        {0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
         0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0},
        {0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1,
         0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0},
        {1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
         0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0},
        {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0,
         1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0},
        {1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1,
         0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1},
        {0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0,
         0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0},
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0,
         0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0},
        {0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
         0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0},
        {0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
         0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1},
        {0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
         0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
        {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1,
         0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0},
        {1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1,
         0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0}
      };

    mybool **a;
    a = new (mybool *)[a_size];
    for (m = 0; m < a_size; ++m) {
      a[m] = new mybool[a_size];
      for (n = 0; n < a_size; ++n) {
        a[m][n] = (mybool)a_const[m][n];
      }
    }
//     print_c_boolmatrix(a, a_size, a_size);

    if (inverse_mat_mod2(a, a_size, TRUE)) {
      std::cout << "right inversion" << std::endl;
//       print_c_boolmatrix(a, a_size, a_size);
    }
    else {
      std::cout << "wrong inversion" << std::endl;
      return 0;
    }

    if (inverse_mat_mod2(a, a_size, TRUE)) {
      std::cout << "right inversion" << std::endl;
//       print_c_boolmatrix(a, a_size, a_size);
    }
    else {
      std::cout << "wrong inversion" << std::endl;
      return 0;
    }


    return 1;
  }
  catch (char *s) {
    std::cout << "exception! " << s << std::endl;
    return 0;
  }
  catch (std::exception& e) {
    std::cout << "exception! " << e.what() << std::endl;
  }
  catch (...) {
    std::cout << "unexpected exception!" << std::endl;
    return -1;
  }

}

---------------------------------------------------------------------------


The functions file is compiled with gcc 2.95.4


---------------------------------------------------------------------

/* File gamatrix.cpp */
/* Thanks to Karim Belabas and Bill Allombert for his help with this
   code */

// #include <iostream>                      // for std::cout
#include "include/gamatrix.h"
#include <pari/pari.h>
// #include "/usr/src/pari-2.2.4.alpha/src/basemath/base3.c"
// GEN u_FpM_inv(GEN a, ulong p);

void freeall(void); // thanks, Karim
void pari_sig_init(void (*f)(int));  // thanks, Karim
// extern "C" {
//   void pari_sig_init(void (*f)(int));
// }

/* inverse_mat_mod2: overwrites a matrix modulo 2 with its inverse;
 * input is a mybool matrix, but internally, PARI library types are
 * used */
/* return true if it has inverse and false if not or cannot be
   computed */
/* mybool inverse_mat_mod2(mybool **_a, unsigned _size) { */
/*   return inverse_mat_mod2(_a, _size, true); */
/* } */

mybool inverse_mat_mod2(mybool **_a, unsigned _size, mybool _CHECK) {
  GEN a;           // temp matrix by columns
  unsigned i, j;   // counter
  unsigned sizepari = _size + 1;
  GEN z;           // auxiliary element of the matrix
  GEN v;           // auxiliary element of the matrix
  long paristacksize = 500000;
  long total_size = sizepari * sizepari + _size * _size * 9; // thanks,
Karim
  GEN determinant;
  GEN ainv;

  printf("inverse_mat_mod2: getting into pari function\n");
//   std::cout << "getting into pari function" << std::cout;

  paristacksize = sizeof(long) * total_size * 10; // thanks again
  pari_init(paristacksize, 2); // init PARI stack

  printf("inverse_mat_mod2: stack initialized\n");
//   std::cout << "inverse_mat_mod2: stack initialized" << std::cout;

  /* allocate memory for the temp matrix
   * note: remember that PARI matrices of mxn need (m+1)x(n+1), and
   * that elements go a[1], a[2], ..., while a[0] is for type
   * information */
  printf("inverse_mat_mod2: size = %d\n", _size);
//   std::cout << "inverse_mat_mod2: size = " << _size << std::endl;
  a = cgetg(sizepari, t_MAT); // `a' has _size columns
  for (i = 1; i <= _size; ++i) { // allocate `a' columns of _size elements
    a[i] = lgetg(sizepari, t_COL);
  }
  for (i = 1; i <= _size; ++i) { // allocate space for elements
    v = (GEN)a[i];
    for (j = 1; j <= _size; ++j) {
      v[j] = lgetg(3, t_INTMOD);
      z = (GEN)v[j];
      z[1] = lstoi(2); // modulo 2
      if (_a[j-1][i-1]) z[2] = un; else z[2] = zero;//lstoi(0);
    }
  }
  printf("inverse_mat_mod2: assigment done\n");
//   std::cout << "inverse_mat_mod2: assigment done" << std::endl;

  /* compute determinant to check invertibility */
  determinant = det(a);
  printf("Hi\n");
  if (determinant == gzero) {
    printf("determinant is zero\n");
  printf("By\n");
    freeall();
  printf("By\n");
    pari_sig_init(SIG_DFL);
  printf("By\n");
    return FALSE;
  }

  /* invert matrix */
  ainv = ginv(a);
  printf("inverse_mat_mod2: matrix inverted\n");
//   std::cout << "inverse_mat_mod2: matrix inverted" << std::endl;

  /* check whether the inverse is valid */
  if (_CHECK) {
    GEN a_ainv;
    a_ainv = gmul(a, ainv);
    for (i = 1; i <= _size; ++i) {
      v = (GEN)a_ainv[i];
      for (j = 1; j <= _size; ++j) {
        z = (GEN)v[j];
        if ((i == j) != itos((GEN)z[2])) {
          freeall();
          pari_sig_init(SIG_DFL);
          return FALSE;
        }
      }
    }
  }

  /* for no _CHECK or valid inversion, overwrite input array with
     inverse matrix */
  for (i = 1; i <= _size; ++i) {
    v = (GEN)ainv[i];
    for (j = 1; j <= _size; ++j) {
      z = (GEN)v[j];
      _a[j-1][i-1] = (mybool)itos((GEN)z[2]);
    }
  }

  /* free PARI's memory (thanks Karim) */
  freeall();
  pari_sig_init(SIG_DFL);

  return TRUE;
}

/* print_c_myboolmatrix: print a C myboolean matrix */
void print_c_boolmatrix(mybool **_a, unsigned _m, unsigned _n) {

  unsigned m, n;

  printf("{\n");
  for (m = 0; m < _m; ++m) {
    printf("{");
    for (n = 0; n < _n; ++n) {
      printf("%d ", _a[m][n]);
      if (n == _n - 1) {
        printf("}\n");
      }
      else {
        printf(" ");
      }
    }
  }
  printf("}\n");

  return;
}

/* XOR function */
/* mybool exor(mybool a, mybool b) { return (mybool)(a != b); } */

----------------------------------------------------------------------

And the included headers file just have

-------------------------------------------------------------------

#ifndef MODULO2_H
#define MODULO2_H

typedef enum{FALSE, TRUE} mybool;

/* inverse_mat_mod2: overwrites a matrix modulo 2 with its inverse;
 * input is a bool matrix, but internally, PARI library types are
 * used */
/* mybool inverse_mat_mod2(mybool **_a, unsigned _size); */
mybool inverse_mat_mod2(mybool **_a, unsigned _size, mybool _CHECK);
/* print_c_boolmatrix: print a C boolean matrix */
void print_c_boolmatrix(mybool **_a, unsigned _m, unsigned _n);
/* logic XOR function */
/* mybool exor(mybool a, mybool b); */
#define exor(a, b) ((mybool)(a != b))

#endif

--------------------------------------------------------------------

regards,

Ramón.

-- 
+++ GMX - Mail, Messaging & more  http://www.gmx.net +++
NEU: Mit GMX ins Internet. Rund um die Uhr für 1 ct/ Min. surfen!