- Timestamp:
- Aug 2, 2016, 6:08:41 PM (8 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- fbaeb59f9025699b9055c1091c0deef1925dbb40
- Parents:
- 6a97c544535256c37da0f10bcb6821b592546c68
- git-author:
- Yue <ren@mathematik.uni-kl.de>2016-08-02 18:08:41+02:00
- git-committer:
- Yue <ren@mathematik.uni-kl.de>2016-08-05 18:10:20+02:00
- Location:
- gfanlib
- Files:
-
- 12 added
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
gfanlib/Makefile.am
r6a97c54 r5cea7a 9 9 AM_CPPFLAGS = -I${top_srcdir} -I${top_builddir} -DGMPRATIONAL 10 10 11 SOURCES = gfanlib_zcone.cpp gfanlib_symmetry.cpp gfanlib_symmetriccomplex.cpp gfanlib_polyhedralfan.cpp gfanlib_zfan.cpp gfanlib_polymakefile.cpp 12 libgfan_la_SOURCES = $(SOURCES)11 # forcefully enable exceptions 12 CXXFLAGS+= $(FEXCEPTIONSFRTTI_CXXFLAGS) 13 13 14 libgfan_includedir =$(includedir)/gfanlib 14 AM_CXXFLAGS = -std=c++11 15 15 16 libgfan_include_HEADERS = gfanlib_z.h gfanlib_q.h gfanlib_vector.h gfanlib_matrix.h gfanlib_zcone.h gfanlib.h gfanlib_polyhedralfan.h gfanlib_polymakefile.h gfanlib_symmetriccomplex.h gfanlib_symmetry.h gfanlib_zfan.h 16 17 SOURCES = gfanlib_circuittableint.cpp gfanlib_mixedvolume.cpp gfanlib_paralleltraverser.cpp gfanlib_polyhedralfan.cpp gfanlib_polymakefile.cpp gfanlib_symmetriccomplex.cpp gfanlib_symmetry.cpp gfanlib_traversal.cpp gfanlib_zcone.cpp gfanlib_zfan.cpp 18 libgfan_la_SOURCES = $(SOURCES) 19 20 libgfan_includedir =$(includedir)/gfanlib 21 libgfan_include_HEADERS = config.h gfanlib_mixedvolume.h gfanlib_polymakefile.h gfanlib_symmetry.h gfanlib_vector.h gfanlib_z.h _config.h gfanlib.h gfanlib_paralleltraverser.h gfanlib_q.h gfanlib_traversal.h gfanlib_zcone.h gfanlib_circuittableint.h gfanlib_matrix.h gfanlib_polyhedralfan.h gfanlib_symmetriccomplex.h gfanlib_tropicalhomotopy.h gfanlib_zfan.h 17 22 18 23 DISTCLEANFILES = config.h -
gfanlib/gfanlib.h
r6a97c54 r5cea7a 16 16 #include "gfanlib_polyhedralfan.h" 17 17 #include "gfanlib_zfan.h" 18 #include "gfanlib_mixedvolume.h" 18 19 19 20 #endif /* GFANLIB_H_ */ -
gfanlib/gfanlib_matrix.h
r6a97c54 r5cea7a 9 9 #define LIB_ZMATRIX_H_ 10 10 11 #include <sstream>12 11 #include <vector> 13 12 #include <algorithm> … … 18 17 template <class typ> class Matrix{ 19 18 int width,height; 20 std::vector<Vector<typ> > rows; 19 // std::vector<Vector<typ> > rows; 20 std::vector<typ> data; 21 21 public: 22 22 // rowIterator; … … 25 25 inline int getHeight()const{return height;}; 26 26 inline int getWidth()const{return width;}; 27 Matrix(const Matrix &a):width(a.getWidth()),height(a.getHeight()), rows(a.rows){28 } 29 Matrix(int height_, int width_):width(width_),height(height_), rows(height_){27 Matrix(const Matrix &a):width(a.getWidth()),height(a.getHeight()),data(a.data){ 28 } 29 Matrix(int height_, int width_):width(width_),height(height_),data(width_*height_){ 30 30 assert(height>=0); 31 31 assert(width>=0); 32 for(int i=0;i<getHeight();i++)rows[i]=Vector<typ>(width);33 32 }; 34 33 ~Matrix(){ … … 36 35 Matrix():width(0),height(0){ 37 36 }; 38 Matrix rowVectorMatrix(Vector<typ> const &v)37 static Matrix rowVectorMatrix(Vector<typ> const &v) 39 38 { 40 39 Matrix ret(1,v.size()); 41 for( unsignedi=0;i<v.size();i++)ret[0][i]=v[i];40 for(int i=0;i<v.size();i++)ret[0][i]=v[i]; 42 41 return ret; 43 42 } … … 47 46 assert(i<getWidth()); 48 47 Vector<typ> ret(getHeight()); 49 for(int j=0;j<getHeight();j++)ret[j]= rows[j][i];48 for(int j=0;j<getHeight();j++)ret[j]=(*this)[j][i]; 50 49 return ret; 51 50 } … … 54 53 Matrix ret(getWidth(),getHeight()); 55 54 for(int i=0;i<getWidth();i++) 56 ret.rows[i]=column(i); 55 for(int j=0;j<getHeight();j++) 56 ret[i][j]=(*this)[j][i]; 57 57 return ret; 58 58 } … … 60 60 { 61 61 Matrix m(n,n); 62 for(int i=0;i<n;i++)m .rows[i]=Vector<typ>::standardVector(n,i);62 for(int i=0;i<n;i++)m[i][i]=typ(1); 63 63 return m; 64 64 } 65 65 void append(Matrix const &m) 66 66 { 67 assert(m.getWidth()==width); 68 data.resize((height+m.height)*width); 69 int oldHeight=height; 70 height+=m.height; 67 71 for(int i=0;i<m.height;i++) 68 72 { 69 rows.push_back(m[i]); 73 for(int j=0;j<m.width;j++) 74 (*this)[i+oldHeight][j]=m[i][j]; 70 75 } 71 height+=m.height;72 76 } 73 77 void appendRow(Vector<typ> const &v) 74 { 75 assert((int)v.size()==width); 76 rows.push_back(v); 77 height++; 78 } 79 void prependRow(Vector<typ> const &v) 80 { 81 assert((int)v.size()==width); 82 rows.insert(rows.begin(),v); 83 height++; 84 } 78 { 79 assert(v.size()==width); 80 data.resize((height+1)*width); 81 height++; 82 for(int j=0;j<width;j++) 83 (*this)[height-1][j]=v[j]; 84 } 85 85 void eraseLastRow() 86 86 { 87 assert( rows.size()>0);88 rows.pop_back();87 assert(height>0); 88 data.resize((height-1)*width); 89 89 height--; 90 90 } … … 111 111 { 112 112 Matrix p=q; 113 for(int i=0;i<q.height;i++)p[i]=s*(q[i]); 113 for(int i=0;i<q.height;i++) 114 for(int j=0;j<q.width;j++) 115 p[i][j]=s*(q[i][j]); 114 116 return p; 115 117 } … … 147 149 for(int i=startRow;i<endRow;i++) 148 150 for(int j=startColumn;j<endColumn;j++) 149 ret[i-startRow][j-startColumn]=rows[i][j]; 150 return ret; 151 } 152 const Vector<typ>& operator[](int n)const{assert(n>=0 && n<getHeight());return (rows[n]);} 153 // Bugfix for gcc4.5 (passing assertion to the above operator): 154 Vector<typ>& operator[](int n){if(!(n>=0 && n<getHeight())){(*(const Matrix<typ>*)(this))[n];}return (rows[n]);} 151 ret[i-startRow][j-startColumn]=(*this)[i][j]; 152 return ret; 153 } 154 155 class RowRef; 156 class const_RowRef{ 157 int rowNumM; 158 Matrix const &matrix; 159 friend class RowRef; 160 public: 161 inline const_RowRef(const Matrix &matrix_, int rowNum_)__attribute__((always_inline)): 162 rowNumM(rowNum_*matrix_.width), 163 matrix(matrix_) 164 { 165 } 166 inline typ const &operator[](int j)const __attribute__((always_inline)) 167 { 168 assert(j>=0); 169 assert(j<matrix.width); 170 return matrix.data[rowNumM+j]; 171 } 172 inline typ const &UNCHECKEDACCESS(int j)const __attribute__((always_inline)) 173 { 174 return matrix.data[rowNumM+j]; 175 } 176 const Vector<typ> toVector()const 177 { 178 Vector<typ> ret(matrix.width); 179 for(int j=0;j<matrix.width;j++) 180 ret[j]=matrix.data[rowNumM+j]; 181 return ret; 182 } 183 operator Vector<typ>()const 184 { 185 return toVector(); 186 } 187 bool operator==(Vector<typ> const &b)const 188 { 189 return toVector()==b; 190 } 191 /* typ dot(Vector<typ> const &b)const 192 { 193 return dot(toVector(),b); 194 }*/ 195 Vector<typ> operator-()const 196 { 197 return -toVector(); 198 } 199 }; 200 class RowRef{ 201 int rowNumM; 202 Matrix &matrix; 203 public: 204 inline RowRef(Matrix &matrix_, int rowNum_): 205 rowNumM(rowNum_*matrix_.width), 206 matrix(matrix_) 207 { 208 } 209 inline typ &operator[](int j) __attribute__((always_inline)) 210 { 211 assert(j>=0); 212 assert(j<matrix.width); 213 return matrix.data[rowNumM+j]; 214 } 215 inline typ &UNCHECKEDACCESS(int j) 216 { 217 return matrix.data[rowNumM+j]; 218 } 219 RowRef &operator=(Vector<typ> const &v) 220 { 221 assert(v.size()==matrix.width); 222 for(int j=0;j<matrix.width;j++) 223 matrix.data[rowNumM+j]=v[j]; 224 225 return *this; 226 } 227 RowRef &operator=(RowRef const &v) 228 { 229 assert(v.matrix.width==matrix.width); 230 for(int j=0;j<matrix.width;j++) 231 matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j]; 232 233 return *this; 234 } 235 /* RowRef &operator+=(Vector<typ> const &v) 236 { 237 assert(v.size()==matrix.width); 238 for(int j=0;j<matrix.width;j++) 239 matrix.data[rowNumM+j]+=v.v[j]; 240 241 return *this; 242 }*/ 243 RowRef &operator+=(RowRef const &v) 244 { 245 assert(v.matrix.width==matrix.width); 246 for(int j=0;j<matrix.width;j++) 247 matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j]; 248 249 return *this; 250 } 251 RowRef &operator+=(const_RowRef const &v) 252 { 253 assert(v.matrix.width==matrix.width); 254 for(int j=0;j<matrix.width;j++) 255 matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j]; 256 257 return *this; 258 } 259 RowRef &operator=(const_RowRef const &v) 260 { 261 assert(v.matrix.width==matrix.width); 262 for(int j=0;j<matrix.width;j++) 263 matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j]; 264 265 return *this; 266 } 267 const Vector<typ> toVector()const 268 { 269 Vector<typ> ret(matrix.width); 270 for(int j=0;j<matrix.width;j++) 271 ret[j]=matrix.data[rowNumM+j]; 272 return ret; 273 } 274 operator Vector<typ>()const 275 { 276 return toVector(); 277 } 278 /* typ dot(Vector<typ> const &b)const 279 { 280 return dot(toVector(),b); 281 }*/ 282 bool isZero()const 283 { 284 for(int j=0;j<matrix.width;j++)if(!(matrix.data[rowNumM+j].isZero()))return false; 285 return true; 286 } 287 }; 288 289 290 inline RowRef operator[](int i) __attribute__((always_inline)) 291 { 292 assert(i>=0); 293 assert(i<height); 294 return RowRef(*this,i); 295 } 296 inline const_RowRef operator[](int i)const __attribute__((always_inline)) 297 { 298 assert(i>=0); 299 assert(i<height); 300 return const_RowRef(*this,i); 301 } 302 303 304 const typ& UNCHECKEDACCESS(int i,int j)const __attribute__((always_inline)){ 305 /* assert(i>=0); 306 assert(i<height); 307 assert(j>=0); 308 assert(j<width);*/ 309 return data[i*width+j];} 310 typ& UNCHECKEDACCESS(int i,int j) __attribute__((always_inline)){ 311 /* assert(i>=0); 312 assert(i<height); 313 assert(j>=0); 314 assert(j<width);*/ 315 return data[i*width+j];} 316 155 317 156 318 … … 164 326 for(int i=0;i<getHeight();i++) 165 327 { 166 if((*this)[i] <b[i])return true;167 if(b[i] <(*this)[i])return false;328 if((*this)[i].toVector()<b[i].toVector())return true; 329 if(b[i].toVector()<(*this)[i].toVector())return false; 168 330 } 169 331 return false; … … 180 342 if(!a.isZero()) 181 343 for(int k=0;k<width;k++) 182 if(! rows[i][k].isZero())183 rows[j][k].madd(rows[i][k],a);344 if(!(*this)[i][k].isZero()) 345 (*this)[j][k].madd((*this)[i][k],a); 184 346 } 185 347 … … 189 351 { 190 352 if(i)f<<","<<std::endl; 191 f<<a .rows[i];353 f<<a[i].toVector(); 192 354 } 193 355 f<<"}"<<std::endl; … … 197 359 std::string toString()const 198 360 { 199 200 201 361 std::stringstream f; 362 f<<*this; 363 return f.str(); 202 364 } 203 365 … … 207 369 void swapRows(int i, int j) 208 370 { 209 std::swap(rows[i],rows[j]);371 for(int a=0;a<width;a++)std::swap((*this)[i][a],(*this)[j][a]); 210 372 } 211 373 /** … … 222 384 while(++j<width) 223 385 { 224 if(! rows[i][j].isZero()) return true;386 if(!(*this)[i][j].isZero()) return true; 225 387 } 226 388 return false; … … 295 457 int bestNumberOfNonZero=0; 296 458 for(int i=currentRow;i<height;i++) 297 if(! rows[i][column].isZero())459 if(!(*this)[i][column].isZero()) 298 460 { 299 461 int nz=0; 300 462 for(int k=column+1;k<width;k++) 301 if(! rows[i][k].isZero())nz++;463 if(!(*this)[i][k].isZero())nz++; 302 464 if(best==-1) 303 465 { … … 345 507 if(makePivotsOne) 346 508 {//THE PIVOT SHOULD BE SET TO ONE IF INTEGRAL IS FALSE 347 if( rows[currentRow][i].sign()>=0)retSwaps++;348 typ inverse=typ(1)/ rows[currentRow][i];509 if((*this)[currentRow][i].sign()>=0)retSwaps++; 510 typ inverse=typ(1)/(*this)[currentRow][i]; 349 511 // if(!rows[currentRow][i].isOne()) 350 512 { 351 513 for(int k=0;k<width;k++) 352 if(! rows[currentRow][k].isZero())353 rows[currentRow][k]*=inverse;514 if(!(*this)[currentRow][k].isZero()) 515 (*this)[currentRow][k]*=inverse; 354 516 } 355 517 } … … 357 519 if(integral) 358 520 { 359 if(! rows[j][i].isZero())521 if(!(*this)[j][i].isZero()) 360 522 { 361 523 typ s; 362 524 typ t; 363 525 364 typ g=typ::gcd( rows[currentRow][i],rows[j][i],s,t);365 typ u=- rows[j][i]/g;366 typ v= rows[currentRow][i]/g;526 typ g=typ::gcd((*this)[currentRow][i],(*this)[j][i],s,t); 527 typ u=-(*this)[j][i]/g; 528 typ v=(*this)[currentRow][i]/g; 367 529 /* We want the (s,t) vector to be as small as possible. 368 530 * We are allowed to adjust by multiples of (u,v). … … 378 540 for(int k=0;k<width;k++) 379 541 { 380 typ A= rows[currentRow][k];381 typ B= rows[j][k];382 383 rows[currentRow][k]=s*A+t*B;384 rows[j][k]=u*A+v*B;542 typ A=(*this)[currentRow][k]; 543 typ B=(*this)[j][k]; 544 545 (*this)[currentRow][k]=s*A+t*B; 546 (*this)[j][k]=u*A+v*B; 385 547 } 386 548 } … … 388 550 else 389 551 { 390 if(! rows[j][i].isZero())391 madd(currentRow,- rows[j][i]/rows[currentRow][i],j);552 if(!(*this)[j][i].isZero()) 553 madd(currentRow,-(*this)[j][i]/(*this)[currentRow][i],j); 392 554 } 393 555 currentRow++; … … 412 574 while(nextPivot(pivotI,pivotJ)) 413 575 { 414 415 rows[pivotI]=rows[pivotI]/rows[pivotI][pivotJ];576 if(scalePivotsToOne) 577 (*this)[pivotI]=(*this)[pivotI].toVector()/(*this)[pivotI][pivotJ]; 416 578 for(int i=0;i<pivotI;i++) 417 if(! rows[i][pivotJ].isZero())418 madd(pivotI,- rows[i][pivotJ]/rows[pivotI][pivotJ],i);419 }579 if(!(*this)[i][pivotJ].isZero()) 580 madd(pivotI,-(*this)[i][pivotJ]/(*this)[pivotI][pivotJ],i); 581 } 420 582 return ret; 421 583 } … … 438 600 if(!v[pivotJ].isZero()) 439 601 { 440 typ s=-v[pivotJ]/ rows[pivotI][pivotJ];602 typ s=-v[pivotJ]/(*this)[pivotI][pivotJ]; 441 603 442 604 for(int k=0;k<width;k++) 443 if(! rows[pivotI][k].isZero())444 v[k].madd( rows[pivotI][k],s);605 if(!(*this)[pivotI][k].isZero()) 606 v[k].madd((*this)[pivotI][k],s); 445 607 } 446 608 return v; … … 472 634 while(nextPivot(pivot2I,pivot2J)) 473 635 { 474 ret[k][pivot2J]= rows[pivot2I][j]/rows[pivot2I][pivot2J];636 ret[k][pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J]; 475 637 } 476 638 ret[k][j]=typ(-1); … … 508 670 while(nextPivot(pivot2I,pivot2J)) 509 671 { 510 diagonalProduct*= rows[pivot2I][pivot2J];672 diagonalProduct*=(*this)[pivot2I][pivot2J]; 511 673 } 512 674 } … … 522 684 while(nextPivot(pivot2I,pivot2J)) 523 685 { 524 ret[pivot2J]= rows[pivot2I][j]/rows[pivot2I][pivot2J];686 ret[pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J]; 525 687 lastEntry-=ret[pivot2J]*ret[pivot2J]; 526 688 } … … 546 708 * Sort the rows of the matrix. 547 709 */ 710 struct rowComparer{ 711 bool operator()(std::pair<Matrix*,int> i, std::pair<Matrix*,int> j) {return ((*i.first)[i.second].toVector()<(*j.first)[j.second].toVector());} 712 } theRowComparer; 548 713 void sortRows() 549 714 { 550 std::sort(rows.begin(),rows.end()); 715 std::vector<std::pair<Matrix*,int> > v; 716 for(int i=0;i<height;i++)v.push_back(std::pair<Matrix*,int>(this,i)); 717 std::sort(v.begin(),v.end(),theRowComparer); 718 Matrix result(height,width); 719 for(int i=0;i<height;i++) 720 result[i]=(*this)[v[i].second].toVector(); 721 data=result.data; 551 722 } 552 723 /** … … 560 731 B.appendRow((*this)[0]); 561 732 for(int i=1;i<getHeight();i++) 562 if( rows[i]!=rows[i-1])B.appendRow((*this)[i]);733 if((*this)[i].toVector()!=(*this)[i-1].toVector())B.appendRow((*this)[i].toVector()); 563 734 *this=B; 564 735 } … … 573 744 assert(bottom.getWidth()==top.getWidth()); 574 745 Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth()); 575 for(int i=0;i<top.getHeight();i++)ret .rows[i]=top.rows[i];576 for(int i=0;i<bottom.getHeight();i++)ret .rows[i+top.getHeight()]=bottom.rows[i];746 for(int i=0;i<top.getHeight();i++)ret[i]=top[i]; 747 for(int i=0;i<bottom.getHeight();i++)ret[i+top.getHeight()]=bottom[i]; 577 748 578 749 return ret; … … 590 761 for(int i=0;i<left.getHeight();i++) 591 762 { 592 for(int j=0;j<left.getWidth();j++)ret .rows[i][j]=left.rows[i][j];593 for(int j=0;j<right.getWidth();j++)ret .rows[i][j+left.getWidth()]=right.rows[i][j];763 for(int j=0;j<left.getWidth();j++)ret[i][j]=left[i][j]; 764 for(int j=0;j<right.getWidth();j++)ret[i][j+left.getWidth()]=right[i][j]; 594 765 } 595 766 return ret; -
gfanlib/gfanlib_polyhedralfan.cpp
r6a97c54 r5cea7a 71 71 int PolyhedralFan::getMaxDimension()const 72 72 { 73 if (cones.empty()) 74 return -1; 73 assert(!cones.empty()); 75 74 76 75 return cones.begin()->dimension(); … … 79 78 int PolyhedralFan::getMinDimension()const 80 79 { 81 if (cones.empty()) 82 return -1; 80 assert(!cones.empty()); 83 81 84 82 return cones.rbegin()->dimension(); … … 355 353 bool notAll=false; 356 354 for(unsigned j=0;j<theCone.indices.size();j++) 357 if(dot(rays[theCone.indices[j]] ,facetCandidates[i]).sign()==0)355 if(dot(rays[theCone.indices[j]].toVector(),facetCandidates[i].toVector()).sign()==0) 358 356 indices.insert(theCone.indices[j]); 359 357 else … … 864 862 int PolyhedralFan::dimensionOfLinealitySpace()const 865 863 { 866 if (cones.empty()) 867 return n; 864 assert(cones.size());//slow! 868 865 return cones.begin()->dimensionOfLinealitySpace(); 869 866 } -
gfanlib/gfanlib_q.h
r6a97c54 r5cea7a 11 11 #include <string.h> 12 12 #include <ostream> 13 #include <assert.h> 13 14 #include "gmp.h" 14 15 … … 160 161 return mpq_sgn(value); 161 162 } 162 static Rational gcd(Rational const &a, Rational const /*&b*/, Rational /*&s*/, Rational /*&t*/)163 static Rational gcd(Rational const &a, Rational const &b, Rational &s, Rational &t) 163 164 { 164 165 /* mpz_t r; … … 167 168 Integer ret(r); 168 169 mpz_clear(r);*/ 169 assert(0); 170 return a; 171 } 172 static Rational gcd(Rational const &a, Rational const /*&b*/) 173 { 174 /* mpz_t r; 175 mpz_init(r); 176 mpz_gcdext(r,s.value,t.value,a.value,b.value); 177 Integer ret(r); 178 mpz_clear(r);*/ 179 assert(0); 170 assert(0 && "gcd for Rational not defined"); 180 171 return a; 181 172 } -
gfanlib/gfanlib_symmetriccomplex.cpp
r6a97c54 r5cea7a 185 185 int SymmetricComplex::getLinDim()const 186 186 { 187 187 return linealitySpace.getHeight(); 188 188 } 189 189 … … 595 595 596 596 if(flags&FPF_cones)polymakeFile.writeStringProperty("CONES",toStringJustCones(getMinDim(),getMaxDim(),false,flags&FPF_group, 0,false,flags&FPF_tPlaneSort)); 597 std::stringstream multiplicities; 598 if(flags&FPF_maximalCones)polymakeFile.writeStringProperty("MAXIMAL_CONES",toStringJustCones(getMinDim(),getMaxDim(),true,flags&FPF_group, &multiplicities,false,flags&FPF_tPlaneSort)); 597 if(flags&FPF_maximalCones)polymakeFile.writeStringProperty("MAXIMAL_CONES",toStringJustCones(getMinDim(),getMaxDim(),true,flags&FPF_group, 0,false,flags&FPF_tPlaneSort)); 599 598 if(flags&FPF_conesCompressed)polymakeFile.writeStringProperty("CONES_ORBITS",toStringJustCones(getMinDim(),getMaxDim(),false,flags&FPF_group, 0,true,flags&FPF_tPlaneSort)); 600 599 if((flags&FPF_conesCompressed) && (flags&FPF_maximalCones))polymakeFile.writeStringProperty("MAXIMAL_CONES_ORBITS",toStringJustCones(getMinDim(),getMaxDim(),true,flags&FPF_group, 0,true,flags&FPF_tPlaneSort)); 601 if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES",multiplicities.str());602 600 603 601 if(!sym.isTrivial()) … … 637 635 stringstream multiplicities; 638 636 polymakeFile.writeStringProperty("MAXIMAL_CONES",symCom.toString(symCom.getMinDim(),symCom.getMaxDim(),true,flags&FPF_group, &multiplicities,false,flags&FPF_tPlaneSort)); 637 if(flags&FPF_multiplicities)polymakeFile.writeStringProperty("MULTIPLICITIES",multiplicities.str()); 639 638 // log1 fprintf(Stderr,"Done producing list of maximal cones.\n"); 640 639 } -
gfanlib/gfanlib_vector.h
r6a97c54 r5cea7a 72 72 // Access 73 73 //-------- 74 typ& operator[](int n) 74 typ& operator[](int n)__attribute__((always_inline)) 75 75 { 76 76 if(!(n>=0 && n<(int)v.size()))outOfRange(n,v.size()); 77 77 return (v[n]); 78 78 } 79 const typ& operator[](int n)const{assert(n>=0 && n<(int)v.size());return (v[n]);} 80 const typ& UNCHECKEDACCESS(int n)const{return (v[n]);} 79 const typ& operator[](int n)const __attribute__((always_inline)){assert(n>=0 && n<(int)v.size());return (v[n]);} 80 typ& UNCHECKEDACCESS(int n)__attribute__((always_inline)){return (v[n]);} 81 const typ& UNCHECKEDACCESS(int n)const __attribute__((always_inline)){return (v[n]);} 81 82 82 83 //------------- … … 85 86 unsigned int size()const{return v.size();}; 86 87 void resize(int n){v.resize(n,typ());}; 87 void grow(int i){if( (int)size()<i)resize(i);}88 void grow(int i){if(size()<i)resize(i);} 88 89 void push_back(typ a) 89 90 { … … 293 294 std::string toString()const 294 295 { 295 296 297 296 std::stringstream f; 297 f<<*this; 298 return f.str(); 298 299 } 299 300 300 301 typ gcd()const 301 302 { 303 typ temp1,temp2; 302 304 typ ret(1); 303 305 for(unsigned i=0;i<size();i++) 304 ret=typ::gcd(ret,v[i] );306 ret=typ::gcd(ret,v[i],temp1,temp2); 305 307 return ret; 306 308 } … … 309 311 assert(!typ::isField()); 310 312 return (*this)/gcd(); 311 }312 313 void debugPrint()const314 {315 std::stringstream s;316 s<<"(";317 for(typename std::vector<typ>::const_iterator i=this->v.begin();i!=this->v.end();i++)318 {319 if(i!=this->v.begin()) s<<",";320 s<<*i;321 }322 s<<")"<<std::endl;323 std::cout << s.str();324 return;325 313 } 326 314 }; … … 415 403 416 404 return ret; 417 } ;418 419 } ;405 } 406 407 } 420 408 421 409 #endif /* LIB_ZVECTOR_H_ */ -
gfanlib/gfanlib_z.h
r6a97c54 r5cea7a 10 10 11 11 #include <string.h> 12 #include <iostream>13 12 #include <ostream> 14 13 … … 41 40 mpz_init_set(value,value_.value); 42 41 } 43 Integer(mpz_t value_)42 Integer(mpz_t const value_) 44 43 { 45 44 mpz_init_set(value,value_); … … 60 59 bool isZero()const{ 61 60 return mpz_sgn(value)==0; 62 }63 bool isOne()const{64 return mpz_cmp_si(value,1);65 61 } 66 62 friend std::ostream &operator<<(std::ostream &f, Integer const &a) … … 73 69 return f; 74 70 } 75 void debugPrint() const76 {77 void (*freefunc)(void *, size_t);78 mp_get_memory_functions(0,0,&freefunc);79 char *str=mpz_get_str(0,10,value);80 std::cout << str;81 freefunc(str,strlen(str)+1);82 return;83 }84 71 Integer& operator+=(const Integer& a) 85 72 { … … 161 148 mpz_init(r); 162 149 mpz_gcdext(r,s.value,t.value,a.value,b.value); 163 Integer ret(r);164 mpz_clear(r);165 return ret;166 }167 static Integer gcd(Integer const &a, Integer const &b)168 {169 if (a.isOne() || b.isOne())170 return Integer(1);171 mpz_t r;172 mpz_init(r);173 mpz_gcd(r,a.value,b.value);174 150 Integer ret(r); 175 151 mpz_clear(r); … … 354 330 return f; 355 331 } 356 friend void debugPrint(IntegerTemplate const &a)357 {358 std::cout << a << std::endl;359 return;360 }361 332 LimbWord signExtension(LimbWord a) 362 333 { -
gfanlib/gfanlib_zcone.cpp
r6a97c54 r5cea7a 12 12 #include <sstream> 13 13 14 #include <config.h> 15 #ifdef HAVE_CDD_SETOPER_H16 #include <cdd/setoper.h>17 #include <cdd/cdd.h>14 //extern "C"{ 15 #ifdef NOCDDPREFIX 16 #include "setoper.h" 17 #include "cdd.h" 18 18 #else 19 #ifdef HAVE_CDDLIB_SETOPER_H 20 #include <cddlib/setoper.h> 21 #include <cddlib/cdd.h> 22 #else 23 #include <setoper.h> 24 #include <cdd.h> 25 #endif //HAVE_CDDLIB_SETOPER_H 26 #endif //HAVE_CDD_SETOPER_H 19 #include "cdd/setoper.h" 20 #include "cdd/cdd.h" 21 #endif 22 //} 27 23 28 24 namespace gfan{ 29 30 static void cddinitGmp() 31 { 32 static bool initialized; 33 if(!initialized) 25 bool isCddlibRequired() 26 { 27 return true; 28 } 29 void initializeCddlibIfRequired() // calling this frequently will cause memory leaks because deinitialisation is not possible with old versions of cddlib. 30 { 31 dd_set_global_constants(); 32 } 33 void deinitializeCddlibIfRequired() 34 { 35 // dd_free_global_constants(); 36 } 37 static void ensureCddInitialisation() 38 { 39 // A more complicated initialisation than the following (meaning attempts to count the number of times 40 // cddlib was requested to be initialised) would require cddlib to be thread aware. 41 // The error below is implemented with an assert(0) because throwing an exception may leave the impression that 42 // it is possible to recover from this error. While that may be true, it would not work in full generality, 43 // as the following if statement cannot test whether dd_free_global_constants() has also been called. 44 // Moverover, in multithreaded environments it would be quite difficult to decide if cddlib was initialised. 45 if(!dd_one[0]._mp_num._mp_d) 46 { 47 std::cerr<<"CDDLIB HAS NOT BEEN INITIALISED!\n" 48 "\n" 49 "Fix this problem by calling the following function in your initialisation code:\n" 50 "dd_set_global_constants();\n" 51 "(after possibly setting the gmp allocators) and\n" 52 "dd_free_global_constants()\n" 53 "in your deinitialisation code (only available for cddlib version>=094d).\n" 54 "This requires the header includes:\n" 55 "#include \"cdd/setoper.h\"\n" 56 "#include \"cdd/cdd.h\"\n" 57 "\n" 58 "Alternatively, you may call gfan:initializeCddlibIfRequired() and deinitializeCddlibIfRequired()\n" 59 "if gfanlib is the only code using cddlib. If at some point cddlib is no longer required by gfanlib\n" 60 "these functions may do nothing.\n" 61 "Because deinitialisation is not possible in cddlib <094d, the functions may leak memory and should not be called often.\n" 62 "\n" 63 "This error message will never appear if the initialisation was done properly, and therefore never appear in a shipping version of your software.\n"; 64 assert(0); 65 } 66 /* 67 static bool initialized; 68 if(!initialized) 34 69 { 35 dd_set_global_constants(); /* First, this must be called. */ 36 37 } 70 dd_set_global_constants(); 71 initialized=true; 72 }*/ 38 73 } 39 74 … … 85 120 int numberOfRows=g.getHeight(); 86 121 dd_MatrixPtr A=NULL; 87 cddinitGmp();122 ensureCddInitialisation(); 88 123 A=ZMatrix2MatrixGmp(g, err); 89 124 for(int i=numberOfInequalities;i<numberOfRows;i++) … … 121 156 FILE *reading=NULL; 122 157 123 cddinitGmp();158 ensureCddInitialisation(); 124 159 125 160 M=ZMatrix2MatrixGmp(g, &err); … … 128 163 // d=M->colsize; 129 164 130 static dd_Arow temp; 165 //static dd_Arow temp; 166 dd_Arow temp; 131 167 dd_InitializeArow(g.getWidth()+1,&temp); 132 168 … … 268 304 for(int i=0;i<a.getHeight();i++) 269 305 { 270 assert(!(a[i]. isZero()));271 b.insert(a[i]. normalized());306 assert(!(a[i].toVector().isZero())); 307 b.insert(a[i].toVector().normalized()); 272 308 } 273 309 … … 302 338 for(int i=0;i!=original.getHeight();i++) 303 339 for(int j=0;j!=a.getHeight();j++) 304 if(!dependent(original[i] ,a[j]))340 if(!dependent(original[i].toVector(),a[j].toVector())) 305 341 { 306 342 ZVector const &I=original[i]; … … 337 373 void removeRedundantRows(ZMatrix &inequalities, ZMatrix &equations, bool removeInequalityRedundancies) 338 374 { 339 cddinitGmp();375 ensureCddInitialisation(); 340 376 341 377 int numberOfEqualities=equations.getHeight(); … … 403 439 { 404 440 QVector retUnscaled(inequalities.getWidth()); 405 cddinitGmp();441 ensureCddInitialisation(); 406 442 int numberOfEqualities=equations.getHeight(); 407 443 int numberOfInequalities=inequalities.getHeight(); … … 459 495 dd_ErrorType err=dd_NoError; 460 496 461 cddinitGmp();497 ensureCddInitialisation(); 462 498 463 499 A=ZMatrix2MatrixGmp(inequalities, equations, &err); … … 567 603 dd_ErrorType err=dd_NoError; 568 604 569 cddinitGmp();605 ensureCddInitialisation(); 570 606 A=ZMatrix2MatrixGmp(inequalities, &err); 571 607 … … 748 784 std::string ZCone::toString()const 749 785 { 750 std::stringstream f; 751 f<<*this; 752 return f.str(); 753 // ======= 754 // std::stringstream s; 755 // s<<"AMBIENT_DIM"<<std::endl; 756 // s<<this->ambientDimension()<<std::endl; 757 758 // gfan::ZMatrix i=this->getInequalities(); 759 // if (this->areFacetsKnown()) 760 // s<<"FACETS"<<std::endl; 761 // else 762 // s<<"INEQUALITIES"<<std::endl; 763 // s<<i<<std::endl; 764 765 // gfan::ZMatrix e=this->getEquations(); 766 // if (this->areImpliedEquationsKnown()) 767 // s<<"LINEAR_SPAN"<<std::endl; 768 // else 769 // s<<"EQUATIONS"<<std::endl; 770 // s<<e<<std::endl; 771 772 // gfan::ZMatrix r=this->extremeRays(); 773 // s<<"RAYS"<<std::endl; 774 // s<<r<<std::endl; 775 776 // gfan::ZMatrix l=this->generatorsOfLinealitySpace(); 777 // s<<"LINEALITY_SPACE"<<std::endl; 778 // s<<l<<std::endl; 779 780 // std::cout << s.str(); 781 // return; 782 // >>>>>>> chg: status update 17.07. 786 std::stringstream f; 787 f<<*this; 788 return f.str(); 783 789 } 784 790 … … 983 989 } 984 990 */ 985 // 986 991 // ZCone dual(newGenerators,linealitySpace); 992 ZCone dual(generators,linealitySpace); 987 993 // dual.findFacets(); 988 994 // dual.canonicalize(); … … 1238 1244 1239 1245 for(int i=0;i<M.getHeight();i++) 1240 if(M[i]. subvector(0,a).isZero()&&!M[i].subvector(a,a+b).isZero())1246 if(M[i].toVector().subvector(0,a).isZero()&&!M[i].toVector().subvector(a,a+b).isZero()) 1241 1247 { 1242 ret.appendRow(M[i]. subvector(a+b,a+b+n));1248 ret.appendRow(M[i].toVector().subvector(a+b,a+b+n)); 1243 1249 } 1244 1250 … … 1252 1258 assert(temp.getHeight()==1); 1253 1259 for(int i=0;i<inequalities.getHeight();i++) 1254 if(dot(temp[0] ,inequalities[i]).sign()<0)1260 if(dot(temp[0].toVector(),inequalities[i].toVector()).sign()<0) 1255 1261 { 1256 temp[0]=-temp[0] ;1262 temp[0]=-temp[0].toVector(); 1257 1263 break; 1258 1264 } … … 1338 1344 } 1339 1345 1340 } ;1346 } -
gfanlib/gfanlib_zcone.h
r6a97c54 r5cea7a 10 10 11 11 #include "gfanlib_matrix.h" 12 #include <sstream>13 12 14 13 namespace gfan{ 15 14 /** 15 * Returns true if cddlib is needed for the ZCone implementation. 16 */ 17 bool isCddlibRequired(); 18 /** 19 * Only call this function if gfanlib is the only code in your program using cddlib. 20 * Should be paired with a deinitializeCddlibIfRequired() call. 21 * Calling the function repeatedly may cause memory leaks even if deinitializeCddlibIfRequired() is also called. 22 */ 23 void initializeCddlibIfRequired(); 24 /** 25 * This function may do nothing. 26 */ 27 void deinitializeCddlibIfRequired(); 16 28 17 29 /** … … 358 370 }; 359 371 360 }; 372 } 373 374 361 375 362 376
Note: See TracChangeset
for help on using the changeset viewer.