Skip to content
Snippets Groups Projects
Commit e7903fa3 authored by Matti Pellika's avatar Matti Pellika
Browse files

Forming ChainComplex class. Coreduction update. Added copy_gmp_matrix 

method to kbipack.
parent 9d00fec6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/CellComplex.cpp
+ 153
− 114
View file @ e7903fa3
......@@ -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
......
Geo/CellComplex.h
+ 54
− 7
View file @ e7903fa3
......@@ -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();
};
......
contrib/kbipack/gmp_matrix.c
+ 84
− 1
View file @ e7903fa3
......@@ -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)
{
......
contrib/kbipack/gmp_matrix.h
+ 10
− 1
View file @ e7903fa3
......@@ -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 *);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment