Actually added the files mentioned in previous commit just now.

Also added an example file celldriver.cpp and transformer.geo
in /utils/misc

Geo/CellComplex.cpp

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka, 16.3.2009.
+
+#include "CellComplex.h"
+
+CellComplex::CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain ){
+  
+  _domain = domain;
+  _subdomain = subdomain;
+  
+  // subdomain need to be inserted first!
+  insertCells(true);
+  insertCells(false);
+  
+  
+}
+void CellComplex::insertCells(bool subdomain){  
+  
+  std::vector<GEntity*> domain;
+  
+  if(subdomain) domain = _subdomain;
+  else domain = _domain;
+  
+  std::vector<int> vertices;
+  int vertex = 0;
+  
+  for(unsigned int j=0; j < domain.size(); j++) {  
+    for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){ 
+      vertices.clear();
+      for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){              
+        MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
+        vertices.push_back(vertex->getNum()); 
+        //_cells[0].insert(new ZeroSimplex(vertex->getNum(), subdomain, vertex->x(), vertex->y(), vertex->z() )); // Add vertices
+        _cells[0].insert(new ZeroSimplex(vertex->getNum(), subdomain));
+      } 
+      if(domain.at(j)->getMeshElement(i)->getDim() == 3){
+        _cells[3].insert(new ThreeSimplex(vertices, subdomain)); // Add volumes
+      }
+      
+      for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumFaces(); k++){
+        vertices.clear();
+        for(int l=0; l < domain.at(j)->getMeshElement(i)->getFace(k).getNumVertices(); l++){
+          vertex = domain.at(j)->getMeshElement(i)->getFace(k).getVertex(l)->getNum();
+          vertices.push_back(vertex); 
+        } 
+        _cells[2].insert(new TwoSimplex(vertices, subdomain)); // Add faces
+      }
+      
+      for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumEdges(); k++){
+        vertices.clear();
+        for(int l=0; l < domain.at(j)->getMeshElement(i)->getEdge(k).getNumVertices(); l++){
+          vertex = domain.at(j)->getMeshElement(i)->getEdge(k).getVertex(l)->getNum();
+          vertices.push_back(vertex); 
+        }
+        _cells[1].insert(new OneSimplex(vertices, subdomain)); // Add edges
+      }
+              
+    }               
+  }
+  
+}
+
+int Simplex::kappa(Cell* tau) const{
+  for(int i=0; i < tau->getNumVertices(); i++){
+    if( !(this->hasVertex(tau->getVertex(i))) ) return 0;
+  }
+  
+  if(this->getDim() - tau->getDim() != 1) return 0;
+  
+  int value=1;
+  for(int i=0; i < tau->getNumVertices(); i++){
+    if(this->getVertex(i) != tau->getVertex(i)) return value;
+    value = value*-1;
+  }
+  
+  return 1;  
+}
+std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, std::vector<int>& vertices){
+  if(dim == 0) return _cells[dim].find(new ZeroSimplex(vertices.at(0)));
+  if(dim == 1) return _cells[dim].find(new OneSimplex(vertices));
+  if(dim == 2) return _cells[dim].find(new TwoSimplex(vertices));
+  return _cells[3].find(new ThreeSimplex(vertices));
+}
+
+std::set<Cell*, Less_Cell>::iterator CellComplex::findCell(int dim, int vertex, int dummy){
+  if(dim == 0) return _cells[dim].lower_bound(new ZeroSimplex(vertex));
+  if(dim == 1) return _cells[dim].lower_bound(new OneSimplex(vertex, dummy));
+  if(dim == 2) return _cells[dim].lower_bound(new TwoSimplex(vertex, dummy));
+  return _cells[3].lower_bound(new ThreeSimplex(vertex, dummy));
+}
+
+
+std::vector<Cell*> CellComplex::bd(Cell* sigma){  
+  std::vector<Cell*> boundary;
+  int dim = sigma->getDim();
+  if(dim < 1) return boundary;
+  
+  
+  citer start = findCell(dim-1, sigma->getVertex(0), 0);
+  if(start == lastCell(dim-1)) return boundary;
+  
+  citer end = findCell(dim-1, sigma->getVertex(sigma->getNumVertices()-1), sigma->getVertex(sigma->getNumVertices()-1));
+  if(end != lastCell(dim-1)) end++;
+
+  //for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
+  for(citer cit = start; cit != end; cit++){
+    if(kappa(sigma, *cit) != 0){
+      boundary.push_back(*cit);
+      if(boundary.size() == sigma->getBdMaxSize()){
+        return boundary;
+      }
+    }
+  }
+  return boundary;    
+}
+std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::bdIt(Cell* sigma){
+  
+  int dim = sigma->getDim();  
+  std::vector< std::set<Cell*, Less_Cell>::iterator >  boundary;
+  if(dim < 1){
+    return boundary;
+  }
+  
+  citer start = findCell(dim-1, sigma->getVertex(0), 0);
+  if(start == lastCell(dim-1)) return boundary;
+  
+  citer end = findCell(dim-1, sigma->getVertex(sigma->getNumVertices()-1), sigma->getVertex(sigma->getNumVertices()-1));
+  if(end != lastCell(dim-1)) end++;
+    
+  for(citer cit = start; cit != end; cit++){
+    if(kappa(sigma, *cit) != 0){
+      boundary.push_back(cit);
+      if(boundary.size() == sigma->getBdMaxSize()){
+        return boundary;
+      }
+    }
+  }
+  
+  return boundary;
+}
+
+
+
+std::vector<Cell*> CellComplex::cbd(Cell* tau){  
+  std::vector<Cell*> coboundary;
+  int dim = tau->getDim();
+  if(dim > 2){
+    return coboundary;
+  }
+  
+  for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
+    if(kappa(*cit, tau) != 0){
+      coboundary.push_back(*cit);
+      if(coboundary.size() == tau->getCbdMaxSize()){
+        return coboundary;
+      }
+    }
+  }
+  return coboundary;    
+}
+std::vector< std::set<Cell*, Less_Cell>::iterator > CellComplex::cbdIt(Cell* tau){
+  
+  std::vector< std::set<Cell*, Less_Cell>::iterator >  coboundary;
+  int dim = tau->getDim();
+  if(dim > 2){
+    return coboundary;
+  }
+  
+  for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
+    if(kappa(*cit, tau) != 0){
+      coboundary.push_back(cit);
+      if(coboundary.size() == tau->getCbdMaxSize()){
+        return coboundary;
+      }
+    }
+  }
+  
+  return coboundary;
+}
+
+
+void CellComplex::removeCell(Cell* cell){
+    _cells[cell->getDim()].erase(cell);
+  delete cell;
+}
+void CellComplex::removeCellIt(std::set<Cell*, Less_Cell>::iterator cell){
+  Cell* c = *cell;
+  int dim = c->getDim();
+  _cells[dim].erase(cell);
+  delete c;
+}
+
+void CellComplex::removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset){
+   Qset.erase(cell);
+  delete cell;
+}
+
+void CellComplex::enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells, 
+                               std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
+  for(unsigned int i=0; i < cells.size(); i++){
+    Cell* cell = *cells.at(i);
+    citer cit = Qset.find(cell);
+    if(cit == Qset.end()){
+      Qset.insert(cell);
+      Q.push(cell);
+    } 
+  }
+}
+void CellComplex::enqueueCells(std::vector<Cell*>& cells, std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset){
+  for(unsigned int i=0; i < cells.size(); i++){
+    citer cit = Qset.find(cells.at(i));
+    if(cit == Qset.end()){
+      Qset.insert(cells.at(i));
+      Q.push(cells.at(i));
+    }
+  }
+}
+
+
+int CellComplex::coreductionMrozek(Cell* generator){
+  
+  int coreductions = 0;
+  
+  std::queue<Cell*> Q;
+  std::set<Cell*, Less_Cell> Qset;
+  
+  Q.push(generator);
+  Qset.insert(generator);
+  //removeCell(generator);
+  
+  std::vector<Cell*> bd_s;
+  //std::vector< std::set<Cell*, Less_Cell>::iterator > bd_s;
+  std::vector<Cell*> cbd_c;
+  //std::vector< std::set<Cell*, Less_Cell>::iterator > cbd_c;
+  Cell* s;
+  int round = 0;
+  while( !Q.empty() ){
+    round++;
+    //printf("%d ", round);
+    s = Q.front();
+    Q.pop();
+    removeCellQset(s, Qset);
+    
+    //bd_s = bdIt(s);
+    bd_s = bd(s);
+    if( bd_s.size() == 1 && inSameDomain(s, bd_s.at(0)) ){
+      removeCell(s);
+      //cbd_c = cbdIt(*bd_s.at(0));
+      cbd_c = cbd(bd_s.at(0));
+      //enqueueCellsIt(cbd_c, Q, Qset);
+      enqueueCells(cbd_c, Q, Qset);
+      //removeCellIt(bd_s.at(0));
+      removeCell(bd_s.at(0));
+      coreductions++;
+    }
+    else if(bd_s.empty()){
+      //cbd_c = cbdIt(s);
+      //enqueueCellsIt(cbd_c, Q, Qset);
+      cbd_c = cbd(s);
+      enqueueCells(cbd_c, Q, Qset);
+    }
+    
+  
+  }
+  printf("Coreduction: %d loops with %d coreductions\n", round, coreductions);
+  return coreductions;
+}
+
+
+int CellComplex::coreduction(int dim){
+  
+  if(dim < 0 || dim > 2) return 0;
+  
+  std::vector<Cell*> bd_c;
+  int count = 0;
+  
+  bool coreduced = true;
+  while (coreduced){
+    coreduced = false;
+    for(citer cit = firstCell(dim+1); cit != lastCell(dim+1); cit++){
+      Cell* cell = *cit;
+      
+      bd_c = bd(cell);
+      if(bd_c.size() == 1 && inSameDomain(cell, bd_c.at(0)) ){
+        removeCell(cell);
+        //cit = firstCell(dim-1);
+        removeCell(bd_c.at(0));
+        count++;
+        coreduced =true;
+      }
+      
+    }
+  }
+  
+  return count;
+    
+}
+
+
+int CellComplex::reduction(int dim){
+  if(dim < 1 || dim > 3) return 0;
+  
+  std::vector<Cell*> cbd_c;
+  int count = 0;
+  
+  bool reduced = true;
+  while (reduced){
+    reduced = false;
+    for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
+      Cell* cell = *cit;
+      //if(reductionMrozek(cell) !=0) cit = firstCell(dim-1);   
+      cbd_c = cbd(cell);
+      if(cbd_c.size() == 1 && inSameDomain(cell, cbd_c.at(0)) ){
+        removeCell(cell);
+        //cit = firstCell(dim-1);
+        removeCell(cbd_c.at(0));
+        count++;
+        reduced =true;
+      }
+      
+    }
+  }
+  return count;  
+}
+
+void CellComplex::reduceComplex(){
+  printf("Cell complex before reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+         getSize(3), getSize(2), getSize(1), getSize(0));
+  reduction(3);
+  reduction(2);
+  reduction(1);
+  printf("Cell complex after reduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+         getSize(3), getSize(2), getSize(1), getSize(0));
+}
+void CellComplex::coreduceComplex(){
+  printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+         getSize(3), getSize(2), getSize(1), getSize(0));
+  while (getSize(0) != 0){
+    citer cit = firstCell(0);
+    Cell* cell = *cit;
+    removeCell(cell);
+    //complex.coreductionMrozek(kaka);
+    coreduction(0);
+    coreduction(1);
+    coreduction(2);
+    
+   }
+  
+  printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+         getSize(3), getSize(2), getSize(1), getSize(0));
+}
+
+void CellComplex::constructHMatrix(int dim){
+  
+  // h: C_dim -> C_(dim-1)
+  
+  if(dim > 3 || dim < 1){
+    return;
+  }
+  
+  if( _cells[dim].size() == 0 ){
+    _HMatrix[dim] = create_gmp_matrix_zero(1, 1);
+    return;
+  }
+  unsigned int cols = _cells[dim].size(); 
+  
+  if( _cells[dim-1].size() == 0){ 
+    _HMatrix[dim] = create_gmp_matrix_zero(1, cols);
+    return;
+  }
+  unsigned int rows = _cells[dim-1].size();
+  
+  mpz_t elems[rows*cols];
+  mpz_array_init( elems[0], rows*cols, 512 );  
+  
+  citer high = firstCell(dim);
+  citer low = firstCell(dim-1);
+  
+  //printf("laa %d %d %d \n", rows, cols, rows*cols);
+  for(unsigned int i=0; i < rows*cols; i++){
+    //printf("%d, ", i);
+    if(low == lastCell(dim-1)){
+      //printf("rowfull %d", i);
+      high++;
+      low = firstCell(dim-1);
+    }
+    mpz_set_ui(elems[i], kappa(*high, *low));
+    low++;
+  }
+  
+  _HMatrix[dim] = create_gmp_matrix(rows, cols, elems);
+
+  return; 
+}
+
+
+void CellComplex::getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, bool subdomain){
+
+  std::vector<GEntity*> domain;
+  if(subdomain) domain = _subdomain;
+  else domain = _domain;
+  
+  for(unsigned int j=0; j < domain.size(); j++) {
+    for(unsigned int i=0; i < domain.at(j)->getNumMeshElements(); i++){
+      for(int k=0; k < domain.at(j)->getMeshElement(i)->getNumVertices(); k++){
+        MVertex* vertex = domain.at(j)->getMeshElement(i)->getVertex(k);
+        domainVertices.insert(vertex);
+      }
+    }
+  }
+      
+/*  for(unsigned int j=0; j < _domain.size(); j++) { 
+    for(unsigned int i=0; i < _domain.at(j)->getNumMeshVertices(); i++){
+      MVertex* vertex =  _domain.at(j)->getMeshVertex(i);
+      domainVertices.insert(vertex);
+    }
+    
+    std::list<GFace*> faces = _domain.at(j)->faces();
+    for(std::list<GFace*>::iterator fit = faces.begin(); fit != faces.end(); fit++){
+      GFace* face = *fit;
+      for(unsigned int i=0; i < face->getNumMeshVertices(); i++){
+        MVertex* vertex =  face->getMeshVertex(i);
+        domainVertices.insert(vertex);
+      }
+    }
+    std::list<GEdge*> edges = _domain.at(j)->edges();
+    for(std::list<GEdge*>::iterator eit = edges.begin(); eit != edges.end(); eit++){
+      GEdge* edge = *eit;
+      for(unsigned int i=0; i < edge->getNumMeshVertices(); i++){
+        MVertex* vertex =  edge->getMeshVertex(i);
+        domainVertices.insert(vertex);
+      }
+    }
+    std::list<GVertex*> vertices = _domain.at(j)->vertices();
+    for(std::list<GVertex*>::iterator vit = vertices.begin(); vit != vertices.end(); vit++){
+      GVertex* vertex = *vit;
+      for(unsigned int i=0; i < vertex->getNumMeshVertices(); i++){
+        MVertex* mvertex =  vertex->getMeshVertex(i);
+        domainVertices.insert(mvertex);
+      }
+    }
+    
+  }
+ */
+  return;
+}
+
+int CellComplex::writeComplexMSH(const std::string &name){
+  
+    
+  FILE *fp = fopen(name.c_str(), "w");
+  if(!fp){
+    Msg::Error("Unable to open file '%s'", name.c_str());
+    return 0;
+  }
+  
+  fprintf(fp, "$NOD\n");
+  
+  std::set<MVertex*, Less_MVertex> domainVertices;
+  getDomainVertices(domainVertices, true);
+  getDomainVertices(domainVertices, false);
+  
+  fprintf(fp, "%d\n", domainVertices.size());
+  
+  for(std::set<MVertex*, Less_MVertex>::iterator vit = domainVertices.begin(); vit != domainVertices.end(); vit++){
+    MVertex* vertex = *vit;
+    fprintf(fp, "%d %.16g %.16g %.16g\n", vertex->getNum(), vertex->x(), vertex->y(), vertex->z() );
+  }
+  
+      
+  fprintf(fp, "$ENDNOD\n");
+  fprintf(fp, "$ELM\n");
+
+  fprintf(fp, "%d\n", _cells[0].size() + _cells[1].size() + _cells[2].size() + _cells[3].size());
+  
+  int index = 1;
+  
+  for(citer cit = firstCell(0); cit != lastCell(0); cit++) {
+    Cell* vertex = *cit;
+    fprintf(fp, "%d %d %d %d %d %d\n", index, 15, 0, 1, 1, vertex->getVertex(0));
+    index++;
+  }
+  
+  
+  for(citer cit = firstCell(1); cit != lastCell(1); cit++) {
+    Cell* edge = *cit;
+    fprintf(fp, "%d %d %d %d %d %d %d\n", index, 1, 0, 1, 2, edge->getVertex(0), edge->getVertex(1));
+    index++;
+  }
+  
+  for(citer cit = firstCell(2); cit != lastCell(2); cit++) {
+    Cell* face = *cit;
+    fprintf(fp, "%d %d %d %d %d %d %d %d\n", index, 2, 0, 1, 3, face->getVertex(0), face->getVertex(1), face->getVertex(2));
+    index++;
+  }
+  for(citer cit = firstCell(3); cit != lastCell(3); cit++) {
+    Cell* volume = *cit;
+    fprintf(fp, "%d %d %d %d %d %d %d %d %d\n", index, 4, 0, 1, 4, volume->getVertex(0), volume->getVertex(1), volume->getVertex(2), volume->getVertex(3));
+    index++;
+  }
+    
+  fprintf(fp, "$ENDELM\n");
+  
+  return 1;
+}
+
+
+void CellComplex::printComplex(int dim){
+  for (citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+    Cell* cell = *cit;
+    for(int i = 0; i < cell->getNumVertices(); i++){
+      printf("%d ", cell->getVertex(i));
+    }
+    printf("\n");
+  }
+}

