diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7b117e1cef3cb1c8af1000b2a7edceb390625a91
--- /dev/null
+++ b/Geo/CellComplex.cpp
@@ -0,0 +1,520 @@
+// 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");
+ }
+}
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
new file mode 100644
index 0000000000000000000000000000000000000000..55aa4bf4f035688d281e94128a0186820c12b7e4
--- /dev/null
+++ b/Geo/CellComplex.h
@@ -0,0 +1,419 @@
+// 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
diff --git a/contrib/kbipack/Makefile b/contrib/kbipack/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..50e6f52bae3bfea37f0d8d918ac9f2ebea24e9e5
--- /dev/null
+++ b/contrib/kbipack/Makefile
@@ -0,0 +1,77 @@
+## 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
diff --git a/contrib/kbipack/Makefile.bak b/contrib/kbipack/Makefile.bak
new file mode 100755
index 0000000000000000000000000000000000000000..edcef1ca84b0499b5cd1602e7271ea3a2aa12f5c
--- /dev/null
+++ b/contrib/kbipack/Makefile.bak
@@ -0,0 +1,77 @@
+#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
diff --git a/contrib/kbipack/Makefile~ b/contrib/kbipack/Makefile~
new file mode 100755
index 0000000000000000000000000000000000000000..3c361030c38340dcf7d3884a72c814f6c9f1ad52
--- /dev/null
+++ b/contrib/kbipack/Makefile~
@@ -0,0 +1,75 @@
+#FLAGS = -O3 -funroll-all-loops -ansi -pedantic -Wall -U DEBUG
+# Flags for different diagnostics:
+#FLAGS = -g -ansi -pedantic -Wall -fprofile-arcs -ftest-coverage -pg
+#FLAGS = -g -ansi -pedantic -Wall -pg
+
+# If gmp.h is not found de facto (e.g. when it is privately intalled),
+# its path should be put here
+#INCL =
+
+#LIBS = -lgmp
+#LIBDIR = .
+#SRC = gmp_normal_form.c gmp_matrix_io.c gmp_matrix.c gmp_blas.c
+
+# The compiler options are for gcc. If it is not available,
+# suitable flags must be set for the compiler
+#CC = gcc
+#RM = rm -f
+#RANLIB = ranlib
+
+#libkbi.a: $(SRC)
+# $(RM) *.o
+# $(CC) -c $(FLAGS) $(SRC)
+# ar r libkbi.a *.o
+# $(RM) *.o
+# $(RANLIB) libkbi.a
+
+#compute_normal_form: libkbi compute_normal_form.c
+# $(CC) $(FLAGS) -L $(LIBDIR) $(INCL) compute_normal_form.c -o compute_normal_form -lkbi $(LIBS)
+#clean:
+# $(RM) compute_normal_form
+# $(RM) *.o
+
+##
+
+include ../../variables
+
+LIB = ../../lib/libGmshKbi${LIBEXT}
+
+INC = ${DASH}I.
+
+CFLAGS = ${OPTIM} ${FLAGS} ${INC} ${SYSINCLUDE} -lgmp
+
+SRC = gmp_normal_form.c gmp_matrix_io.c gmp_matrix.c gmp_blas.c
+
+OBJ = ${SRC:.c=${OBJEXT}}
+
+.SUFFIXES: ${OBJEXT} .c
+
+${LIB}: ${OBJ}
+ ${AR} ${ARFLAGS}${LIB} ${OBJ}
+ ${RANLIB} ${LIB}
+
+cpobj: ${OBJ}
+ cp -f ${OBJ} ../../lib/
+
+.c${OBJEXT}:
+ ${CC} ${CFLAGS} ${DASH}c $<
+
+clean:
+ ${RM} *.o *.obj
+
+depend:
+ (sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+ ${CXX} -MM ${CFLAGS} ${SRC} | sed 's/.o:/$${OBJEXT}:/g' \
+ ) > Makefile.new
+ cp Makefile Makefile.bak
+ cp Makefile.new Makefile
+ rm -f Makefile.new
+
+# DO NOT DELETE THIS LINE
+gmp_normal_form${OBJEXT}: gmp_normal_form.c gmp_blas.h gmp_matrix.h \
+ gmp_normal_form.h
+gmp_matrix_io${OBJEXT}: gmp_matrix_io.c gmp_matrix_io.h gmp_matrix.h gmp_blas.h
+gmp_matrix${OBJEXT}: gmp_matrix.c gmp_matrix.h gmp_blas.h
+gmp_blas${OBJEXT}: gmp_blas.c gmp_blas.h
diff --git a/contrib/kbipack/compute_normal_form.c b/contrib/kbipack/compute_normal_form.c
new file mode 100755
index 0000000000000000000000000000000000000000..f86e043cc27e27144738acbffcbe7eafd6e7a6b3
--- /dev/null
+++ b/contrib/kbipack/compute_normal_form.c
@@ -0,0 +1,355 @@
+/*
+ COMPUTE_NORMAL_FORM
+
+ Hermite and Smith normal forms of integer matrices.
+
+ Implementation: Dense matrix with GMP-library's mpz_t elements to
+ hold huge integers.
+
+ Algorithm: Kannan - Bachem algorithm with improvement by
+ Chou and Collins. Expects a large number of unit invariant
+ factors and takes advantage of them as they appear.
+
+ References:
+ [1] Ravindran Kannan, Achim Bachem:
+ "Polynomial algorithms for computing the Smith and Hermite normal
+ forms of an integer matrix",
+ SIAM J. Comput., vol. 8, no. 5, pp. 499-507, 1979.
+ [2] Tsu-Wu J.Chou, George E. Collins:
+ "Algorithms for the solution of systems of linear Diophantine
+ equations",
+ SIAM J. Comput., vol. 11, no. 4, pp. 687-708, 1982.
+ [3] GMP homepage http://www.swox.com/gmp/
+ [4] GNU gmp page http://www.gnu.org/software/gmp/
+
+ Copyright (C) 3.11.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: compute_normal_form.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+
+#include<stdlib.h>
+#include<stdio.h>
+#include<string.h>
+#include"gmp_matrix_io.h"
+#include"gmp_normal_form.h"
+
+typedef enum {HERMITE, SMITH} nf_flag;
+typedef enum {NOT_WANTED, WANTED} w_flag;
+
+static int nf_parse_options(int arg_counter, char ** args,
+ char ** p_infile_name, char ** p_outfile_name,
+ nf_flag * p_normal_form_requested,
+ inverted_flag * p_left_inverted,
+ inverted_flag * p_right_inverted,
+ w_flag * p_help_wanted,
+ w_flag * p_verbose_wanted)
+{
+ size_t ind1, ind2;
+ size_t outfilename_length;
+
+ /* If no infile is specified, the error is indisputable.
+ => Print help and exit.
+ Infile should be the last argument. */
+
+ if(strncmp(args[arg_counter-1], "-", 1) == 0)
+ {
+ *(p_help_wanted) = WANTED;
+ return EXIT_SUCCESS;
+ }
+
+ /* Should I find arguments? */
+ if(arg_counter > 2)
+ {
+ for(ind1 = 1; ind1 < arg_counter-1; ind1++)
+ {
+ /* No multiple infile names ! */
+ if(strncmp(args[ind1], "-", 1) != 0)
+ {
+ *(p_help_wanted) = WANTED;
+ return EXIT_SUCCESS;
+ }
+
+ /* Parse the option letters after '-' one by one */
+ for(ind2 = 1; ind2 < strlen(args[ind1]); ind2++)
+ {
+ switch(args[ind1][ind2])
+ {
+ case 'h':
+ *(p_help_wanted) = WANTED;
+ return EXIT_SUCCESS;
+ break;
+ case 'H':
+ *(p_normal_form_requested) = HERMITE;
+ break;
+ case 'S':
+ *(p_normal_form_requested) = SMITH;
+ break;
+ case 'l':
+ *(p_left_inverted) = INVERTED;
+ break;
+ case 'r':
+ *(p_right_inverted) = INVERTED;
+ break;
+ case 'v':
+ *(p_verbose_wanted) = WANTED;
+ break;
+ }
+ }
+ }
+ }
+
+ /* Get the infile name */
+ * p_infile_name = (char *) calloc(strlen(args[arg_counter-1]) + 1,
+ sizeof(char));
+ if(* p_infile_name == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ strcpy(* p_infile_name, args[arg_counter-1]);
+
+ /* Create the outfile name */
+ outfilename_length = strcspn( * p_infile_name, ".");
+ * p_outfile_name = (char *) calloc(outfilename_length + 1,
+ sizeof(char));
+ if( * p_outfile_name == NULL)
+ {
+ free( * p_infile_name);
+ return EXIT_FAILURE;
+ }
+
+ strncpy(* p_outfile_name, * p_infile_name, outfilename_length);
+ (* p_outfile_name)[outfilename_length] = '\0';
+
+ return EXIT_SUCCESS;
+}
+
+static int nf_print_help(void)
+{
+ fprintf(stdout, "Usage: compute_normal_form [-hHSlrv] infile\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Description: Compute the Hermite/Smith normal form of an integer matrix.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Input: \"infile\" is an ascii file, which contains an integer matrix\n");
+ fprintf(stdout, " in the sparse coordinate format:\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, " # All comments before the data\n");
+ fprintf(stdout, " rows columns nonzeros\n");
+ fprintf(stdout, " row1 col1 value1\n");
+ fprintf(stdout, " row2 col2 value2\n");
+ fprintf(stdout, " ...\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, " No multiple infiles, please.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Options: -h This help\n");
+ fprintf(stdout, " -H Compute the Hermite normal form (default)\n");
+ fprintf(stdout, " -S Compute the Smith normal form\n");
+ fprintf(stdout, " -l Left factor inverted (A = P^(T)LU instead of PA = LU for Hermite, \n");
+ fprintf(stdout, " U^(-1)A = SV instead of A = USV for Smith)\n");
+ fprintf(stdout, " -r Right factor inverted (PAU^(-1) = L instead of PA = LU for Hermite, \n");
+ fprintf(stdout, " AV^(-1) = US instead of A = USV for Smith)\n");
+ fprintf(stdout, " -v Verbose. Print the matrix and the decomposition to sdtout, \n");
+ fprintf(stdout, " so that the user may check that the problem is read and written\n");
+ fprintf(stdout, " correctly.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Output: For inpur matrix in e.g. \"infile.coord\", writes ascii files\n");
+ fprintf(stdout, " \"infile.left\", \"infile.can\", and \"infile.right\", which \n");
+ fprintf(stdout, " contain the left factor, canonical form, and the right factor, \n");
+ fprintf(stdout, " respectively. The matrices are written in the sparse coordinate\n");
+ fprintf(stdout, " format.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Note: The return values may be big enough to not fit the standard integer types.\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Copyright (C) 3.11.2003 Saku Suuriniemi, Tampere University of Technology/Electromagnetics\n");
+ fprintf(stdout, " P.O.Box 692, FIN-33101 Tampere, Finland. saku.suuriniemi@tut.fi\n");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "Licenced under GNU General Public Licence (GPL) version 2 or later, ABSOLUTELY NO WARRANTY.\n");
+
+ return EXIT_SUCCESS;
+}
+
+
+static char* append_suffix(const char * file_base_name, const char * suffix)
+{
+ char * filename;
+
+ if((file_base_name == NULL) || (suffix == NULL))
+ {
+ return NULL;
+ }
+
+ filename =
+ (char*) calloc(strlen(file_base_name)+strlen(suffix)+2,
+ sizeof(char));
+
+ if(filename == NULL)
+ {
+ return NULL;
+ }
+
+ strcpy(filename, file_base_name);
+ strcat(filename, ".");
+ strcat(filename, suffix);
+
+ return filename;
+}
+
+
+
+int main(int argc, char * argv[])
+{
+ char * infile_name;
+ char * file_base_name;
+ char * outfile_name;
+ gmp_matrix * A;
+ gmp_normal_form * normal_form;
+ w_flag help_wanted;
+ w_flag verbose_wanted;
+ nf_flag normal_form_requested;
+ inverted_flag left_inverted;
+ inverted_flag right_inverted;
+
+ if(argc == 1)
+ {
+ fprintf(stdout, "Please try \"compute_normal_form -h\" for help.\n");
+ return EXIT_SUCCESS;
+ }
+
+ /* Default settings before option parsing */
+ normal_form_requested = HERMITE;
+ left_inverted = NOT_INVERTED;
+ right_inverted = NOT_INVERTED;
+ help_wanted = NOT_WANTED;
+ verbose_wanted = NOT_WANTED;
+
+ if( nf_parse_options(argc, argv,
+ &infile_name, &file_base_name,
+ &normal_form_requested,
+ &left_inverted, &right_inverted,
+ &help_wanted,
+ &verbose_wanted)
+ == EXIT_FAILURE )
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* In this case, no dynamical memory has been allocated by
+ nf_parse_options ! */
+ if(help_wanted == WANTED)
+ {
+ return nf_print_help();
+ }
+
+ if(verbose_wanted == WANTED)
+ {
+ fprintf(stdout, "Computing ");
+ if(normal_form_requested == HERMITE)
+ {
+ fprintf(stdout, "Hermite");
+ }
+ else
+ {
+ fprintf(stdout, "Smith");
+ }
+ fprintf(stdout, " normal form\n");
+ if(left_inverted == INVERTED)
+ {
+ fprintf(stdout, "The left factor will be inverted.\n");
+ }
+ if(right_inverted == INVERTED)
+ {
+ fprintf(stdout, "The right factor will be inverted.\n");
+ }
+ fprintf(stdout, "Matrix read from file %s\n", infile_name);
+ }
+
+ /*************************************************************/
+ /* We have all we need to get to the computational business: */
+ /*************************************************************/
+
+ /* Read the matrix */
+ A = gmp_matrix_read_coord(infile_name);
+ free(infile_name);
+
+ if(verbose_wanted == WANTED)
+ {
+ fprintf(stdout, "Input matrix:\n");
+ gmp_matrix_fprintf(stdout, A);
+ }
+
+ /* Compute the requested normal form */
+ switch(normal_form_requested)
+ {
+ case SMITH:
+ normal_form = create_gmp_Smith_normal_form(A, left_inverted, right_inverted);
+ break;
+
+ default:
+ normal_form = create_gmp_Hermite_normal_form(A, left_inverted, right_inverted);
+ break;
+ }
+
+
+ if(verbose_wanted == WANTED)
+ {
+ fprintf(stdout, "Left factor:\n");
+ gmp_matrix_fprintf(stdout, normal_form->left);
+ fprintf(stdout, "Canonical form:\n");
+ gmp_matrix_fprintf(stdout, normal_form->canonical);
+ fprintf(stdout, "Right factor:\n");
+ gmp_matrix_fprintf(stdout, normal_form->right);
+ }
+
+ /* Save the result to the output files */
+ outfile_name = append_suffix(file_base_name, "left");
+ if( gmp_matrix_write_coord(outfile_name, normal_form -> left)
+ == EXIT_SUCCESS)
+ {
+ if(verbose_wanted == WANTED)
+ fprintf(stdout, "Left factor successfully writen to file %s\n", outfile_name);
+ }
+
+ free(outfile_name);
+
+ outfile_name = append_suffix(file_base_name, "can");
+ if( gmp_matrix_write_coord(outfile_name, normal_form -> canonical)
+ == EXIT_SUCCESS)
+ {
+ if(verbose_wanted == WANTED)
+ fprintf(stdout, "Canonical form successfully writen to file %s\n", outfile_name);
+ }
+ free(outfile_name);
+
+ outfile_name = append_suffix(file_base_name, "right");
+ if( gmp_matrix_write_coord(outfile_name, normal_form -> right)
+ == EXIT_SUCCESS)
+ {
+ if(verbose_wanted == WANTED)
+ fprintf(stdout, "Right factor successfully writen to file %s\n", outfile_name);
+ }
+ free(outfile_name);
+
+ free(file_base_name);
+ destroy_gmp_normal_form(normal_form);
+
+ return EXIT_SUCCESS;
+}
diff --git a/contrib/kbipack/gmp_blas.c b/contrib/kbipack/gmp_blas.c
new file mode 100755
index 0000000000000000000000000000000000000000..501c4eba543f02980b4138c039bbbb8e920913db
--- /dev/null
+++ b/contrib/kbipack/gmp_blas.c
@@ -0,0 +1,227 @@
+/*
+ Source file for integer-oriented elementary versions of some
+ subroutines of BLAS. There routines are chosen to facilitate the
+ computation of Hermite and Smith normal forms of matrices of
+ modest size.
+
+ Copyright (C) 28.10.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: gmp_blas.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+/* #include<stdio.h> */
+#include"gmp_blas.h"
+
+void gmp_blas_swap(size_t n, mpz_t* x, size_t incx, mpz_t* y, size_t incy)
+{
+ mpz_t elbow;
+ size_t ind;
+
+ mpz_init(elbow);
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_set(elbow , x[ind*incx]);
+ mpz_set(x[ind*incx], y[ind*incy]);
+ mpz_set(y[ind*incy], elbow);
+ }
+ mpz_clear(elbow);
+
+ return;
+}
+
+void gmp_blas_scal(size_t n, mpz_t a, mpz_t* x, size_t incx)
+{
+ size_t ind;
+
+ for(ind = 0; ind < n; ind++)
+ {
+ mpz_mul (x[ind*incx], x[ind*incx], a);
+ }
+
+ return;
+}
+
+void gmp_blas_copy(size_t n, 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;
+}
diff --git a/contrib/kbipack/gmp_blas.h b/contrib/kbipack/gmp_blas.h
new file mode 100755
index 0000000000000000000000000000000000000000..6c9c88b2f8454539f9ba694239d52232b33e34df
--- /dev/null
+++ b/contrib/kbipack/gmp_blas.h
@@ -0,0 +1,103 @@
+/*
+ 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
diff --git a/contrib/kbipack/gmp_blas.o b/contrib/kbipack/gmp_blas.o
new file mode 100755
index 0000000000000000000000000000000000000000..02064581725915c72fb204793a3f9436e14c5cfc
Binary files /dev/null and b/contrib/kbipack/gmp_blas.o differ
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
new file mode 100755
index 0000000000000000000000000000000000000000..23d648c89eff75d10167c303d32132c6a5c9401b
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix.c
@@ -0,0 +1,673 @@
+/*
+ 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;
+}
diff --git a/contrib/kbipack/gmp_matrix.c~ b/contrib/kbipack/gmp_matrix.c~
new file mode 100755
index 0000000000000000000000000000000000000000..2afe8d1dd2d6b34312a4859e01d672a1480bf947
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix.c~
@@ -0,0 +1,678 @@
+/*
+ Source file for integer-oriented matrix, relying on the arbitrary
+ precision integers from the GNU gmp package.
+
+ Copyright (C) 28.10.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: gmp_matrix.c~,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+
+#include<stdlib.h>
+#include"gmp_matrix.h"
+
+int gmp_kaka(int laalaa){
+ printf("kissa %d", laalaa);
+ return laalaa;
+}
+
+gmp_matrix *
+create_gmp_matrix(size_t r, size_t c,
+ const mpz_t * e)
+{
+ gmp_matrix * new_matrix;
+ size_t ind;
+
+ new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix));
+ if(new_matrix == NULL)
+ {
+ return NULL;
+ }
+
+ new_matrix -> storage = (mpz_t *) calloc(r*c, sizeof(mpz_t));
+ if(new_matrix -> storage == NULL)
+ {
+ free(new_matrix);
+ return NULL;
+ }
+
+ new_matrix -> rows = r;
+ new_matrix -> cols = c;
+
+ for(ind = 0; ind < r*c; ind ++)
+ {
+ mpz_init (new_matrix -> storage[ind]);
+ mpz_set (new_matrix -> storage[ind], e[ind]);
+ }
+
+ return new_matrix;
+}
+
+
+gmp_matrix *
+create_gmp_matrix_identity(size_t dim)
+{
+ gmp_matrix * new_matrix;
+ size_t ind;
+
+ new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix));
+ if(new_matrix == NULL)
+ {
+ return NULL;
+ }
+
+ new_matrix -> storage = (mpz_t *) calloc(dim*dim, sizeof(mpz_t));
+ if(new_matrix -> storage == NULL)
+ {
+ free(new_matrix);
+ return NULL;
+ }
+
+ new_matrix -> rows = dim;
+ new_matrix -> cols = dim;
+
+ for(ind = 0; ind < dim*dim; ind ++)
+ {
+ mpz_init_set_si (new_matrix -> storage[ind], 0);
+ }
+
+ for(ind = 0; ind < dim; ind ++)
+ {
+ mpz_set_ui(new_matrix -> storage[ind*(dim+1)], 1);
+ }
+
+ return new_matrix;
+}
+
+
+gmp_matrix *
+create_gmp_matrix_zero(size_t rows, size_t cols)
+{
+ gmp_matrix * new_matrix;
+ size_t ind;
+
+ new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix));
+ if(new_matrix == NULL)
+ {
+ return NULL;
+ }
+
+ new_matrix -> storage = (mpz_t *) calloc(rows*cols, sizeof(mpz_t));
+ if(new_matrix -> storage == NULL)
+ {
+ free(new_matrix);
+ return NULL;
+ }
+
+ new_matrix -> rows = rows;
+ new_matrix -> cols = cols;
+
+ for(ind = 0; ind < rows*cols; ind ++)
+ {
+ mpz_init_set_si (new_matrix -> storage[ind], 0);
+ }
+
+ return new_matrix;
+}
+
+int
+destroy_gmp_matrix(gmp_matrix * m)
+{
+ size_t ind, nmb_storage;;
+
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ if(m -> storage == NULL)
+ {
+ free(m);
+ return EXIT_FAILURE;
+ }
+
+ nmb_storage = (m -> rows)*(m -> cols);
+
+ for(ind = 0; ind < nmb_storage; ind++)
+ {
+ mpz_clear(m -> storage[ind]);
+ }
+
+ free(m -> storage);
+ free(m);
+ return EXIT_SUCCESS;
+}
+
+
+size_t
+gmp_matrix_rows(const gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return 0;
+ }
+
+ return m -> rows;
+}
+
+
+size_t
+gmp_matrix_cols(const gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return 0;
+ }
+
+ return m -> cols;
+}
+
+/* elem <- (matrix(row, col)) */
+int
+gmp_matrix_get_elem(mpz_t elem, size_t row, size_t col,
+ const gmp_matrix * m)
+{
+#ifdef DEBUG
+
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+#endif
+
+ mpz_set(elem, m -> storage[(col-1)*(m -> rows)+row-1]);
+ return EXIT_SUCCESS;
+}
+
+int
+gmp_matrix_swap_rows(size_t row1, size_t row2, gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((row1 < 1) || (row1 > m->rows) || (row2 < 1) || (row2 > m->rows))
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* printf("Swapping rows %i %i\n", row1, row2); */
+ gmp_blas_swap(m -> cols,
+ &(m -> storage[row1-1]), m -> rows,
+ &(m -> storage[row2-1]), m -> rows);
+
+ return EXIT_SUCCESS;
+}
+
+int
+gmp_matrix_swap_cols(size_t col1, size_t col2, gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((col1 < 1) || (col1 > m->cols) || (col2 < 1) || (col2 > m->cols))
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* printf("Swapping cols %i %i\n", col1, col2); */
+ gmp_blas_swap(m -> rows,
+ &(m -> storage[(m->rows)*(col1-1)]), 1,
+ &(m -> storage[(m->rows)*(col2-1)]), 1);
+
+ return EXIT_SUCCESS;
+}
+
+int
+gmp_matrix_negate_row(size_t row, gmp_matrix * m)
+{
+ mpz_t minus_one;
+
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((row < 1) || (row > m->rows))
+ {
+ return EXIT_FAILURE;
+ }
+
+ mpz_init(minus_one);
+ mpz_set_si(minus_one, -1);
+ gmp_blas_scal(m -> cols, minus_one, (&m -> storage[row-1]), m -> rows);
+ mpz_clear(minus_one);
+ return EXIT_SUCCESS;
+}
+
+int
+gmp_matrix_negate_col(size_t col, gmp_matrix * m)
+{
+ mpz_t minus_one;
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((col < 1) || (col > m->cols))
+ {
+ return EXIT_FAILURE;
+ }
+
+ mpz_init(minus_one);
+ mpz_set_si(minus_one, -1);
+ gmp_blas_scal(m -> rows, minus_one,
+ &(m -> storage[(m->rows)*(col-1)]), 1);
+ mpz_clear(minus_one);
+ return EXIT_SUCCESS;
+}
+
+/* row2 <- a*row1 + row2*/
+int
+gmp_matrix_add_row(mpz_t a, size_t row1, size_t row2,
+ gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((row1 < 1) || (row1 > m->rows) || (row2 < 1) || (row2 > m->rows))
+ {
+ return EXIT_FAILURE;
+ }
+
+ gmp_blas_axpy(m->cols, a,
+ (const mpz_t *) &(m->storage[row1-1]), m->rows,
+ &(m->storage[row2-1]), m->rows);
+
+ return EXIT_SUCCESS;
+}
+
+int
+gmp_matrix_add_col(mpz_t a, size_t col1, size_t col2,
+ gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((col1 < 1) || (col1 > m->cols) || (col2 < 1) || (col2 > m->cols))
+ {
+ return EXIT_FAILURE;
+ }
+
+ gmp_blas_axpy(m->rows, a,
+ (const mpz_t *) &(m -> storage[(m->rows)*(col1-1)]), 1,
+ &(m -> storage[(m->rows)*(col2-1)]), 1);
+
+ return EXIT_SUCCESS;
+}
+
+/* row1 <- a*row1 + b*row2
+ row2 <- c*row1 + d*row2 */
+int
+gmp_matrix_row_rot(mpz_t a, mpz_t b, size_t row1,
+ mpz_t c, mpz_t d, size_t row2,
+ gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((row1 < 1) || (row1 > m->rows) || (row2 < 1) || (row2 > m->rows))
+ {
+ return EXIT_FAILURE;
+ }
+
+ gmp_blas_rot(m->cols,
+ a, b, &(m->storage[row1-1]), m->rows,
+ c, d, &(m->storage[row2-1]), m->rows);
+
+ return EXIT_SUCCESS;
+}
+int
+gmp_matrix_col_rot(mpz_t a, mpz_t b, size_t col1,
+ mpz_t c, mpz_t d, size_t col2,
+ gmp_matrix * m)
+{
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+ if((col1 < 1) || (col1 > m->cols) || (col2 < 1) || (col2 > m->cols))
+ {
+ return EXIT_FAILURE;
+ }
+
+/* printf("a: %i b: %i c: %i d:%i col1: %i col2: %i\n", */
+/* mpz_get_si(a), */
+/* mpz_get_si(b), */
+/* mpz_get_si(c), */
+/* mpz_get_si(d), */
+/* col1, col2); */
+
+ gmp_blas_rot(m->rows,
+ a, b, &(m -> storage[(m->rows)*(col1-1)]), 1,
+ c, d, &(m -> storage[(m->rows)*(col2-1)]), 1);
+
+ return EXIT_SUCCESS;
+}
+
+size_t
+gmp_matrix_col_inz (size_t r1, size_t r2, size_t c,
+ gmp_matrix * m)
+{
+ size_t result;
+
+ if(m == NULL)
+ {
+ return 0;
+ }
+ if((r1 < 1) || (r1 > m->rows) ||
+ (r2 < 1) || (r2 > m->rows) ||
+ (r2 < r1) || (c < 1) || (c > m->cols))
+ {
+ return 0;
+ }
+
+ result = gmp_blas_inz(r2-r1+1,
+ (const mpz_t *) &(m->storage[(c-1)*(m->rows)+r1-1]),
+ 1);
+
+ if(result > r2-r1+1)
+ {
+ return 0;
+ }
+
+ return result;
+}
+
+size_t
+gmp_matrix_row_inz (size_t r, size_t c1, size_t c2,
+ gmp_matrix * m)
+{
+ size_t result;
+
+ if(m == NULL)
+ {
+ return 0;
+ }
+ if((r < 1) || (r > m->rows) ||
+ (c1 < 1) || (c1 > m->cols) ||
+ (c2 < c1) || (c2 < 1) || (c2 > m->cols))
+ {
+ return 0;
+ }
+
+ result = gmp_blas_inz(c2-c1+1,
+ (const mpz_t *) &(m->storage[(c1-1)*(m->rows)+r-1]),
+ m->rows);
+
+ if(result > c2-c1+1)
+ {
+ return 0;
+ }
+
+ return result;
+}
+
+
+int
+gmp_matrix_is_diagonal(const gmp_matrix * M)
+{
+ size_t i,j;
+ size_t rows, cols;
+
+ if(M == NULL)
+ {
+ return 0;
+ }
+
+ rows = M->rows;
+ cols = M->cols;
+
+ for(j = 1; j <= cols; j ++)
+ {
+ for(i = 1; i <= rows; i ++)
+ {
+ if((mpz_cmp_si(M->storage[(i-1)+(j-1)*rows], 0) != 0) &&
+ (i != j))
+ {
+ return 0;
+ }
+ }
+ }
+
+ return 1;
+}
+
+
+int
+gmp_matrix_transp(gmp_matrix * M)
+{
+ mpz_t * new_storage;
+ size_t i,j;
+ size_t rows, cols;
+
+ if(M == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ rows = M->rows;
+ cols = M->cols;
+
+ new_storage = (mpz_t *) calloc(rows*cols, sizeof(mpz_t));
+ if(new_storage == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ for(i = 1; i <= rows; i++)
+ {
+ for(j = 1; j <= cols; j++)
+ {
+ mpz_init_set(new_storage[(j-1)+(i-1)*cols],
+ M-> storage[(i-1)+(j-1)*rows]);
+ mpz_clear(M-> storage[(i-1)+(j-1)*rows]);
+ }
+ }
+ free(M->storage);
+
+ M -> storage = new_storage;
+ M -> rows = cols;
+ M -> cols = rows;
+
+ return EXIT_SUCCESS;
+}
+
+
+int
+gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B)
+{
+ mpz_t * new_storage;
+ size_t i,j;
+ size_t rows_A, cols_A, rows_B, cols_B;
+
+ if((A == NULL) || (B == NULL))
+ {
+ return EXIT_FAILURE;
+ }
+
+ rows_A = A->rows;
+ cols_A = A->cols;
+ rows_B = B->rows;
+ cols_B = B->cols;
+
+ if(cols_A != rows_B)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* Create new storage for the product */
+ new_storage = (mpz_t *) calloc(rows_A*cols_B, sizeof(mpz_t));
+ if(new_storage == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* Compute the product to the storage */
+ for(j = 1; j <= cols_B; j++)
+ {
+ for(i = 1; i <= rows_A; i++)
+ {
+ mpz_init (new_storage[(i-1)+(j-1)*rows_A]);
+ gmp_blas_dot(&(new_storage[(i-1)+(j-1)*rows_A]),
+ cols_A,
+ (const mpz_t *) &(A->storage[i-1]), rows_A,
+ (const mpz_t *) &(B->storage[(j-1)*rows_B]), 1);
+ }
+ }
+
+ /* Get rid of the old storage */
+ for(i = 1; i <= rows_A*cols_A; i++)
+ {
+ mpz_clear (A->storage[i-1]);
+ }
+ free(A->storage);
+
+ /* Update A */
+ A -> storage = new_storage;
+ A -> cols = cols_B;
+
+ return EXIT_SUCCESS;
+}
+
+int
+gmp_matrix_left_mult(const gmp_matrix * A, gmp_matrix * B)
+{
+ mpz_t * new_storage;
+ size_t i,j;
+ size_t rows_A, cols_A, rows_B, cols_B;
+
+ if((A == NULL) || (B == NULL))
+ {
+ return EXIT_FAILURE;
+ }
+
+ rows_A = A->rows;
+ cols_A = A->cols;
+ rows_B = B->rows;
+ cols_B = B->cols;
+
+ if(cols_A != rows_B)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* Create new storage for the product */
+ new_storage = (mpz_t *) calloc(rows_A*cols_B, sizeof(mpz_t));
+ if(new_storage == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* Compute the product to the storage */
+ for(j = 1; j <= cols_B; j++)
+ {
+ for(i = 1; i <= rows_A; i++)
+ {
+ mpz_init (new_storage[(i-1)+(j-1)*rows_A]);
+ gmp_blas_dot(&(new_storage[(i-1)+(j-1)*rows_A]),
+ cols_A,
+ (const mpz_t *) &(A->storage[i-1]), rows_A,
+ (const mpz_t *) &(B->storage[(j-1)*rows_B]), 1);
+ }
+ }
+
+ /* Get rid of the old storage */
+ for(i = 1; i <= rows_B*cols_B; i++)
+ {
+ mpz_clear (B->storage[i-1]);
+ }
+ free(B->storage);
+
+ /* Update A */
+ B -> storage = new_storage;
+ B -> rows = rows_A;
+
+ return EXIT_SUCCESS;
+}
+
+
+int gmp_matrix_printf(const gmp_matrix * m)
+{
+
+ mpz_t outz;
+ size_t rows, cols, i, j;
+
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ rows = m->rows;
+ cols = m->cols;
+ mpz_init(outz);
+ for(i = 1; i <= rows ; i++)
+ {
+ for(j = 1; j <= cols ; j++)
+ {
+ gmp_matrix_get_elem(outz, i, j, m);
+ mpz_out_str (stdout, 10, outz);
+ printf(" ");
+ }
+ printf("\n");
+ }
+
+ mpz_clear(outz);
+
+ return EXIT_SUCCESS;
+}
+
+int gmp_matrix_fprintf(FILE* stream, const gmp_matrix * m)
+{
+ mpz_t outz;
+ size_t rows, cols, i, j;
+
+ if(m == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ rows = m->rows;
+ cols = m->cols;
+ mpz_init(outz);
+ for(i = 1; i <= rows ; i++)
+ {
+ for(j = 1; j <= cols ; j++)
+ {
+ gmp_matrix_get_elem(outz, i, j, m);
+ mpz_out_str (stream, 10, outz);
+ printf(" ");
+ }
+ printf("\n");
+ }
+
+ mpz_clear(outz);
+
+ return EXIT_SUCCESS;
+}
diff --git a/contrib/kbipack/gmp_matrix.h b/contrib/kbipack/gmp_matrix.h
new file mode 100755
index 0000000000000000000000000000000000000000..6e32583db36d840637cc64d42aac9236c449e3c4
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix.h
@@ -0,0 +1,142 @@
+/*
+ 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
diff --git a/contrib/kbipack/gmp_matrix.h~ b/contrib/kbipack/gmp_matrix.h~
new file mode 100755
index 0000000000000000000000000000000000000000..217b5fe7e1d1279cf617579cdc2abefaedfc9c31
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix.h~
@@ -0,0 +1,144 @@
+/*
+ Header file for integer-oriented matrix, relying on the arbitrary
+ precision integers from the GNU gmp package.
+
+ Copyright (C) 28.10.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: gmp_matrix.h~,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+#ifndef __GMP_MATRIX_H__
+#define __GMP_MATRIX_H__
+
+#include"gmp_blas.h"
+
+typedef struct
+{
+ size_t rows;
+ size_t cols;
+ mpz_t * storage;
+} gmp_matrix;
+
+int gmp_kaka(int laalaa);
+
+/* Sets the values of "elems" column by column. The user is
+ responsible for sufficient supply. */
+gmp_matrix *
+create_gmp_matrix(size_t rows, size_t cols,
+ const mpz_t * elems);
+
+gmp_matrix *
+create_gmp_matrix_identity(size_t dim);
+
+gmp_matrix *
+create_gmp_matrix_zero(size_t rows, size_t cols);
+
+int
+destroy_gmp_matrix(gmp_matrix *);
+
+size_t
+gmp_matrix_rows(const gmp_matrix *);
+
+size_t
+gmp_matrix_cols(const gmp_matrix *);
+
+/* elem <- (matrix(row, col)) */
+int
+gmp_matrix_get_elem(mpz_t elem, size_t row, size_t col,
+ const gmp_matrix *);
+
+int
+gmp_matrix_swap_rows(size_t row1, size_t row2, gmp_matrix *);
+
+int
+gmp_matrix_swap_cols(size_t col1, size_t col2, gmp_matrix *);
+
+int
+gmp_matrix_negate_row(size_t row, gmp_matrix *);
+
+int
+gmp_matrix_negate_col(size_t col, gmp_matrix *);
+
+/* row2 <- a*row1 + row2*/
+int
+gmp_matrix_add_row(mpz_t a, size_t row1, size_t row2,
+ gmp_matrix *);
+int
+gmp_matrix_add_col(mpz_t a, size_t col1, size_t col2,
+ gmp_matrix *);
+
+/* row1 <- a*row1 + b*row2
+ row2 <- c*row1 + d*row2 */
+int
+gmp_matrix_row_rot(mpz_t a, mpz_t b, size_t row1,
+ mpz_t c, mpz_t d, size_t row2,
+ gmp_matrix *);
+int
+gmp_matrix_col_rot(mpz_t a, mpz_t b, size_t col1,
+ mpz_t c, mpz_t d, size_t col2,
+ gmp_matrix *);
+
+/* 0 for no, 1 for yes */
+int
+gmp_matrix_is_diagonal(const gmp_matrix * M);
+
+/* Finds a nonzero in a subcolumn M(r1:r2,c). */
+/* Returns zero if no nonzeros found. */
+size_t
+gmp_matrix_col_inz (size_t r1, size_t r2, size_t c,
+ gmp_matrix * M);
+
+/* Finds a nonzero in a subrow M(r,c1:c2). */
+/* Returns zero if no nonzeros found. */
+size_t
+gmp_matrix_row_inz (size_t r, size_t c1, size_t c2,
+ gmp_matrix * M);
+
+int
+gmp_matrix_transp(gmp_matrix * M);
+
+/* A <- A*B */
+int
+gmp_matrix_right_mult(gmp_matrix * A, const gmp_matrix * B);
+
+/* B <- A*B */
+int
+gmp_matrix_left_mult(const gmp_matrix * A, gmp_matrix * B);
+
+/* (TBD ?)Implement the BLAS style GEMM? Place it here at all? */
+/* (T) (T) */
+/* A <- B * C */
+/* Give 1 if the transpose of the matrix is to be used, else 0 */
+/* int */
+/* gmp_matrix_mult(gmp_matrix * A, */
+/* const gmp_matrix * B, int transpose_B, */
+/* const gmp_matrix * C, int transpose_C); */
+
+/*
+ Mainly for diagnostics ...
+ --------------------------
+*/
+
+int gmp_matrix_printf(const gmp_matrix *);
+int gmp_matrix_fprintf(FILE*, const gmp_matrix *);
+
+
+#endif
diff --git a/contrib/kbipack/gmp_matrix.o b/contrib/kbipack/gmp_matrix.o
new file mode 100755
index 0000000000000000000000000000000000000000..b151210c86fa837e700afce35afff9326bbec789
Binary files /dev/null and b/contrib/kbipack/gmp_matrix.o differ
diff --git a/contrib/kbipack/gmp_matrix_io.c b/contrib/kbipack/gmp_matrix_io.c
new file mode 100755
index 0000000000000000000000000000000000000000..f85639304a34e500096a2d3196a23044c8a212e8
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix_io.c
@@ -0,0 +1,135 @@
+/*
+ IO-routines for gmp integer matrices in coordinate format.
+
+ Copyright (C) 19.11.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: gmp_matrix_io.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+ */
+
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+#include "gmp_matrix_io.h"
+
+gmp_matrix * gmp_matrix_read_coord(char* filename)
+{
+ FILE * p_file;
+ size_t read_values;
+ size_t row, col, nrows, ncols, dummy;
+ int val;
+ char buffer[1000];
+ gmp_matrix * new_matrix;
+
+ p_file = fopen(filename, "r");
+ if(p_file == NULL)
+ return NULL;
+
+ /* Matlab and Octave typically include comments on ascii files. */
+ /* They start with # */
+ fgets(buffer, 999, p_file);
+ while(strncmp(buffer, "#",1) == 0)
+ {
+ fgets(buffer, 999, p_file);
+ }
+
+ /* First read the size of the matrix */
+ read_values = sscanf(buffer, "%u %u %u", &nrows, &ncols, &dummy);
+
+ /* Create the matrix */
+ new_matrix = create_gmp_matrix_zero(nrows, ncols);
+ if(new_matrix == NULL)
+ {
+ fclose(p_file);
+ return NULL;
+ }
+
+ /* Read the values to the matrix */
+ while(read_values != EOF)
+ {
+ read_values = fscanf(p_file, "%u %u %i\n", &row, &col, &val);
+ if( (row <= nrows) && (row > 0) && (col <= ncols) && (col > 0) )
+ {
+ mpz_set_si ( new_matrix->storage[(row-1)+(col-1)*nrows], val );
+ }
+ }
+ fclose(p_file);
+
+ return new_matrix;
+}
+
+int gmp_matrix_write_coord(char* filename, 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;
+}
diff --git a/contrib/kbipack/gmp_matrix_io.h b/contrib/kbipack/gmp_matrix_io.h
new file mode 100755
index 0000000000000000000000000000000000000000..74df3f2ad138c63068ac277064a72875d93e0c59
--- /dev/null
+++ b/contrib/kbipack/gmp_matrix_io.h
@@ -0,0 +1,37 @@
+/*
+ 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
diff --git a/contrib/kbipack/gmp_matrix_io.o b/contrib/kbipack/gmp_matrix_io.o
new file mode 100755
index 0000000000000000000000000000000000000000..c7231192ab4007c70d68c083fb81cab8a3fddd0f
Binary files /dev/null and b/contrib/kbipack/gmp_matrix_io.o differ
diff --git a/contrib/kbipack/gmp_normal_form.c b/contrib/kbipack/gmp_normal_form.c
new file mode 100755
index 0000000000000000000000000000000000000000..42785337e7e54cf85e1bc82b1dc0ea52b14ed925
--- /dev/null
+++ b/contrib/kbipack/gmp_normal_form.c
@@ -0,0 +1,783 @@
+/*
+ Implementation for integer computation of Hermite and Smith normal
+ forms of matrices of modest size.
+
+ Implementation: Dense matrix with GMP-library's mpz_t elements to
+ hold huge integers.
+
+ Algorithm: Kannan - Bachem algorithm with improvement by
+ Chou and Collins. Expects a large number of unit invariant
+ factors and takes advantage of them as they appear.
+
+ References:
+ [1] Ravindran Kannan, Achim Bachem:
+ "Polynomial algorithms for computing the Smith and Hermite normal
+ forms of an integer matrix",
+ SIAM J. Comput., vol. 8, no. 5, pp. 499-507, 1979.
+ [2] Tsu-Wu J.Chou, George E. Collins:
+ "Algorithms for the solution of systems of linear Diophantine
+ equations",
+ SIAM J. Comput., vol. 11, no. 4, pp. 687-708, 1982.
+ [3] GMP homepage http://www.swox.com/gmp/
+ [4] GNU gmp page http://www.gnu.org/software/gmp/
+
+ Copyright (C) 30.10.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: gmp_normal_form.c,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+#include<stdlib.h>
+#include"gmp_blas.h"
+#include"gmp_matrix.h"
+#include"gmp_normal_form.h"
+
+
+/* The unaltered matrix and identity left and right factors */
+static gmp_normal_form *
+create_gmp_trivial_normal_form(gmp_matrix * A,
+ inverted_flag left_inverted,
+ inverted_flag right_inverted)
+{
+ size_t rows, cols;
+
+ gmp_normal_form * new_nf;
+ if(A == NULL)
+ {
+ return NULL;
+ }
+
+ new_nf = (gmp_normal_form *) malloc(sizeof(gmp_normal_form));
+ if(new_nf == NULL)
+ {
+ destroy_gmp_matrix(A);
+ return NULL;
+ }
+
+ rows = A -> rows;
+ cols = A -> cols;
+
+ if((rows == 0) || (cols == 0))
+ {
+ destroy_gmp_matrix(A);
+ return NULL;
+ }
+
+ new_nf->left = create_gmp_matrix_identity(rows);
+ if(new_nf->left == NULL)
+ {
+ destroy_gmp_matrix(A);
+ free(new_nf);
+ return NULL;
+ }
+
+ new_nf->right = create_gmp_matrix_identity(cols);
+ if(new_nf->right == NULL)
+ {
+ destroy_gmp_matrix(A);
+ destroy_gmp_matrix(new_nf->left);
+ free(new_nf);
+ return NULL;
+ }
+
+ new_nf -> canonical = A;
+ new_nf -> left_inverted = left_inverted;
+ new_nf -> right_inverted = right_inverted;
+
+ return new_nf;
+}
+
+static int
+gmp_Hermite_eliminate_step(gmp_matrix * L, gmp_matrix * U,
+ size_t col, inverted_flag right_inverted)
+{
+ size_t ind, row_limit;
+ mpz_t pivot, elem;
+ mpz_t bez1, bez2, gc_div;
+ mpz_t cff1, cff2;
+
+ mpz_init(pivot);
+ mpz_init(elem);
+ mpz_init(bez1);
+ mpz_init(bez2);
+ mpz_init(cff1);
+ mpz_init(cff2);
+ mpz_init(gc_div);
+
+ row_limit = (L->rows >= col) ?
+ col-1 :
+ L->rows;
+
+ for(ind = 1; ind <= row_limit; ind++)
+ {
+ gmp_matrix_get_elem(elem, ind, col, L);
+
+ /* Eliminate only if nonzero */
+ if(mpz_sgn (elem) != 0)
+ {
+ gmp_matrix_get_elem(pivot, ind, ind, L);
+
+ /* Extended Euclid's:
+ bez1*pivot+bez2*elem = gc_div */
+ gmp_blas_rotg(gc_div, bez1, bez2, pivot, elem);
+
+ /* Make cff1 = -elem/gc_div */
+ mpz_divexact(cff1, elem, gc_div);
+ mpz_neg (cff1, cff1);
+ /* Make cff2 = pivot/gc_div */
+ mpz_divexact(cff2, pivot, gc_div);
+
+ /* Update the HNF canonical matrix */
+ gmp_matrix_col_rot(bez1, bez2, ind,
+ cff1, cff2, col,
+ L);
+
+ /* Update the unimodular factor matrix */
+ if(right_inverted == INVERTED)
+ {
+ gmp_matrix_col_rot(bez1, bez2, ind,
+ cff1, cff2, col,
+ U);
+ }
+ else
+ {
+
+ /* [bez1, -elem/gc_div] [pivot/gc_div, elem/gc_div]
+ [bez2, pivot/gc_div] [-bez2, bez1 ] = I_2 */
+ mpz_neg(cff1, cff1);
+ mpz_neg(bez2, bez2);
+ gmp_matrix_row_rot(cff2, cff1, ind,
+ bez2, bez1, col,
+ U);
+ }
+ }
+ }
+ mpz_clear(pivot);
+ mpz_clear(elem);
+ mpz_clear(bez1);
+ mpz_clear(bez2);
+ mpz_clear(cff1);
+ mpz_clear(cff2);
+ mpz_clear(gc_div);
+
+ return EXIT_SUCCESS;
+}
+
+
+
+static int
+gmp_Hermite_reduce_step(gmp_matrix * L, gmp_matrix * U,
+ size_t col, inverted_flag right_inverted)
+{
+
+ size_t i, j;
+ mpz_t pivot, elem;
+ mpz_t cff;
+
+ mpz_init(pivot);
+ mpz_init(elem);
+ mpz_init(cff);
+
+ if(col > (L->rows))
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* printf("Col = %i\n", col);
+ printf("L before\n");
+ gmp_matrix_printf(L);*/
+
+ for(j = col-1; j >= 1; j--)
+ {
+ for(i = j+1; i <= col; i++)
+ {
+ gmp_matrix_get_elem(pivot, i, i, L);
+ gmp_matrix_get_elem(elem, i, j, L);
+ /* printf(" i %i j %i\n",i,j); */
+
+ if((mpz_cmpabs(elem, pivot) >= 0) || (mpz_sgn(elem) < 0))
+ {
+ /* The objective:
+ 0 <= elem + k*pivot < pivot
+ -elem <= k*pivot < pivot - elem
+ -elem/pivot <= k < - elem/pivot + 1
+
+ Solution:
+ k = ceil(-elem/pivot)
+ */
+
+ /* Make cff = -elem */
+ mpz_neg (cff, elem);
+ mpz_cdiv_q (cff, cff, pivot);
+
+ /* printf("col %i j %i\n", i, j);
+ printf("elem %i k %i pivot %i\n",
+ mpz_get_si(elem),
+ mpz_get_si(cff),
+ mpz_get_si(pivot));*/
+
+
+ gmp_matrix_add_col(cff, i, j, L);
+
+ /* Update the unimodular factor matrix */
+ if(right_inverted == INVERTED)
+ {
+ gmp_matrix_add_col(cff, i, j, U);
+ }
+
+ /* [1 0] [ 1 0] = [1 0]
+ [cff 1] [-cff 1] [0 1] */
+ else
+ {
+ mpz_neg (cff, cff);
+ gmp_matrix_add_row(cff, j, i, U);
+ }
+
+ /* printf("\n");
+ gmp_matrix_printf(L);*/
+
+ }
+ }
+ }
+
+/* printf("L after\n"); */
+/* gmp_matrix_printf(L); */
+
+ mpz_clear(pivot);
+ mpz_clear(elem);
+ mpz_clear(cff);
+ return EXIT_SUCCESS;
+}
+
+
+static int
+gmp_normal_form_make_Hermite(gmp_normal_form * nf)
+{
+ size_t rows, cols;
+ size_t pivot_ind;
+ size_t schur, colind;
+ mpz_t pivot;
+ gmp_matrix * canonical;
+ gmp_matrix * left;
+ gmp_matrix * right;
+
+ if(nf == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* OK, it's safe to get to business */
+ canonical = nf->canonical;
+ left = nf->left;
+ right = nf->right;
+ rows = canonical -> rows;
+ cols = canonical -> cols;
+
+ mpz_init(pivot);
+
+ /* "schur" denotes where the present Schur complement starts */
+ schur = 1;
+
+ while((schur <= rows) && (schur <= cols))
+ {
+ /* Eliminate a column */
+ if (schur > 1)
+ {
+ gmp_Hermite_eliminate_step(canonical, right,
+ schur, nf->right_inverted);
+ }
+
+ /* Find a pivot */
+ pivot_ind = gmp_matrix_col_inz(schur, rows, schur, canonical);
+
+
+ /* If no nonzeros was found, the column is all zero, hence
+ settled with. Permute it to the end and decrement cols. */
+ if(pivot_ind == 0)
+ {
+ gmp_matrix_swap_cols(schur, cols, canonical);
+
+ if(nf -> right_inverted == INVERTED)
+ {
+ gmp_matrix_swap_cols(schur, cols, right);
+ }
+ else
+ {
+ gmp_matrix_swap_rows(schur, cols, right);
+ }
+
+ cols--;
+
+ /* When the whole column was zeroed, the diagonal
+ elements may have got reduced. Reduce the sub-
+ diagonals as well*/
+
+ if(schur > 1)
+ {
+ gmp_Hermite_reduce_step (canonical, right, schur-1,
+ nf -> right_inverted);
+ }
+ }
+
+ /* A nonzero pivot was found. Permute it to the diagonal position,
+ make it positive, and reduce the off-diagonals.
+ The schur complement now starts from the next diagonal. */
+ else
+ {
+ pivot_ind = schur+pivot_ind-1;
+ gmp_matrix_swap_rows(schur, pivot_ind, canonical);
+
+ if(nf->left_inverted == INVERTED)
+ {
+ gmp_matrix_swap_rows(schur, pivot_ind, left);
+ }
+ else
+ {
+ gmp_matrix_swap_cols(schur, pivot_ind, left);
+ }
+
+ /* Make the pivot positive */
+ gmp_matrix_get_elem(pivot, schur, schur, canonical);
+
+ if(mpz_cmp_si(pivot, 0) < 0)
+ {
+ gmp_matrix_negate_col(schur, canonical);
+
+ if(nf->right_inverted == INVERTED)
+ {
+ gmp_matrix_negate_col(schur, right);
+ }
+ else
+ {
+ gmp_matrix_negate_row(schur, right);
+ }
+ }
+
+ /* Reduce off-diagonals */
+ gmp_Hermite_reduce_step (canonical, right, schur, nf -> right_inverted);
+
+ schur++;
+ }
+ }
+
+ /* The Schur complement is now empty. There may still be uneliminated
+ columns left (in case of a wide matrix) */
+
+ colind = schur;
+ while (colind <= cols)
+ {
+ gmp_Hermite_eliminate_step (canonical, right, colind, nf->right_inverted);
+ gmp_Hermite_reduce_step (canonical, right, schur-1, nf->right_inverted);
+ colind++;
+ }
+
+ mpz_clear(pivot);
+
+ return EXIT_SUCCESS;
+}
+
+
+
+gmp_normal_form *
+create_gmp_Hermite_normal_form(gmp_matrix * A,
+ inverted_flag left_inverted,
+ inverted_flag right_inverted)
+{
+ gmp_normal_form * new_nf;
+
+ if(A == NULL)
+ {
+ return NULL;
+ }
+
+ new_nf =
+ create_gmp_trivial_normal_form(A, left_inverted, right_inverted);
+
+ if(new_nf == NULL)
+ {
+ /* A has been destroyed already */
+ return NULL;
+ }
+
+ if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ return NULL;
+ }
+
+ return new_nf;
+}
+
+
+gmp_normal_form *
+create_gmp_Smith_normal_form(gmp_matrix * A,
+ inverted_flag left_inverted,
+ inverted_flag right_inverted)
+{
+
+ gmp_normal_form * new_nf;
+ gmp_matrix * canonical;
+ gmp_matrix * elbow;
+ inverted_flag elbow_flag;
+ size_t rows, cols;
+ size_t i, j;
+ size_t subd_ind, row_undivisible;
+ size_t last_ready_row, first_ready_col, ind;
+ mpz_t pivot;
+ mpz_t remainder;
+
+ if(A == NULL)
+ {
+ return NULL;
+ }
+
+ new_nf =
+ create_gmp_trivial_normal_form(A, left_inverted, right_inverted);
+
+ if(new_nf == NULL)
+ {
+ /* A has been destroyed already */
+ return NULL;
+ }
+
+ canonical = new_nf -> canonical;
+ mpz_init(pivot);
+ mpz_init(remainder);
+ rows = canonical->rows;
+ cols = canonical->cols;
+ last_ready_row = 0;
+ first_ready_col = (cols < rows) ? (cols+1) : (rows+1);
+
+ while(first_ready_col > last_ready_row + 1)
+ {
+ /*******/
+ /* HNF */
+ /*******/
+
+ if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+
+#ifdef DEBUG
+ printf("HNF\n");
+ gmp_matrix_printf (new_nf -> left);
+ gmp_matrix_printf (new_nf -> canonical);
+ gmp_matrix_printf (new_nf -> right);
+#endif
+
+ /**********************/
+ /* Find ready columns */
+ /**********************/
+
+ /* If a diagonal entry is zero, so is the corresponding
+ column. The zero diagonals always reside in the end.
+ Seek until zero diagonal encountered, but stay within the matrix! */
+ ind = 1;
+ while ( (ind < first_ready_col) &&
+ (mpz_cmp_si(canonical -> storage[(ind-1)+(ind-1)*rows], 0) != 0)
+ )
+ {
+ ind ++;
+ }
+ first_ready_col = ind;
+
+ /* Note: The number of ready cols is settled after the first HNF,
+ but the check is cheap. */
+
+ /**********************************************/
+ /* Permute unit diagonals such that they lead */
+ /**********************************************/
+
+ /* If the recently computed HNF has ones on the diagonal, their
+ corresponding rows are all zero (except the diagonal).
+ They are then settled, because the next LHNF kills the elements
+ on their columns. */
+
+ ind = last_ready_row+1;
+
+ /* Stay within the nonzero cols of the matrix */
+ while (ind < first_ready_col)
+ {
+ /* Unit diagonal encountered */
+ if(mpz_cmp_si ( canonical->storage[(ind-1) + (ind-1)*rows],
+ 1) == 0)
+ {
+ /* If not in the beginning, permute to extend the leading minor
+ with unit diagonals */
+ if(ind != last_ready_row+1)
+ {
+ gmp_matrix_swap_rows(last_ready_row+1, ind,
+ new_nf -> canonical);
+
+ if(left_inverted == INVERTED)
+ {
+ gmp_matrix_swap_rows(last_ready_row+1, ind,
+ new_nf -> left);
+ }
+ else
+ {
+ gmp_matrix_swap_cols(last_ready_row+1, ind,
+ new_nf -> left);
+ }
+
+ gmp_matrix_swap_cols(last_ready_row+1, ind,
+ new_nf -> canonical);
+
+ if(right_inverted == INVERTED)
+ {
+ gmp_matrix_swap_cols(last_ready_row+1, ind,
+ new_nf -> right);
+ }
+ else
+ {
+ gmp_matrix_swap_rows(last_ready_row+1, ind,
+ new_nf -> right);
+ }
+ }
+ last_ready_row ++;
+ }
+ ind++;
+ }
+
+#ifdef DEBUG
+ printf("Leading units\n");
+ gmp_matrix_printf (new_nf -> left);
+ gmp_matrix_printf (new_nf -> canonical);
+ gmp_matrix_printf (new_nf -> right);
+#endif
+
+ /********/
+ /* LHNF */
+ /********/
+
+ /* Extravagant strategy: Compute HNF of an explicit tranpose. */
+
+ /* 1. Transpose everything in the normal form */
+ if(gmp_matrix_transp(canonical) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->left) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->right) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+
+ elbow = new_nf ->left;
+ new_nf ->left = new_nf ->right;
+ new_nf ->right = elbow;
+ elbow_flag = new_nf ->left_inverted;
+ new_nf ->left_inverted = new_nf ->right_inverted;
+ new_nf ->right_inverted = elbow_flag;
+
+ /* 2. Compute HNF of the transpose */
+ if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+
+ /* 3. Transpose everything back */
+ if(gmp_matrix_transp(canonical) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->left) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->right) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ elbow = new_nf ->left;
+ new_nf ->left = new_nf ->right;
+ new_nf ->right = elbow;
+ elbow_flag = new_nf ->left_inverted;
+ new_nf ->left_inverted = new_nf ->right_inverted;
+ new_nf ->right_inverted = elbow_flag;
+
+
+#ifdef DEBUG
+ printf("LHNF\n");
+ gmp_matrix_printf (new_nf -> left);
+ gmp_matrix_printf (new_nf -> canonical);
+ gmp_matrix_printf (new_nf -> right);
+#endif
+
+
+ /*****************************************************/
+ /* Check if more of the leading normal form is ready */
+ /*****************************************************/
+
+ /* The matrix is in LHNF, i.e. it is upper triangular.
+ If the row trailing the top left diagonal element is
+ zero, the diagonal element may be an invariant factor
+ on its final position, and the stage may be ready.
+
+ The stage may not still be ready: The leading diagonal element
+ of D may not divide the rest of the Schur complement. */
+
+ subd_ind = 0;
+ row_undivisible = 0;
+
+ /* Explanation of loop conditions:
+ 1.) No relative primes found from Schur complement
+ 2.) Stay within the Schur complement
+ 3.) If a nonzero is found from the trailing row, the stage is
+ definitely not ready */
+ while (row_undivisible == 0 &&
+ last_ready_row + 1 < first_ready_col &&
+ subd_ind == 0)
+ {
+ subd_ind =
+ gmp_matrix_row_inz(last_ready_row+1,
+ last_ready_row+2, cols,
+ canonical);
+
+ /* printf("subd_ind %i\n", subd_ind);
+ printf("last_ready_row %i\n", last_ready_row); */
+
+ /* No nonzeros found, the stage may be ready */
+ if (subd_ind == 0)
+ {
+ mpz_set (pivot,
+ canonical->storage[(last_ready_row)+
+ (last_ready_row)*rows]);
+
+ /* Check whether the pivot divides all elements in the Schur
+ complement */
+ row_undivisible = 0;
+ for(j = last_ready_row+2; j < first_ready_col; j++)
+ {
+ for(i = last_ready_row+2; i <= j; i++)
+ {
+ mpz_tdiv_r (remainder,
+ canonical->storage[(i-1)+
+ (j-1)*rows],
+ pivot);
+
+ if(mpz_cmp_si(remainder, 0) !=0)
+ {
+ row_undivisible = i;
+ /* col_undivisible = j; */
+ }
+ }
+ }
+
+ /* printf("Row undivisible %i\n", row_undivisible); */
+
+ /* If a relative prime was found from the Schur complement,
+ add that row to the first row of the Schur complement */
+ if(row_undivisible != 0)
+ {
+ mpz_set_si (remainder, 1);
+ gmp_matrix_add_row(remainder,
+ row_undivisible, last_ready_row+1,
+ canonical);
+
+ /* [ 1 0] [1 0] = [1 0]
+ [-1 1] [1 1] [0 1] */
+
+ if(left_inverted == INVERTED)
+ {
+ gmp_matrix_add_row(remainder,
+ row_undivisible, last_ready_row+1,
+ new_nf->left);
+ }
+ else
+ {
+ mpz_neg (remainder, remainder);
+ gmp_matrix_add_col(remainder,
+ last_ready_row+1, row_undivisible,
+ new_nf->left);
+ }
+ }
+ /* The Schur complement is smaller by one again */
+ else
+ {
+ last_ready_row++;
+ }
+ }
+ }
+ } /* The main loop ends here */
+
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return new_nf;
+}
+
+int destroy_gmp_normal_form(gmp_normal_form * nf)
+{
+ int status;
+
+ if(nf == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ status = EXIT_SUCCESS;
+ if(destroy_gmp_matrix(nf -> left) == EXIT_FAILURE)
+ {
+ status = EXIT_FAILURE;
+ }
+
+ if(destroy_gmp_matrix(nf -> canonical) == EXIT_FAILURE)
+ {
+ status = EXIT_FAILURE;
+ }
+
+ if(destroy_gmp_matrix(nf -> right) == EXIT_FAILURE)
+ {
+ status = EXIT_FAILURE;
+ }
+ free(nf);
+
+ return status;
+}
+
+
diff --git a/contrib/kbipack/gmp_normal_form.c~ b/contrib/kbipack/gmp_normal_form.c~
new file mode 100755
index 0000000000000000000000000000000000000000..5f695a53d00542bd82baea3889dcebefac08aea4
--- /dev/null
+++ b/contrib/kbipack/gmp_normal_form.c~
@@ -0,0 +1,784 @@
+/*
+ Implementation for integer computation of Hermite and Smith normal
+ forms of matrices of modest size.
+
+ Implementation: Dense matrix with GMP-library's mpz_t elements to
+ hold huge integers.
+
+ Algorithm: Kannan - Bachem algorithm with improvement by
+ Chou and Collins. Expects a large number of unit invariant
+ factors and takes advantage of them as they appear.
+
+ References:
+ [1] Ravindran Kannan, Achim Bachem:
+ "Polynomial algorithms for computing the Smith and Hermite normal
+ forms of an integer matrix",
+ SIAM J. Comput., vol. 8, no. 5, pp. 499-507, 1979.
+ [2] Tsu-Wu J.Chou, George E. Collins:
+ "Algorithms for the solution of systems of linear Diophantine
+ equations",
+ SIAM J. Comput., vol. 11, no. 4, pp. 687-708, 1982.
+ [3] GMP homepage http://www.swox.com/gmp/
+ [4] GNU gmp page http://www.gnu.org/software/gmp/
+
+ Copyright (C) 30.10.2003 Saku Suuriniemi TUT/CEM
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+ Saku Suuriniemi, TUT/Electromagetics
+ P.O.Box 692, FIN-33101 Tampere, Finland
+ saku.suuriniemi@tut.fi
+
+ $Id: gmp_normal_form.c~,v 1.1 2009-03-30 14:10:57 matti Exp $
+*/
+
+#include<stdlib.h>
+#include"gmp_blas.h"
+#include"gmp_matrix.h"
+#include"gmp_normal_form.h"
+
+
+/* The unaltered matrix and identity left and right factors */
+static gmp_normal_form *
+create_gmp_trivial_normal_form(gmp_matrix * A,
+ inverted_flag left_inverted,
+ inverted_flag right_inverted)
+{
+ size_t rows, cols;
+
+ gmp_normal_form * new_nf;
+ if(A == NULL)
+ {
+ return NULL;
+ }
+
+ new_nf = (gmp_normal_form *) malloc(sizeof(gmp_normal_form));
+ if(new_nf == NULL)
+ {
+ destroy_gmp_matrix(A);
+ return NULL;
+ }
+
+ rows = A -> rows;
+ cols = A -> cols;
+
+ if((rows == 0) || (cols == 0))
+ {
+ destroy_gmp_matrix(A);
+ return NULL;
+ }
+
+ new_nf->left = create_gmp_matrix_identity(rows);
+ if(new_nf->left == NULL)
+ {
+ destroy_gmp_matrix(A);
+ free(new_nf);
+ return NULL;
+ }
+
+ new_nf->right = create_gmp_matrix_identity(cols);
+ if(new_nf->right == NULL)
+ {
+ destroy_gmp_matrix(A);
+ destroy_gmp_matrix(new_nf->left);
+ free(new_nf);
+ return NULL;
+ }
+
+ new_nf -> canonical = A;
+ new_nf -> left_inverted = left_inverted;
+ new_nf -> right_inverted = right_inverted;
+
+ return new_nf;
+}
+
+static int
+gmp_Hermite_eliminate_step(gmp_matrix * L, gmp_matrix * U,
+ size_t col, inverted_flag right_inverted)
+{
+ size_t ind, row_limit;
+ mpz_t pivot, elem;
+ mpz_t bez1, bez2, gc_div;
+ mpz_t cff1, cff2;
+
+ mpz_init(pivot);
+ mpz_init(elem);
+ mpz_init(bez1);
+ mpz_init(bez2);
+ mpz_init(cff1);
+ mpz_init(cff2);
+ mpz_init(gc_div);
+
+ row_limit = (L->rows >= col) ?
+ col-1 :
+ L->rows;
+
+ for(ind = 1; ind <= row_limit; ind++)
+ {
+ gmp_matrix_get_elem(elem, ind, col, L);
+
+ /* Eliminate only if nonzero */
+ if(mpz_sgn (elem) != 0)
+ {
+ gmp_matrix_get_elem(pivot, ind, ind, L);
+
+ /* Extended Euclid's:
+ bez1*pivot+bez2*elem = gc_div */
+ gmp_blas_rotg(gc_div, bez1, bez2, pivot, elem);
+
+ /* Make cff1 = -elem/gc_div */
+ mpz_divexact(cff1, elem, gc_div);
+ mpz_neg (cff1, cff1);
+ /* Make cff2 = pivot/gc_div */
+ mpz_divexact(cff2, pivot, gc_div);
+
+ /* Update the HNF canonical matrix */
+ gmp_matrix_col_rot(bez1, bez2, ind,
+ cff1, cff2, col,
+ L);
+
+ /* Update the unimodular factor matrix */
+ if(right_inverted == INVERTED)
+ {
+ gmp_matrix_col_rot(bez1, bez2, ind,
+ cff1, cff2, col,
+ U);
+ }
+ else
+ {
+
+ /* [bez1, -elem/gc_div] [pivot/gc_div, elem/gc_div]
+ [bez2, pivot/gc_div] [-bez2, bez1 ] = I_2 */
+ mpz_neg(cff1, cff1);
+ mpz_neg(bez2, bez2);
+ gmp_matrix_row_rot(cff2, cff1, ind,
+ bez2, bez1, col,
+ U);
+ }
+ }
+ }
+ mpz_clear(pivot);
+ mpz_clear(elem);
+ mpz_clear(bez1);
+ mpz_clear(bez2);
+ mpz_clear(cff1);
+ mpz_clear(cff2);
+ mpz_clear(gc_div);
+
+ return EXIT_SUCCESS;
+}
+
+
+
+static int
+gmp_Hermite_reduce_step(gmp_matrix * L, gmp_matrix * U,
+ size_t col, inverted_flag right_inverted)
+{
+
+ size_t i, j;
+ mpz_t pivot, elem;
+ mpz_t cff;
+
+ mpz_init(pivot);
+ mpz_init(elem);
+ mpz_init(cff);
+
+ if(col > (L->rows))
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* printf("Col = %i\n", col);
+ printf("L before\n");
+ gmp_matrix_printf(L);*/
+
+ for(j = col-1; j >= 1; j--)
+ {
+ for(i = j+1; i <= col; i++)
+ {
+ gmp_matrix_get_elem(pivot, i, i, L);
+ gmp_matrix_get_elem(elem, i, j, L);
+ /* printf(" i %i j %i\n",i,j); */
+
+ if((mpz_cmpabs(elem, pivot) >= 0) || (mpz_sgn(elem) < 0))
+ {
+ /* The objective:
+ 0 <= elem + k*pivot < pivot
+ -elem <= k*pivot < pivot - elem
+ -elem/pivot <= k < - elem/pivot + 1
+
+ Solution:
+ k = ceil(-elem/pivot)
+ */
+
+ /* Make cff = -elem */
+ mpz_neg (cff, elem);
+ mpz_cdiv_q (cff, cff, pivot);
+
+ /* printf("col %i j %i\n", i, j);
+ printf("elem %i k %i pivot %i\n",
+ mpz_get_si(elem),
+ mpz_get_si(cff),
+ mpz_get_si(pivot));*/
+
+
+ gmp_matrix_add_col(cff, i, j, L);
+
+ /* Update the unimodular factor matrix */
+ if(right_inverted == INVERTED)
+ {
+ gmp_matrix_add_col(cff, i, j, U);
+ }
+
+ /* [1 0] [ 1 0] = [1 0]
+ [cff 1] [-cff 1] [0 1] */
+ else
+ {
+ mpz_neg (cff, cff);
+ gmp_matrix_add_row(cff, j, i, U);
+ }
+
+ /* printf("\n");
+ gmp_matrix_printf(L);*/
+
+ }
+ }
+ }
+
+/* printf("L after\n"); */
+/* gmp_matrix_printf(L); */
+
+ mpz_clear(pivot);
+ mpz_clear(elem);
+ mpz_clear(cff);
+ return EXIT_SUCCESS;
+}
+
+
+static int
+gmp_normal_form_make_Hermite(gmp_normal_form * nf)
+{
+ size_t rows, cols;
+ size_t pivot_ind;
+ size_t schur, colind;
+ mpz_t pivot;
+ gmp_matrix * canonical;
+ gmp_matrix * left;
+ gmp_matrix * right;
+
+ if(nf == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* OK, it's safe to get to business */
+ canonical = nf->canonical;
+ left = nf->left;
+ right = nf->right;
+ rows = canonical -> rows;
+ cols = canonical -> cols;
+
+ mpz_init(pivot);
+
+ /* "schur" denotes where the present Schur complement starts */
+ schur = 1;
+
+ while((schur <= rows) && (schur <= cols))
+ {
+ /* Eliminate a column */
+ if (schur > 1)
+ {
+ gmp_Hermite_eliminate_step(canonical, right,
+ schur, nf->right_inverted);
+ }
+
+ /* Find a pivot */
+ pivot_ind = gmp_matrix_col_inz(schur, rows, schur, canonical);
+
+
+ /* If no nonzeros was found, the column is all zero, hence
+ settled with. Permute it to the end and decrement cols. */
+ if(pivot_ind == 0)
+ {
+ gmp_matrix_swap_cols(schur, cols, canonical);
+
+ if(nf -> right_inverted == INVERTED)
+ {
+ gmp_matrix_swap_cols(schur, cols, right);
+ }
+ else
+ {
+ gmp_matrix_swap_rows(schur, cols, right);
+ }
+
+ cols--;
+
+ /* When the whole column was zeroed, the diagonal
+ elements may have got reduced. Reduce the sub-
+ diagonals as well*/
+
+ if(schur > 1)
+ {
+ gmp_Hermite_reduce_step (canonical, right, schur-1,
+ nf -> right_inverted);
+ }
+ }
+
+ /* A nonzero pivot was found. Permute it to the diagonal position,
+ make it positive, and reduce the off-diagonals.
+ The schur complement now starts from the next diagonal. */
+ else
+ {
+ pivot_ind = schur+pivot_ind-1;
+ gmp_matrix_swap_rows(schur, pivot_ind, canonical);
+
+ if(nf->left_inverted == INVERTED)
+ {
+ gmp_matrix_swap_rows(schur, pivot_ind, left);
+ }
+ else
+ {
+ gmp_matrix_swap_cols(schur, pivot_ind, left);
+ }
+
+ /* Make the pivot positive */
+ gmp_matrix_get_elem(pivot, schur, schur, canonical);
+
+ if(mpz_cmp_si(pivot, 0) < 0)
+ {
+ gmp_matrix_negate_col(schur, canonical);
+
+ if(nf->right_inverted == INVERTED)
+ {
+ gmp_matrix_negate_col(schur, right);
+ }
+ else
+ {
+ gmp_matrix_negate_row(schur, right);
+ }
+ }
+
+ /* Reduce off-diagonals */
+ gmp_Hermite_reduce_step (canonical, right, schur, nf -> right_inverted);
+
+ schur++;
+ }
+ }
+
+ /* The Schur complement is now empty. There may still be uneliminated
+ columns left (in case of a wide matrix) */
+
+ colind = schur;
+ while (colind <= cols)
+ {
+ gmp_Hermite_eliminate_step (canonical, right, colind, nf->right_inverted);
+ gmp_Hermite_reduce_step (canonical, right, schur-1, nf->right_inverted);
+ colind++;
+ }
+
+ mpz_clear(pivot);
+
+ return EXIT_SUCCESS;
+}
+
+
+
+gmp_normal_form *
+create_gmp_Hermite_normal_form(gmp_matrix * A,
+ inverted_flag left_inverted,
+ inverted_flag right_inverted)
+{
+ gmp_normal_form * new_nf;
+
+ if(A == NULL)
+ {
+ return NULL;
+ }
+
+ new_nf =
+ create_gmp_trivial_normal_form(A, left_inverted, right_inverted);
+
+ if(new_nf == NULL)
+ {
+ /* A has been destroyed already */
+ return NULL;
+ }
+
+ if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ return NULL;
+ }
+
+ return new_nf;
+}
+
+
+gmp_normal_form *
+create_gmp_Smith_normal_form(gmp_matrix * A,
+ inverted_flag left_inverted,
+ inverted_flag right_inverted)
+{
+
+ gmp_normal_form * new_nf;
+ gmp_matrix * canonical;
+ gmp_matrix * elbow;
+ inverted_flag elbow_flag;
+ size_t rows, cols;
+ size_t i, j;
+ size_t subd_ind, row_undivisible;
+ size_t last_ready_row, first_ready_col, ind;
+ mpz_t pivot;
+ mpz_t remainder;
+
+ if(A == NULL)
+ {
+ return NULL;
+ }
+
+ new_nf =
+ create_gmp_trivial_normal_form(A, left_inverted, right_inverted);
+
+ if(new_nf == NULL)
+ {
+ /* A has been destroyed already */
+ return NULL;
+ }
+
+ canonical = new_nf -> canonical;
+ mpz_init(pivot);
+ mpz_init(remainder);
+ rows = canonical->rows;
+ cols = canonical->cols;
+ last_ready_row = 0;
+ first_ready_col = (cols < rows) ? (cols+1) : (rows+1);
+
+ while(first_ready_col > last_ready_row + 1)
+ {
+ /*******/
+ /* HNF */
+ /*******/
+
+ if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+
+#ifdef DEBUG
+ printf("HNF\n");
+ gmp_matrix_printf (new_nf -> left);
+ gmp_matrix_printf (new_nf -> canonical);
+ gmp_matrix_printf (new_nf -> right);
+#endif
+
+ /**********************/
+ /* Find ready columns */
+ /**********************/
+
+ /* If a diagonal entry is zero, so is the corresponding
+ column. The zero diagonals always reside in the end.
+ Seek until zero diagonal encountered, but stay within the matrix! */
+ ind = 1;
+ while ( (mpz_cmp_si(canonical -> storage[(ind-1)+(ind-1)*rows], 0)
+ != 0)
+ &&
+ (ind < first_ready_col) )
+ {
+ ind ++;
+ }
+ first_ready_col = ind;
+
+ /* Note: The number of ready cols is settled after the first HNF,
+ but the check is cheap. */
+
+ /**********************************************/
+ /* Permute unit diagonals such that they lead */
+ /**********************************************/
+
+ /* If the recently computed HNF has ones on the diagonal, their
+ corresponding rows are all zero (except the diagonal).
+ They are then settled, because the next LHNF kills the elements
+ on their columns. */
+
+ ind = last_ready_row+1;
+
+ /* Stay within the nonzero cols of the matrix */
+ while (ind < first_ready_col)
+ {
+ /* Unit diagonal encountered */
+ if(mpz_cmp_si ( canonical->storage[(ind-1) + (ind-1)*rows],
+ 1) == 0)
+ {
+ /* If not in the beginning, permute to extend the leading minor
+ with unit diagonals */
+ if(ind != last_ready_row+1)
+ {
+ gmp_matrix_swap_rows(last_ready_row+1, ind,
+ new_nf -> canonical);
+
+ if(left_inverted == INVERTED)
+ {
+ gmp_matrix_swap_rows(last_ready_row+1, ind,
+ new_nf -> left);
+ }
+ else
+ {
+ gmp_matrix_swap_cols(last_ready_row+1, ind,
+ new_nf -> left);
+ }
+
+ gmp_matrix_swap_cols(last_ready_row+1, ind,
+ new_nf -> canonical);
+
+ if(right_inverted == INVERTED)
+ {
+ gmp_matrix_swap_cols(last_ready_row+1, ind,
+ new_nf -> right);
+ }
+ else
+ {
+ gmp_matrix_swap_rows(last_ready_row+1, ind,
+ new_nf -> right);
+ }
+ }
+ last_ready_row ++;
+ }
+ ind++;
+ }
+
+#ifdef DEBUG
+ printf("Leading units\n");
+ gmp_matrix_printf (new_nf -> left);
+ gmp_matrix_printf (new_nf -> canonical);
+ gmp_matrix_printf (new_nf -> right);
+#endif
+
+ /********/
+ /* LHNF */
+ /********/
+
+ /* Extravagant strategy: Compute HNF of an explicit tranpose. */
+
+ /* 1. Transpose everything in the normal form */
+ if(gmp_matrix_transp(canonical) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->left) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->right) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+
+ elbow = new_nf ->left;
+ new_nf ->left = new_nf ->right;
+ new_nf ->right = elbow;
+ elbow_flag = new_nf ->left_inverted;
+ new_nf ->left_inverted = new_nf ->right_inverted;
+ new_nf ->right_inverted = elbow_flag;
+
+ /* 2. Compute HNF of the transpose */
+ if(gmp_normal_form_make_Hermite(new_nf) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+
+ /* 3. Transpose everything back */
+ if(gmp_matrix_transp(canonical) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->left) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ if(gmp_matrix_transp(new_nf->right) != EXIT_SUCCESS)
+ {
+ destroy_gmp_normal_form(new_nf);
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return NULL;
+ }
+ elbow = new_nf ->left;
+ new_nf ->left = new_nf ->right;
+ new_nf ->right = elbow;
+ elbow_flag = new_nf ->left_inverted;
+ new_nf ->left_inverted = new_nf ->right_inverted;
+ new_nf ->right_inverted = elbow_flag;
+
+
+#ifdef DEBUG
+ printf("LHNF\n");
+ gmp_matrix_printf (new_nf -> left);
+ gmp_matrix_printf (new_nf -> canonical);
+ gmp_matrix_printf (new_nf -> right);
+#endif
+
+
+ /*****************************************************/
+ /* Check if more of the leading normal form is ready */
+ /*****************************************************/
+
+ /* The matrix is in LHNF, i.e. it is upper triangular.
+ If the row trailing the top left diagonal element is
+ zero, the diagonal element may be an invariant factor
+ on its final position, and the stage may be ready.
+
+ The stage may not still be ready: The leading diagonal element
+ of D may not divide the rest of the Schur complement. */
+
+ subd_ind = 0;
+ row_undivisible = 0;
+
+ /* Explanation of loop conditions:
+ 1.) No relative primes found from Schur complement
+ 2.) Stay within the Schur complement
+ 3.) If a nonzero is found from the trailing row, the stage is
+ definitely not ready */
+ while (row_undivisible == 0 &&
+ last_ready_row + 1 < first_ready_col &&
+ subd_ind == 0)
+ {
+ subd_ind =
+ gmp_matrix_row_inz(last_ready_row+1,
+ last_ready_row+2, cols,
+ canonical);
+
+ /* printf("subd_ind %i\n", subd_ind);
+ printf("last_ready_row %i\n", last_ready_row); */
+
+ /* No nonzeros found, the stage may be ready */
+ if (subd_ind == 0)
+ {
+ mpz_set (pivot,
+ canonical->storage[(last_ready_row)+
+ (last_ready_row)*rows]);
+
+ /* Check whether the pivot divides all elements in the Schur
+ complement */
+ row_undivisible = 0;
+ for(j = last_ready_row+2; j < first_ready_col; j++)
+ {
+ for(i = last_ready_row+2; i <= j; i++)
+ {
+ mpz_tdiv_r (remainder,
+ canonical->storage[(i-1)+
+ (j-1)*rows],
+ pivot);
+
+ if(mpz_cmp_si(remainder, 0) !=0)
+ {
+ row_undivisible = i;
+ /* col_undivisible = j; */
+ }
+ }
+ }
+
+ /* printf("Row undivisible %i\n", row_undivisible); */
+
+ /* If a relative prime was found from the Schur complement,
+ add that row to the first row of the Schur complement */
+ if(row_undivisible != 0)
+ {
+ mpz_set_si (remainder, 1);
+ gmp_matrix_add_row(remainder,
+ row_undivisible, last_ready_row+1,
+ canonical);
+
+ /* [ 1 0] [1 0] = [1 0]
+ [-1 1] [1 1] [0 1] */
+
+ if(left_inverted == INVERTED)
+ {
+ gmp_matrix_add_row(remainder,
+ row_undivisible, last_ready_row+1,
+ new_nf->left);
+ }
+ else
+ {
+ mpz_neg (remainder, remainder);
+ gmp_matrix_add_col(remainder,
+ last_ready_row+1, row_undivisible,
+ new_nf->left);
+ }
+ }
+ /* The Schur complement is smaller by one again */
+ else
+ {
+ last_ready_row++;
+ }
+ }
+ }
+ } /* The main loop ends here */
+
+ mpz_clear(pivot);
+ mpz_clear(remainder);
+ return new_nf;
+}
+
+int destroy_gmp_normal_form(gmp_normal_form * nf)
+{
+ int status;
+
+ if(nf == NULL)
+ {
+ return EXIT_FAILURE;
+ }
+
+ status = EXIT_SUCCESS;
+ if(destroy_gmp_matrix(nf -> left) == EXIT_FAILURE)
+ {
+ status = EXIT_FAILURE;
+ }
+
+ if(destroy_gmp_matrix(nf -> canonical) == EXIT_FAILURE)
+ {
+ status = EXIT_FAILURE;
+ }
+
+ if(destroy_gmp_matrix(nf -> right) == EXIT_FAILURE)
+ {
+ status = EXIT_FAILURE;
+ }
+ free(nf);
+
+ return status;
+}
+
+
diff --git a/contrib/kbipack/gmp_normal_form.h b/contrib/kbipack/gmp_normal_form.h
new file mode 100755
index 0000000000000000000000000000000000000000..3250646a6505e1cda3f0b293132ebfeb5d2afaf7
--- /dev/null
+++ b/contrib/kbipack/gmp_normal_form.h
@@ -0,0 +1,154 @@
+/*
+ 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
diff --git a/contrib/kbipack/gmp_normal_form.o b/contrib/kbipack/gmp_normal_form.o
new file mode 100755
index 0000000000000000000000000000000000000000..f165a904233fc099325f1bb789f948c5207669de
Binary files /dev/null and b/contrib/kbipack/gmp_normal_form.o differ
diff --git a/contrib/kbipack/kbihnf.m b/contrib/kbipack/kbihnf.m
new file mode 100755
index 0000000000000000000000000000000000000000..f04977ca950ecadde14a1a85a42e00419e717485
--- /dev/null
+++ b/contrib/kbipack/kbihnf.m
@@ -0,0 +1,58 @@
+function [P, L, U] = kbihnf(A)
+
+% [P, L, U] = kbihnf(A)
+%
+% Computes the integer Hermite normal form L such that
+%
+% P*A = L*U
+%
+% Uses external Kannan-Bachem implementation
+
+% 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
+
+ [i,j,v] = find(A);
+
+ fid = fopen('hnf_temp.crd', 'w');
+ fprintf(fid, '%g %g %g\n', size(A,1), size(A,2), length(i));
+ if ~isempty(v)
+ fprintf(fid, '%g %g %g\n', [i'; j'; v']);
+ end
+ fclose(fid);
+
+ system ("./compute_normal_form -Hl hnf_temp.crd");
+
+ load hnf_temp.left;
+ P = zeros(hnf_temp(1,1), hnf_temp(1,2));
+ for i = 2:size(hnf_temp,1)
+ P(hnf_temp(i,1), hnf_temp(i,2)) = hnf_temp(i,3);
+ end
+
+ load hnf_temp.can;
+ L = zeros(hnf_temp(1,1), hnf_temp(1,2));
+ for i = 2:size(hnf_temp,1)
+ L(hnf_temp(i,1), hnf_temp(i,2)) = hnf_temp(i,3);
+ end
+
+ load hnf_temp.right;
+ U = zeros(hnf_temp(1,1), hnf_temp(1,2));
+ for i = 2:size(hnf_temp,1)
+ U(hnf_temp(i,1), hnf_temp(i,2)) = hnf_temp(i,3);
+ end
diff --git a/contrib/kbipack/kbisnf.m b/contrib/kbipack/kbisnf.m
new file mode 100755
index 0000000000000000000000000000000000000000..066509dceb1ce1ac06c3a9c78b38a798801de212
--- /dev/null
+++ b/contrib/kbipack/kbisnf.m
@@ -0,0 +1,58 @@
+function [U, S, V] = kbisnf(A)
+
+% [U, S, V] = kbisnf(A)
+%
+% Computes the integer Smith normal form S such that
+%
+% A = U*S*V
+%
+% Uses external Kannan-Bachem implementation
+
+% Copyright (C) 4.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
+
+ [i,j,v] = find(A);
+
+ fid = fopen('snf_temp.crd', 'w');
+ fprintf(fid, '%g %g %g\n', size(A,1), size(A,2), length(i));
+ if ~isempty(v)
+ fprintf(fid, '%g %g %g\n', [i'; j'; v']);
+ end
+ fclose(fid);
+
+ system ('./compute_normal_form -S snf_temp.crd');
+
+ load snf_temp.left;
+ U = zeros(snf_temp(1,1), snf_temp(1,2));
+ for i = 2:size(snf_temp,1)
+ U(snf_temp(i,1), snf_temp(i,2)) = snf_temp(i,3);
+ end
+
+ load snf_temp.can;
+ S = zeros(snf_temp(1,1), snf_temp(1,2));
+ for i = 2:size(snf_temp,1)
+ S(snf_temp(i,1), snf_temp(i,2)) = snf_temp(i,3);
+ end
+
+ load snf_temp.right;
+ V = zeros(snf_temp(1,1), snf_temp(1,2));
+ for i = 2:size(snf_temp,1)
+ V(snf_temp(i,1), snf_temp(i,2)) = snf_temp(i,3);
+ end
diff --git a/contrib/kbipack/licence b/contrib/kbipack/licence
new file mode 100755
index 0000000000000000000000000000000000000000..5b6e7c66c276e7610d4a73c70ec1a1f7c1003259
--- /dev/null
+++ b/contrib/kbipack/licence
@@ -0,0 +1,340 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ 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
+ (at your option) 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
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Library General
+Public License instead of this License.
diff --git a/contrib/kbipack/readme b/contrib/kbipack/readme
new file mode 100755
index 0000000000000000000000000000000000000000..3fbf20be8b076ec4289724b3646d011d0412cda7
--- /dev/null
+++ b/contrib/kbipack/readme
@@ -0,0 +1,106 @@
+-----------
+KBIPack 1.0
+-----------
+
+This is an ANSI C-implementation of Kannan-Bachem algorithms for Hermite
+and Smith normal forms for integer matrices. I wrote it to be the core of
+an integer homology group solver. You are welcome to improve this package,
+but I am unfortunately unable to invest much more time to it myself. The
+license is GNU GPL (see below).
+
+The algorithms are Kannan - Bachem algorithms with improvement by Chou and
+Collins. The Smith normal form routine expects a large number of unit
+invariant factors (this occurs in my application) 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/
+
+Technicalities:
+---------------
+The package relies on GNU multiple precision library gmp to represent large
+integers. You must have gmp installed before you compile.
+
+It should run smoothly on UNIX systems with working standard I/O, and
+in Cygwin (http://www.cygwin.com/).
+
+You can use this package in at least three ways:
+1) You can call the functions bkihnf and bkisnf in respective .m-files
+ from MatLab or Octave. Note that the numbers in the results may be too
+ large to be exactly read back by them.
+2) You can use the executable "compute normal form" just like any command.
+ Use the -h flag to get info.
+3) You can use just the core routines from your own application.
+
+Contents:
+---------
+bkihnf.m Integer Hermite normal form for Octave/MatLab
+bkisnf.m Integer Smith normal form for Octave/MatLab
+compute_normal_form.c A main program for a comand-line version
+gmp_blas.c Low-level routines for integer arrays
+gmp_blas.h
+gmp_matrix.c Integer matrix type
+gmp_matrix.h
+gmp_matrix_io.c I/O for integer matrices
+gmp_matrix_io.h
+gmp_normal_form.c The core computation routine
+gmp_normal_form.h
+LICENCE The GNU General Public License version 2
+Makefile A rudimentary Makefile (just type "make")
+README This file
+
+Installation:
+-------------
+First install gmp (http://www.swox.com/gmp/ or a binary package).
+If you use gcc, just type "make". If not, alter the Makefile as necessary.
+
+Licensing:
+----------
+Copyright (C) 2005 Saku Suuriniemi
+
+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
+http://www.em.tut.fi/Eng/
+
+Background:
+-----------
+The package is a part of my PhD project, so if you find this package useful
+and wish to refer to my thesis (also available on my homepage), here is the
+BiBTeX entry:
+
+@phdthesis{Suuriniemi:PhD2004,
+author="Saku Suuriniemi",
+school="Tampere University of Technology",
+title="Homological computations in electromagnetic modeling",
+address="Tampere",
+year=2004,
+isbn="952-15-1237-7"
+}
+
+Thanks:
+-------
+Thanks are due to Chistophe Geuzaine for testing and comments.
diff --git a/utils/misc/celldriver.cpp b/utils/misc/celldriver.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..c8b91e03a96f84cd6470f3f13df467e9ddb7dc36
--- /dev/null
+++ b/utils/misc/celldriver.cpp
@@ -0,0 +1,63 @@
+// configure, compile and install Gmsh as a library with
+//
+// ./configure --disable-gui --disable-netgen --disable-chaco
+// --disable-metis --disable-tetgen --prefix=/usr/local/
+// make install-lib
+//
+// Then compile this driver with "g++ kaka.cpp -lGmsh -llapack -lblas -lgmp"
+
+#include <stdio.h>
+#include <gmsh/Gmsh.h>
+#include <gmsh/GModel.h>
+#include <gmsh/MElement.h>
+#include <gmsh/CellComplex.h>
+
+
+int main(int argc, char **argv)
+{
+ GmshInitialize(argc, argv);
+ GModel *m = new GModel();
+ m->readGEO("transformer.geo");
+ m->mesh(3);
+ m->writeMSH("transformer.msh");
+
+ printf("This model has: %d GRegions, %d GFaces, %d GEdges and %d GVertices. \n" , m->getNumRegions(), m->getNumFaces(), m->getNumEdges(), m->getNumVertices());
+
+ std::vector<GEntity*> domain;
+ std::vector<GEntity*> subdomain;
+
+ // whole model
+ for(GModel::riter rit = m->firstRegion(); rit != m->lastRegion(); rit++){
+ GEntity *r = *rit;
+ domain.push_back(r);
+ }
+
+ // the ports
+ GModel::fiter sub = m->firstFace();
+ GEntity *s = *sub;
+ subdomain.push_back(s);
+ s = *(++sub);
+ subdomain.push_back(s);
+ s =*(++sub);
+ subdomain.push_back(s);
+ s= *(++sub);
+ subdomain.push_back(s);
+
+ CellComplex complex = CellComplex(domain, subdomain);
+
+ printf("Cell complex of this model has: %d volumes, %d faces, %d edges and %d vertices\n",
+ complex.getSize(3), complex.getSize(2), complex.getSize(1), complex.getSize(0));
+
+ complex.reduceComplex();
+ complex.writeComplexMSH("reduced_complex.msh");
+ complex.coreduceComplex();
+ complex.writeComplexMSH("coreduced_complex.msh");
+
+ delete m;
+ GmshFinalize();
+
+}
+
+
+
+
diff --git a/utils/misc/transformer.geo b/utils/misc/transformer.geo
new file mode 100755
index 0000000000000000000000000000000000000000..9021d10d66be327e345dc8d8fec4347d35a15fbc
--- /dev/null
+++ b/utils/misc/transformer.geo
@@ -0,0 +1,99 @@
+// Gmsh project created on Thu Mar 26 10:01:45 2009
+
+m=0.5;
+
+Point(newp) = {0, 0, 0, m};
+Point(newp) = {10, 0, 0, m};
+Point(newp) = {10, 10, 0, m};
+Point(newp) = {0, 10, 0, m};
+Point(newp) = {4, 4, 0, m};
+Point(newp) = {6, 4, 0, m};
+Point(newp) = {6, 6, 0, m};
+Point(newp) = {4, 6, 0, m};
+
+Point(newp) = {2, 0, 0, m};
+Point(newp) = {8, 0, 0, m};
+Point(newp) = {2, 10, 0, m};
+Point(newp) = {8, 10, 0, m};
+
+Point(newp) = {0, 0, 1, m};
+Point(newp) = {10, 0, 1, m};
+Point(newp) = {10, 10, 1, m};
+Point(newp) = {0, 10, 1, m};
+Point(newp) = {4, 4, 1, m};
+Point(newp) = {6, 4, 1, m};
+Point(newp) = {6, 6, 1, m};
+Point(newp) = {4, 6, 1, m};
+
+Point(newp) = {2, 0, 1, m};
+Point(newp) = {8, 0, 1, m};
+Point(newp) = {2, 10, 1, m};
+Point(newp) = {8, 10, 1, m};
+Line(1) = {16, 23};
+Line(2) = {23, 11};
+Line(3) = {11, 4};
+Line(4) = {4, 16};
+Line(5) = {24, 12};
+Line(6) = {12, 3};
+Line(7) = {3, 15};
+Line(8) = {15, 24};
+Line(9) = {10, 2};
+Line(10) = {2, 14};
+Line(11) = {14, 22};
+Line(12) = {22, 10};
+Line(13) = {21, 9};
+Line(14) = {9, 1};
+Line(15) = {1, 13};
+Line(16) = {13, 21};
+Line Loop(17) = {3, 4, 1, 2};
+Ruled Surface(18) = {17};
+Line Loop(19) = {6, 7, 8, 5};
+Ruled Surface(20) = {19};
+Line Loop(21) = {9, 10, 11, 12};
+Ruled Surface(22) = {21};
+Line Loop(23) = {14, 15, 16, 13};
+Ruled Surface(24) = {23};
+Line(25) = {16, 13};
+Line(26) = {1, 4};
+Line(27) = {11, 12};
+Line(28) = {24, 23};
+Line(29) = {21, 22};
+Line(30) = {10, 9};
+Line(31) = {2, 3};
+Line(32) = {15, 14};
+Line(33) = {20, 19};
+Line(34) = {19, 18};
+Line(35) = {18, 17};
+Line(36) = {17, 20};
+Line(37) = {8, 7};
+Line(38) = {7, 6};
+Line(39) = {6, 18};
+Line(40) = {5, 6};
+Line(41) = {5, 8};
+Line(42) = {20, 8};
+Line(43) = {17, 5};
+Line(44) = {19, 7};
+Line Loop(45) = {27, -5, 28, 2};
+Ruled Surface(46) = {45};
+Line Loop(47) = {25, -15, 26, 4};
+Ruled Surface(48) = {47};
+Line Loop(49) = {29, 12, 30, -13};
+Ruled Surface(50) = {49};
+Line Loop(51) = {32, -10, 31, 7};
+Ruled Surface(52) = {51};
+Line Loop(53) = {41, -42, -36, 43};
+Ruled Surface(54) = {53};
+Line Loop(55) = {35, 43, 40, 39};
+Ruled Surface(56) = {55};
+Line Loop(57) = {34, -39, -38, -44};
+Ruled Surface(58) = {57};
+Line Loop(59) = {33, 44, -37, -42};
+Ruled Surface(60) = {59};
+Line Loop(61) = {28, -1, 25, 16, 29, -11, -32, 8};
+Line Loop(62) = {33, 34, 35, 36};
+Ruled Surface(63) = {61, 62};
+Line Loop(64) = {3, -26, -14, -30, 9, 31, -6, -27};
+Line Loop(65) = {37, 38, -40, 41};
+Ruled Surface(66) = {64, 65};
+Surface Loop(67) = {63, 46, 66, 18, 48, 24, 50, 22, 52, 20, 54, 60, 58, 56};
+Volume(68) = {67};