diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 2980645832fd5aa74df3a127a4a145ac30431d46..ce80a05b9b53d6fe7390286a4dec9bf3b7174ad1 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -290,7 +290,7 @@ int CellComplex::coreduction(int dim){
removeCell(cell);
removeCell(bd_c.at(0));
count++;
- coreduced =true;
+ coreduced = true;
}
}
@@ -317,7 +317,7 @@ int CellComplex::reduction(int dim){
removeCell(cell);
removeCell(cbd_c.at(0));
count++;
- reduced =true;
+ reduced = true;
}
}
@@ -335,130 +335,86 @@ void CellComplex::reduceComplex(){
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){
+ std::set<Cell*, Less_Cell> removedCells[4];
- // h_dim: C_dim -> C_(dim-1)
-
-
- if(dim > 3 || dim < 1){
- return;
- }
-
- destroy_gmp_matrix(_HMatrix[dim]);
-
- if( _cells[dim].size() == 0 ){
- _HMatrix[dim] = create_gmp_matrix_zero(1, 1);
- //gmp_matrix_printf(_HMatrix[dim]);
- return;
- }
- unsigned int cols = _cells[dim].size();
-
- if( _cells[dim-1].size() == 0){
- _HMatrix[dim] = create_gmp_matrix_zero(1, cols);
- //gmp_matrix_printf(_HMatrix[dim]);
- return;
- }
- unsigned int rows = _cells[dim-1].size();
-
- mpz_t elems[rows*cols];
- mpz_array_init( elems[0], rows*cols, sizeof(mpz_t) );
-
- 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);
+ printf("Cell complex before coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+ getSize(3), getSize(2), getSize(1), getSize(0));
+ for(int i = 0; i < 3; i++){
+ while (getSize(i) != 0){
+ citer cit = firstCell(i);
+ Cell* cell = *cit;
+ removedCells[i].insert(cell);
+ _cells[i].erase(cell);
+ //removeCell(cell);
+ //complex.coreductionMrozek(kaka);
+ coreduction(0);
+ coreduction(1);
+ coreduction(2);
}
- mpz_set_ui(elems[i], kappa(*high, *low));
- low++;
+
}
- */
- /*
- unsigned int i = 0;
- while(i < rows*cols){
- while(low != lastCell(dim-1)){
- mpz_set_si(elems[i], kappa(*high, *low));
- //printf(" %d %d, ",i, kappa(*high, *low));
- i++;
- low++;
+ for(int i = 0; i < 3; i++){
+ for(citer cit = removedCells[i].begin(); cit != removedCells[i].end(); cit++){
+ Cell* cell = *cit;
+ _cells[i].insert(cell);
}
- low = firstCell(dim-1);
- high++;
}
-
- _HMatrix[dim] = create_gmp_matrix(rows, cols,(const mpz_t *) elems);
-
- //gmp_matrix_printf(_HMatrix[dim]);
-
- return;
+ printf("Cell complex after coreduction: %d volumes, %d faces, %d edges and %d vertices.\n",
+ getSize(3), getSize(2), getSize(1), getSize(0));
}
-*/
+
std::vector<gmp_matrix*> CellComplex::constructHMatrices(){
// h_dim: C_dim -> C_(dim-1)
std::vector<gmp_matrix*> HMatrix;
-
- HMatrix.push_back(create_gmp_matrix_zero(1, 1));
-
+
+ HMatrix.push_back(NULL);
+
for(int dim = 1; dim < 4; dim++){
- if( _cells[dim].size() == 0 ){
- HMatrix.push_back( create_gmp_matrix_zero(1, 1));
- //gmp_matrix_printf(HMatrix[dim]);
- break;
- }
- unsigned int cols = _cells[dim].size();
+ unsigned int cols = _cells[dim].size();
+ unsigned int rows = _cells[dim-1].size();
- if( _cells[dim-1].size() == 0){
- HMatrix.push_back(create_gmp_matrix_zero(1, cols));
- //gmp_matrix_printf(HMatrix[dim]);
- break;
+ for(citer cit = firstCell(dim); cit != lastCell(dim); cit++){
+ Cell* cell = *cit;
+ if(cell->inSubdomain()) cols--;
+ }
+ for(citer cit = firstCell(dim-1); cit != lastCell(dim-1); cit++){
+ Cell* cell = *cit;
+ if(cell->inSubdomain()) rows--;
}
- unsigned int rows = _cells[dim-1].size();
- mpz_t elems[rows*cols];
- mpz_array_init( elems[0], rows*cols, sizeof(mpz_t) );
- citer high = firstCell(dim);
- citer low = firstCell(dim-1);
-
- unsigned int i = 0;
- while(i < rows*cols){
- while(low != lastCell(dim-1)){
- mpz_set_si(elems[i], kappa(*high, *low));
- //printf(" %d %d, ",i, kappa(*high, *low));
- i++;
- low++;
+ if( cols == 0 ){
+ HMatrix.push_back(NULL);
+ }
+ else if( rows == 0){
+ HMatrix.push_back(create_gmp_matrix_zero(1, cols));
+ }
+ else{
+ long int elems[rows*cols];
+
+ citer high = firstCell(dim);
+ citer low = firstCell(dim-1);
+
+ unsigned int i = 0;
+ while(i < rows*cols){
+ while(low != lastCell(dim-1)){
+ Cell* lowcell = *low;
+ Cell* highcell = *high;
+ if(!(highcell->inSubdomain() || lowcell->inSubdomain())){
+ elems[i] = kappa(*high, *low);
+ i++;
+ }
+ low++;
+ }
+ low = firstCell(dim-1);
+ high++;
}
- low = firstCell(dim-1);
- high++;
+ HMatrix.push_back(create_gmp_matrix_int(rows, cols, elems));
+
}
- HMatrix.push_back(create_gmp_matrix(rows, cols,(const mpz_t *) elems));
-
- //gmp_matrix_printf(_HMatrix[dim]);
}
return HMatrix;
}
@@ -575,29 +531,112 @@ int CellComplex::writeComplexMSH(const std::string &name){
}
-void CellComplex::printComplex(int dim){
+void CellComplex::printComplex(int dim, bool subcomplex){
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));
+ if(!subcomplex && !cell->inSubdomain()) printf("%d ", cell->getVertex(i));
}
printf("\n");
}
}
-gmp_matrix* ChainComplex::ker(gmp_matrix* HMatrix){
+void ChainComplex::KerCod(int dim){
+
+ if(dim < 1 || dim > 3 || _HMatrix[dim] == NULL) return;
+
+ gmp_matrix* HMatrix = copy_gmp_matrix(_HMatrix[dim], 1, 1,
+ gmp_matrix_rows(_HMatrix[dim]), gmp_matrix_cols(_HMatrix[dim]));
+
+ gmp_normal_form* normalForm = create_gmp_Hermite_normal_form(HMatrix, NOT_INVERTED, INVERTED);
+ //printMatrix(normalForm->left);
+ //printMatrix(normalForm->canonical);
+ //printMatrix(normalForm->right);
+
+ int minRowCol = std::min(gmp_matrix_rows(normalForm->canonical), gmp_matrix_cols(normalForm->canonical));
+ int rank = 0;
+ mpz_t elem;
+ mpz_init(elem);
+
+ while(rank < minRowCol){
+ gmp_matrix_get_elem(elem, rank+1, rank+1, normalForm->canonical);
+ if(mpz_cmp_si(elem,0) == 0) break;
+ rank++;
+ }
- // H = USV -> WH = US, W = inv(V)
- gmp_normal_form* W_USform = create_gmp_Smith_normal_form(HMatrix, NOT_INVERTED, INVERTED);
-
+ if(rank != gmp_matrix_cols(normalForm->canonical)){
+ _kerH[dim] = copy_gmp_matrix(normalForm->right, 1, rank+1,
+ gmp_matrix_rows(normalForm->right), gmp_matrix_cols(normalForm->right));
+ }
+ if(rank > 0){
+ _codH[dim] = copy_gmp_matrix(normalForm->canonical, 1, 1,
+ gmp_matrix_rows(normalForm->canonical), rank);
+ gmp_matrix_left_mult(normalForm->left, _codH[dim]);
+ }
- return W_USform->canonical;
+ mpz_clear(elem);
+ destroy_gmp_normal_form(normalForm);
+ return;
}
+/*
+//j:B->Z
+void ChainComplex::Inclusion(int dim){
+
+ if(dim < 1 || dim > 3 || _kerH[dim] == NULL || _codH[dim] == NULL) return;
+
+ int rows = gmp_matrix_rows(Zbasis);
+ int cols = gmp_matrix_cols(Zbasis);
+ if(rows < cols) return;
+
+ rows = gmp_matrix_rows(Bbasis);
+ cols = gmp_matrix_cols(Bbasis);
+ if(rows < cols) return;
+
+ gmp_matrix* Zbasis = copy_gmp_matrix(_ker[dim], 1, 1,
+ gmp_matrix_rows(_kerH[dim]), gmp_matrix_cols(_kerH[dim]));
+ gmp_matrix* Bbasis = copy_gmp_matrix(_cod[dim], 1, 1,
+ gmp_matrix_rows(_codH[dim]), gmp_matrix_cols(_codH[dim]));
+
+ // A*inv(V) = U*S
+ normalForm = create_gmp_Smith_normal_form(Zbasis, NOT_INVERTED, INVERTED);
+
+ mpz_t elem;
+ mpz_init(elem);
+
+ for(int i = 1; i < cols; i++){
+ gmp_matrix_get_elem(elem, i, i, normalForm->canonical);
+ if(mpz_cmp_si(elem,0) == 0){
+ destroy_gmp_normal_form(normalForm);
+ return;
+ }
+ }
+
+ gmp_matrix_left_mult(normalForm->left, Bbasis);
+
+}
+*/
+void ChainComplex::matrixTest(){
+
+ int rows = 3;
+ int cols = 6;
+
+ long int elems[rows*cols];
+ for(int i = 1; i<=rows*cols; i++) elems[i-1] = i;
+
+ gmp_matrix* matrix = create_gmp_matrix_int(rows, cols, elems);
+
+ gmp_matrix* copymatrix = copy_gmp_matrix(matrix, 3, 2, 3, 5);
+
+ printMatrix(matrix);
+ printMatrix(copymatrix);
+
+ return;
+}
#endif
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index b7f4dcc619eef0a2199f43137fe6485f195fc10c..5f73d241d72811ffdca7ae9cbbb2907a490bba0b 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -44,7 +44,7 @@ class Cell
int _cbdSize;
// whether this cell belongs to a subdomain
- // used to in relative homology computation
+ // used in relative homology computation
bool _inSubdomain;
public:
@@ -294,9 +294,19 @@ class ThreeSimplex : public Simplex
class Less_Cell{
public:
bool operator()(const Cell* c1, const Cell* c2) const {
+
+ // subcomplex cells in the end
+ //if(c1->inSubdomain() != c2->inSubdomain()){
+ // if(c1->inSubdomain()) return false;
+ // else return true;
+ //}
+
+ // cells with fever vertices first
if(c1->getNumVertices() != c2->getNumVertices()){
return (c1->getNumVertices() < c2->getNumVertices());
}
+
+ // "natural number" -like ordering (the number of a vertice corresponds a digit)
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;
@@ -327,7 +337,7 @@ class CellComplex
// sorted containers of unique cells in this cell complex
// one for each dimension
std::set<Cell*, Less_Cell> _cells[4];
-
+
// boundary operator matrices for this chain complex
// h_k: C_k -> C_(k-1)
//gmp_matrix* _HMatrix[4];
@@ -356,7 +366,7 @@ class CellComplex
public:
CellComplex( std::set<Cell*, Less_Cell>* cells ) {
for(int i = 0; i < 4; i++){
- _cells[i] = cells[i];
+ //_cells[i] = cells[i];
}
}
CellComplex( std::vector<GEntity*> domain, std::vector<GEntity*> subdomain );
@@ -409,7 +419,7 @@ class CellComplex
virtual int coreductionMrozek(Cell* generator);
// print the vertices of cells of certain dimension
- virtual void printComplex(int dim);
+ virtual void printComplex(int dim, bool subcomplex);
// write this cell complex in legacy .msh format
virtual int writeComplexMSH(const std::string &name);
@@ -419,39 +429,76 @@ class CellComplex
// get the boundary operator matrix dim->dim-1
//virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
+ // construct boundary operator matrices of this cell complex
+ // used to construct a chain complex
virtual std::vector<gmp_matrix*> constructHMatrices();
};
+// A class representing a chain complex of a cell complex.
+// This should only be constructed for a reduced cell complex because of
+// dense matrix reprentations and great computational complexity in its methods.
class ChainComplex{
private:
// boundary operator matrices for this chain complex
// h_k: C_k -> C_(k-1)
gmp_matrix* _HMatrix[4];
+ // Basis matrices for the kernels and codomains of the boundary operator matrices
+ gmp_matrix* _kerH[4];
+ gmp_matrix* _codH[4];
+
+ gmp_matrix* _JMatrix[4];
+ gmp_matrix* _QMatrix[4];
+ std::vector<int> _torsion[4];
+
public:
ChainComplex( std::vector<gmp_matrix*> HMatrix ){
for(int i = 0; i < HMatrix.size(); i++){
_HMatrix[i] = HMatrix.at(i);
- }
+ }
+ for(int i = 0; i < 4; i++){
+ _kerH[i] = NULL;
+ _codH[i] = NULL;
+ _JMatrix[i] = NULL;
+ _QMatrix[i] = NULL;
+ }
+
}
ChainComplex(){
for(int i = 0; i < 4; i++){
_HMatrix[i] = create_gmp_matrix_zero(1,1);
+ _kerH[i] = NULL;
+ _codH[i] = NULL;
+ _JMatrix[i] = NULL;
+ _QMatrix[i] = NULL;
}
}
virtual ~ChainComplex(){}
// get the boundary operator matrix dim->dim-1
virtual gmp_matrix* getHMatrix(int dim) { return _HMatrix[dim]; }
-
- virtual gmp_matrix* ker(gmp_matrix* HMatrix);
+ virtual gmp_matrix* getKerHMatrix(int dim) { return _kerH[dim]; }
+ virtual gmp_matrix* getCodHMatrix(int dim) { return _codH[dim]; }
+ virtual gmp_matrix* getJMatrix(int dim) { return _JMatrix[dim]; }
+ virtual gmp_matrix* getQMatrix(int dim) { return _QMatrix[dim]; }
+
+ // Compute basis for kernel and codomain of boundary operator matrix of dimension dim (ie. ker(h_dim) and cod(h_dim) )
+ virtual void KerCod(int dim);
+ // Compute matrix representation J for inclusion relation from dim-cells who are boundary of dim+1-cells
+ // to cycles of dim-cells (ie. j: cod(h_(dim+1)) -> ker(h_dim) )
+ virtual void Inclusion(int dim){}
+ // Compute quotient problem for given inclusion relation j to find representatives of homology groups
+ // and possible torsion coeffcients
+ virtual void Quotient(int dim){}
virtual int printMatrix(gmp_matrix* matrix){
printf("%d rows and %d columns\n", gmp_matrix_rows(matrix), gmp_matrix_cols(matrix));
return gmp_matrix_printf(matrix); }
+ // debugging aid
+ virtual void matrixTest();
};
diff --git a/contrib/kbipack/gmp_matrix.c b/contrib/kbipack/gmp_matrix.c
index 23d648c89eff75d10167c303d32132c6a5c9401b..2d07724b0804ac8f598a68b0ca9fc37842758892 100755
--- a/contrib/kbipack/gmp_matrix.c
+++ b/contrib/kbipack/gmp_matrix.c
@@ -22,7 +22,7 @@
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 $
+ $Id: gmp_matrix.c,v 1.2 2009-04-03 11:06:12 matti Exp $
*/
@@ -61,6 +61,38 @@ create_gmp_matrix(size_t r, size_t c,
return new_matrix;
}
+gmp_matrix *
+create_gmp_matrix_int(size_t r, size_t c,
+ const long int * e)
+{
+ gmp_matrix * new_matrix;
+ size_t ind;
+
+ new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix));
+ if(new_matrix == NULL)
+ {
+ return NULL;
+ }
+
+ new_matrix -> storage = (mpz_t *) calloc(r*c, sizeof(mpz_t));
+ if(new_matrix -> storage == NULL)
+ {
+ free(new_matrix);
+ return NULL;
+ }
+
+ new_matrix -> rows = r;
+ new_matrix -> cols = c;
+
+ for(ind = 0; ind < r*c; ind ++)
+ {
+ mpz_init (new_matrix -> storage[ind]);
+ mpz_set_si (new_matrix -> storage[ind], e[ind]);
+ }
+
+ return new_matrix;
+}
+
gmp_matrix *
create_gmp_matrix_identity(size_t dim)
@@ -128,6 +160,57 @@ create_gmp_matrix_zero(size_t rows, size_t cols)
return new_matrix;
}
+gmp_matrix *
+copy_gmp_matrix(const gmp_matrix * matrix,
+ const size_t start_row, const size_t start_col,
+ const size_t end_row, const size_t end_col)
+{
+ gmp_matrix * new_matrix;
+ size_t ind;
+ size_t r;
+ size_t c;
+ size_t old_rows;
+ size_t old_cols;
+ size_t i;
+ size_t j;
+
+ new_matrix = (gmp_matrix * ) malloc(sizeof(gmp_matrix));
+ if(new_matrix == NULL)
+ {
+ return NULL;
+ }
+
+ r = end_row-start_row+1;
+ c = end_col-start_col+1;
+ new_matrix -> storage = (mpz_t *) calloc(r*c, sizeof(mpz_t));
+ if(new_matrix -> storage == NULL)
+ {
+ free(new_matrix);
+ return NULL;
+ }
+
+ new_matrix -> rows = r;
+ new_matrix -> cols = c;
+
+ old_rows = matrix -> rows;
+ old_cols = matrix -> cols;
+
+ ind = 0;
+ for(j = 1; j <= old_cols; j++){
+ if(j >= start_col && j <= end_col){
+ for(i = 1; i <= old_rows; i++){
+ if(i >= start_row && i <= end_row){
+ mpz_init (new_matrix -> storage[ind]);
+ mpz_set (new_matrix -> storage[ind], matrix -> storage[(j-1)*old_rows+(i-1)]);
+ ind++;
+ }
+ }
+ }
+ }
+
+ return new_matrix;
+}
+
int
destroy_gmp_matrix(gmp_matrix * m)
{
diff --git a/contrib/kbipack/gmp_matrix.h b/contrib/kbipack/gmp_matrix.h
index 6e32583db36d840637cc64d42aac9236c449e3c4..784a87c20609f436bacda294851d6887107914b3 100755
--- a/contrib/kbipack/gmp_matrix.h
+++ b/contrib/kbipack/gmp_matrix.h
@@ -22,7 +22,7 @@
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 $
+ $Id: gmp_matrix.h,v 1.2 2009-04-03 11:06:13 matti Exp $
*/
#ifndef __GMP_MATRIX_H__
@@ -42,6 +42,9 @@ typedef struct
gmp_matrix *
create_gmp_matrix(size_t rows, size_t cols,
const mpz_t * elems);
+gmp_matrix *
+create_gmp_matrix_int(size_t rows, size_t cols,
+ const long int * elems);
gmp_matrix *
create_gmp_matrix_identity(size_t dim);
@@ -49,6 +52,12 @@ create_gmp_matrix_identity(size_t dim);
gmp_matrix *
create_gmp_matrix_zero(size_t rows, size_t cols);
+/* Copies a block of a matrix to another matrix. No resizing (yet). */
+gmp_matrix *
+copy_gmp_matrix(const gmp_matrix * matrix,
+ const size_t start_row, const size_t start_col,
+ const size_t end_row, const size_t end_col);
+
int
destroy_gmp_matrix(gmp_matrix *);