Geo/CellComplex.h

+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka, 16.3.2009.
+
+#ifndef _CELLCOMPLEX_H_
+#define _CELLCOMPLEX_H_
+
+#include <stdio.h>
+#include <string>
+#include <algorithm>
+#include <set>
+#include <queue>
+#include "MElement.h"
+#include "GModel.h"
+#include "GEntity.h"
+#include "GRegion.h"
+#include "GFace.h"
+#include "GVertex.h"
+#include "gmp.h"
+extern "C" {
+  #include "gmp_normal_form.h" // perhaps make c++ headers instead?
+}
+
+// Abstract class representing a cell of a cell complex.
+class Cell
+{  
+  protected:
+   
+   // cell dimension
+   int _dim;
+   // maximum number of lower dimensional cells on the boundary of this cell
+   int _bdMaxSize;
+   // maximum number of higher dimensional cells on the coboundary of this cell
+   int _cbdMaxSize;
+
+   int _bdSize;
+   int _cbdSize;
+   
+   // whether this cell belongs to a subdomain
+   // used to in relative homology computation
+   bool _inSubdomain;
+   
+  public:
+   Cell(){}
+   virtual ~Cell(){}
+   
+   virtual int getDim() const { return _dim; };
+   //virtual int getTag() const { return _tag; };
+   
+   // get the number of vertices this cell has
+   virtual int getNumVertices() const = 0; //{return _vertices.size();}
+   
+   virtual int getVertex(int vertex) const = 0; //{return _vertices.at(vertex);}
+   virtual std::vector<int> getVertexVector() const = 0;
+   
+   // returns 1 or -1 if lower dimensional cell tau is on the boundary of this cell
+   // otherwise returns 0
+   // implementation will vary depending on cell type
+   virtual int kappa(Cell* tau) const = 0;
+   
+   // true if this cell has given vertex
+   virtual bool hasVertex(int vertex) const = 0;
+   
+   virtual unsigned int getBdMaxSize() const { return _bdMaxSize; };
+   virtual unsigned int getCbdMaxSize() const { return _cbdMaxSize; };
+
+   virtual int getBdSize() { return _bdSize; }
+   virtual void setBdSize(int size) { _bdSize = size; }
+   virtual int getCbdSize() { return _cbdSize; }
+   virtual void setCbdSize(int size) { _cbdSize = size; }
+   
+   virtual bool inSubdomain() const { return _inSubdomain; }
+   
+   // print cell vertices
+   virtual void printCell() const = 0;
+   
+   // return the coordinates of this cell (only used on zero-cells if anywhere)
+   virtual inline double x() const { return 0; }
+   virtual inline double y() const { return 0; }
+   virtual inline double z() const { return 0; }
+   
+   bool operator==(const Cell* c2){
+     if(this->getNumVertices() != c2->getNumVertices()){
+       return false;
+     }
+     for(int i=0; i < this->getNumVertices();i++){
+       if(this->getVertex(i) != c2->getVertex(i)){
+         return false;
+       }
+     }
+     return true;
+   }
+   
+   bool operator<(const Cell* c2) const {
+     if(this->getNumVertices() != c2->getNumVertices()){
+       return (this->getNumVertices() < c2->getNumVertices());
+     }
+     for(int i=0; i < this->getNumVertices();i++){
+       if(this->getVertex(i) < c2->getVertex(i)){
+         return true;
+       }
+     }
+     return false;
+   }
+   
+   
+   
+};
+
+// Simplicial cell type.
+class Simplex : public Cell
+{ 
+   
+ public:
+   Simplex(){
+     _cbdSize = 1000; // big number
+     _bdSize = 1000;
+   }
+   virtual ~Simplex(){}  
+  
+   // kappa for simplices
+   // checks whether lower dimensional simplex tau is on the boundary of this cell
+   virtual int kappa(Cell* tau) const; 
+   
+};
+
+// Zero simplex cell type.
+class ZeroSimplex : public Simplex
+{
+ private:
+   // number of the vertex
+   // same as the corresponding vertex in the finite element mesh
+   int _v;
+   // coordinates of this zero simplex, if needed
+   double _x, _y, _z;
+ public:
+   
+   ZeroSimplex(int vertex, bool subdomain=false, double x=0, double y=0, double z=0) : Simplex(){
+     _v = vertex;
+     _dim = 0;
+     _bdMaxSize = 0;
+     _cbdMaxSize = 1000; // big number
+     _x = x;
+     _y = y;
+     _z = z;
+     _inSubdomain = subdomain;
+   }
+   
+   virtual ~ZeroSimplex(){}
+   
+   virtual int getNumVertices() const { return 1; }
+   virtual int getVertex(int vertex) const {return _v; }
+   virtual bool hasVertex(int vertex) const {return (_v == vertex); }
+   
+   virtual std::vector<int> getVertexVector() const { return std::vector<int>(1,_v); }
+
+   virtual inline double x() const { return _x; }
+   virtual inline double y() const { return _y; }
+   virtual inline double z() const { return _z; }
+   
+   virtual void printCell()const { printf("Vertices: %d\n", _v); }
+   
+};
+
+// One simplex cell type.
+class OneSimplex : public Simplex
+{
+  private:
+   // numbers of the vertices of this simplex
+   // same as the corresponding vertices in the finite element mesh
+   int _v[2];
+  public:
+   
+   OneSimplex(std::vector<int> vertices, bool subdomain=false) : Simplex(){
+     sort(vertices.begin(), vertices.end());
+     _v[0] = vertices.at(0);
+     _v[1] = vertices.at(1);
+     _dim = 1;
+     _bdMaxSize = 2;
+     _cbdMaxSize = 1000;
+     _inSubdomain = subdomain;
+   }
+   
+   // constructor for the dummy one simplex
+   // used to find another definite one simplex from the cell complex
+   // first vertex gives the lower bound from where to look
+   OneSimplex(int vertex, int dummy){
+     _v[0] = vertex;
+     _v[1] = dummy;
+   }
+   
+   virtual ~OneSimplex(){}
+   
+   virtual int getNumVertices() const { return 2; }
+   virtual int getVertex(int vertex) const {return _v[vertex]; }
+   virtual bool hasVertex(int vertex) const {return (_v[0] == vertex || _v[1] == vertex); }
+   
+   virtual std::vector<int> getVertexVector() const { 
+     return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+   
+   virtual void printCell() const { printf("Vertices: %d %d\n", _v[0], _v[1]); }
+};
+
+// Two simplex cell type.
+class TwoSimplex : public Simplex
+{
+  private:
+   // numbers of the vertices of this simplex
+   // same as the corresponding vertices in the finite element mesh
+   int _v[3];
+  public:
+   
+   TwoSimplex(std::vector<int> vertices, bool subdomain=false) : Simplex(){
+     sort(vertices.begin(), vertices.end());
+     _v[0] = vertices.at(0);
+     _v[1] = vertices.at(1);
+     _v[2] = vertices.at(2);
+     _dim = 2;
+     _bdMaxSize = 3;
+     _cbdMaxSize = 2;
+     _inSubdomain = subdomain;
+   }
+   // constructor for the dummy one simplex
+   TwoSimplex(int vertex, int dummy){
+     _v[0] = vertex;
+     _v[1] = dummy;
+     _v[2] = dummy;
+   }
+   
+   virtual ~TwoSimplex(){}
+   
+   virtual int getNumVertices() const { return 3; }
+   virtual int getVertex(int vertex) const {return _v[vertex]; }
+   virtual bool hasVertex(int vertex) const {return 
+       (_v[0] == vertex || _v[1] == vertex || _v[2] == vertex); }
+   virtual std::vector<int> getVertexVector() const { 
+     return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+   
+   virtual void printCell() const { printf("Vertices: %d %d %d\n", _v[0], _v[1], _v[2]); }
+};
+
+// Three simplex cell type.
+class ThreeSimplex : public Simplex
+{
+  private:
+   // numbers of the vertices of this simplex
+   // same as the corresponding vertices in the finite element mesh
+   int _v[4];
+  public:
+   
+   ThreeSimplex(std::vector<int> vertices, bool subdomain=false) : Simplex(){
+     sort(vertices.begin(), vertices.end());
+     _v[0] = vertices.at(0);
+     _v[1] = vertices.at(1);
+     _v[2] = vertices.at(2);
+     _v[3] = vertices.at(3);
+     _dim = 3;
+     _bdMaxSize = 4;
+     _cbdMaxSize = 0;
+     _inSubdomain = subdomain;
+   }
+   // constructor for the dummy one simplex
+   ThreeSimplex(int vertex, int dummy){
+     _v[0] = vertex;
+     _v[1] = dummy;
+     _v[2] = dummy;
+     _v[3] = dummy;
+   }
+   
+   virtual ~ThreeSimplex(){}
+   
+   virtual int getNumVertices() const { return 4; }
+   virtual int getVertex(int vertex) const {return _v[vertex]; }
+   virtual bool hasVertex(int vertex) const {return 
+       (_v[0] == vertex || _v[1] == vertex || _v[2] == vertex || _v[3] == vertex); }
+   virtual std::vector<int> getVertexVector() const { 
+     return std::vector<int>(_v, _v + sizeof(_v)/sizeof(int) ); }
+   
+   virtual void printCell() const { printf("Vertices: %d %d %d %d\n", _v[0], _v[1], _v[2], _v[3]); }
+};
+
+// Ordering for the cells.
+class Less_Cell{
+  public:
+   bool operator()(const Cell* c1, const Cell* c2) const {
+     if(c1->getNumVertices() != c2->getNumVertices()){
+       return (c1->getNumVertices() < c2->getNumVertices());
+     }
+     for(int i=0; i < c1->getNumVertices();i++){
+       if(c1->getVertex(i) < c2->getVertex(i)){
+         return true;
+       }
+       else if (c1->getVertex(i) > c2->getVertex(i)){
+         return false;
+       }
+     }
+     return false;
+   }
+};
+
+// Ordering for the finite element mesh vertices
+class Less_MVertex{
+  public:
+   bool operator()(const MVertex* v1, const MVertex* v2) const {
+     return (v1->getNum() < v2->getNum());
+   }
+};
+
+// A class representing a cell complex made out of a finite element mesh.
+class CellComplex
+{
+ private:
+   
+   // the domain in the model which this cell complex covers
+   std::vector<GEntity*> _domain;
+   // a subdomain of the given domain
+   // used in relative homology computation, may be empty
+   std::vector<GEntity*> _subdomain;
+   
+   // sorted container of unique cells in this cell complex 
+   // mapped according to their dimension
+   std::map< int, std::set<Cell*, Less_Cell> >  _cells;
+   
+   // boundary operator matrices for this cell complex
+   // h_k: C_k -> C_(k-1)
+   std::map<int, gmp_matrix*>  _HMatrix; 
+   
+   
+ public:
+   // iterator for the cells of same dimension
+   typedef std::set<Cell*, Less_Cell>::iterator citer;
+   
+  protected: 
+   // enqueue cells in queue if they are not there already
+   virtual void enqueueCells(std::vector<Cell*>& cells, 
+                             std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
+   virtual void enqueueCellsIt(std::vector< std::set<Cell*, Less_Cell>::iterator >& cells,
+                               std::queue<Cell*>& Q, std::set<Cell*, Less_Cell>& Qset);
+     
+   // remove cell from the queue set
+   virtual void removeCellQset(Cell*& cell, std::set<Cell*, Less_Cell>& Qset);
+   
+   // get all the finite element mesh vertices in the domain of this cell complex
+   virtual void getDomainVertices(std::set<MVertex*, Less_MVertex>& domainVertices, 
+                                  bool subdomain);
+   
+   // insert cells into this cell complex
+   virtual void insertCells(bool subdomain);
+   
+  public: 
+   CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
+   virtual ~CellComplex(){}
+   
+   // get the number of certain dimensional cells
+   virtual int getSize(int dim){ return _cells[dim].size(); }
+   
+   // iterators to the first and last cells of certain dimension
+   virtual citer firstCell(int dim) {return _cells[dim].begin(); }
+   virtual citer lastCell(int dim) {return _cells[dim].end(); }
+  
+   // find a cell in this cell complex
+   virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, std::vector<int>& vertices);
+   virtual std::set<Cell*, Less_Cell>::iterator findCell(int dim, int vertex, int dummy=0);
+   
+   
+   // kappa for two cells of this cell complex
+   // implementation will vary depending on cell type
+   virtual int kappa(Cell* sigma, Cell* tau) const { return sigma->kappa(tau); }
+   
+   // cells on the boundary of a cell
+   virtual std::vector<Cell*> bd(Cell* sigma);
+   virtual std::vector< std::set<Cell*, Less_Cell>::iterator > bdIt(Cell* sigma);
+   // cells on the coboundary of a cell
+   virtual std::vector<Cell*> cbd(Cell* tau);
+   virtual std::vector< std::set<Cell*, Less_Cell>::iterator > cbdIt(Cell* tau);
+   
+   // remove a cell from this cell complex
+   virtual void removeCell(Cell* cell);
+   virtual void removeCellIt(std::set<Cell*, Less_Cell>::iterator cell);
+   
+   // check whether two cells both belong to subdomain or if neither one does
+   virtual bool inSameDomain(Cell* c1, Cell* c2) const { return 
+       ( (!c1->inSubdomain() && !c2->inSubdomain()) || (c1->inSubdomain() && c2->inSubdomain()) ); }
+   
+   // coreduction of this cell complex
+   // removes corection pairs of cells of dimension dim and dim+1
+   virtual int coreduction(int dim);
+   // reduction of this cell complex
+   // removes reduction pairs of cell of dimension dim and dim-1
+   virtual int reduction(int dim);
+   
+   // useful functions for (co)reduction of cell complex
+   virtual void reduceComplex();
+   virtual void coreduceComplex();
+   
+   // queued coreduction presented in Mrozek's paper
+   // slower, but produces cleaner result
+   virtual int coreductionMrozek(Cell* generator);
+   
+   
+   // construct boundary operator matrix dim->dim-1
+   virtual void constructHMatrix(int dim);
+   // get the boundary operator matrix dim->dim-1
+   virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
+   
+   // print the vertices of cells of certain dimension
+   virtual void printComplex(int dim);
+   
+   // write this cell complex in legacy .msh format
+   virtual int writeComplexMSH(const std::string &name); 
+   
+};
+
+
+
+#endif

contrib/kbipack/Makefile

+## Original Makefile by Saku Suuriniemi
+
+#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
+
+## Makefile for Gmsh
+
+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

contrib/kbipack/Makefile.bak

+#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_matrix${OBJEXT}: gmp_matrix.c gmp_blas.h gmp_matrix.h
+gmp_blas${OBJEXT}: gmp_blas.c gmp_blas.h
+gmp_matrix_io${OBJEXT}: gmp_matrix_io.c gmp_matrix_io.h
+gmp_normal_form${OBJEXT}: gmp_normal_form.c gmp_blas.h gmp_matrix.h gmp_normal_form.h

contrib/kbipack/Makefile~

+#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

contrib/kbipack/compute_normal_form.c

+/* 
+   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;
+}

contrib/kbipack/gmp_blas.c

+/* 
+   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, const 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, const 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, 
+	     const mpz_t* x, size_t incx, 
+	     const 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, const 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, const 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(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, const 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;
+}

contrib/kbipack/gmp_blas.h

+/* 
+   Header 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.h,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+#ifndef __GMP_BLAS_H__
+#define __GMP_BLAS_H__
+
+#include <stdio.h>
+#include <gmp.h>
+
+/***********/
+/* Level 1 */
+/***********/
+
+/* x <-> y */
+void 
+gmp_blas_swap(size_t n,          mpz_t * x, size_t incx, mpz_t * y, size_t incy);
+
+/* x <- ax */
+void 
+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); 
+
+/* y <- ax + y */
+void 
+gmp_blas_axpy(size_t n, mpz_t a, 
+	      const 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); 
+
+
+/* Givens rotations are obviously impossible. However, the extended Euclid  
+   algorithm produces its close relative: 
+
+      g = gcd(a,b) & ka+lb = g 
+
+   E.g. in the computation of HNF, right-multipliations by the matrix 
+
+     [k -b/g]
+     [l  a/g]
+
+   occur frequently. Hence the name for the "Euclid rotation" function pair: */
+
+/* For a,b, find k,l,g such that ka+lb = g = gcd(a,b) holds */
+void gmp_blas_rotg(mpz_t g, mpz_t k, mpz_t l, mpz_t a, mpz_t b);
+#define gmp_blas_rotg mpz_gcdext
+
+
+/* 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); 
+
+
+/* In integer computations, the first nonzero element, minimal nonzero element 
+   by modulus, and maximal nonzero element by modulus come often handy: */
+
+/* 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);
+
+/* Similarly, x[(k-1)*incx] is maximal by modulus */
+size_t gmp_blas_iamax(size_t n, const 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);
+
+#endif

contrib/kbipack/gmp_matrix.c

+/* 
+   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"
+
+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;
+}

contrib/kbipack/gmp_matrix.c~

+/* 
+   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;
+}

contrib/kbipack/gmp_matrix.h

+/* 
+   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;
+
+/* 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

contrib/kbipack/gmp_matrix.h~

+/* 
+   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

contrib/kbipack/gmp_matrix_io.c

+/*
+  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, const 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;
+}

contrib/kbipack/gmp_matrix_io.h

+/*
+  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.h,v 1.1 2009-03-30 14:10:57 matti Exp $
+ */
+
+#ifndef _GMP_MATRIX_IO_
+#define _GMP_MATRIX_IO_
+
+
+#include "gmp_matrix.h"
+
+gmp_matrix * gmp_matrix_read_coord(char* filename);
+int          gmp_matrix_write_coord(char* filename, const gmp_matrix * M);
+
+
+#endif

contrib/kbipack/gmp_normal_form.c

+/* 
+   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;
+}
+
+

contrib/kbipack/gmp_normal_form.h

+/* 
+   Header file 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.2003Saku 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.h,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+#ifndef _GMP_NORMAL_FORM_
+#define _GMP_NORMAL_FORM_
+
+#include"gmp_blas.h"
+#include"gmp_matrix.h"
+
+typedef enum {NOT_INVERTED, INVERTED} inverted_flag;
+
+typedef struct
+{
+  gmp_matrix * left;
+  gmp_matrix * canonical;
+  gmp_matrix * right;
+  inverted_flag left_inverted;
+  inverted_flag right_inverted;
+} gmp_normal_form;
+
+/* For efficiency, the routines assume responsibility for the input matrices. 
+   Do *not* destroy them yourself! */
+
+/***********************/
+/* Hermite normal form */
+/***********************/
+/* 
+   PA = LU, 
+
+   left:
+   o P permutation matrix 
+     - Use "left_inverted = INVERTED" for left = P, 
+       and "left_inverted = NOT_INVERTED"  for left = P^T, i.e. for the 
+       decomposition 
+                      A = P^T L U.
+   canonical:
+   o L lower triangular, 
+     - diagonal entries positive, 
+     - subdiagonal entries positive 
+     - subdiagonal entries smaller than the diagonal entry of their row 
+
+   right:
+   o U unimodular,
+     - Use "right_inverted = NOT_INVERTED" for right = U, i.e. for the 
+       decomposition 
+                      P A = L U, 
+       and "right_inverted = INVERTED" for right = inv(U), i.e. for the 
+       decompositions 
+                      P A inv(U) = L and  A inv(U) = P^T L. 
+
+   Algorithm: Kannan-Bachem algorithm with improved (by Chou & Collins) 
+              reduction phase.
+
+   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.
+*/
+
+gmp_normal_form * 
+create_gmp_Hermite_normal_form(gmp_matrix * A, 
+			       inverted_flag left_inverted,
+			       inverted_flag right_inverted);
+
+
+
+/*********************/
+/* Smith normal form */
+/*********************/
+/* 
+   A = USV, 
+
+   left:
+   o U unimodular factor matrix 
+     - Use "left_inverted = NOT_INVERTED"  for left = U and 
+       "left_inverted = INVERTED" for left = inv(U), i.e. for 
+       the decomposition 
+                          inv(U) A = S V.
+
+   canonical:
+   o S diagonal, 
+     - first k diagonal entries positive (invariant factors of A), 
+       rest zero 
+     - each positive diagonal entry divisible by the previous one
+
+   right:
+   o V unimodular,
+     - Use "right_inverted = NOT_INVERTED" for right = V, 
+       and "right_inverted = INVERTED" for right = inv(V), i.e. for the 
+       decompositions 
+                      A inv(V) = U S and  inv(U) A inv(V) = S. 
+
+   Algorithm: Successive Hermite normal forms by Kannan-Bachem algorithm 
+              with improved (by Chou & Collins) reduction phase. See 
+	      reference [1].
+*/
+
+gmp_normal_form *
+create_gmp_Smith_normal_form(gmp_matrix * A,
+			     inverted_flag left_inverted,
+			     inverted_flag right_inverted);
+
+int destroy_gmp_normal_form(gmp_normal_form*);
+
+#endif