Changeset 5cea7a in git for gfanlib


Ignore:
Timestamp:
Aug 2, 2016, 6:08:41 PM (8 years ago)
Author:
Yue <ren@…>
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
Message:
chg: new gfanlib version 0.6
Location:
gfanlib
Files:
12 added
10 edited

Legend:

Unmodified
Added
Removed
  • gfanlib/Makefile.am

    r6a97c54 r5cea7a  
    99AM_CPPFLAGS = -I${top_srcdir} -I${top_builddir} -DGMPRATIONAL
    1010
    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
     12CXXFLAGS+= $(FEXCEPTIONSFRTTI_CXXFLAGS)
    1313
    14 libgfan_includedir  =$(includedir)/gfanlib
     14AM_CXXFLAGS = -std=c++11
    1515
    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
     17SOURCES = 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
     18libgfan_la_SOURCES = $(SOURCES)
     19
     20libgfan_includedir =$(includedir)/gfanlib
     21libgfan_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
    1722
    1823DISTCLEANFILES =  config.h
  • gfanlib/gfanlib.h

    r6a97c54 r5cea7a  
    1616#include "gfanlib_polyhedralfan.h"
    1717#include "gfanlib_zfan.h"
     18#include "gfanlib_mixedvolume.h"
    1819
    1920#endif /* GFANLIB_H_ */
  • gfanlib/gfanlib_matrix.h

    r6a97c54 r5cea7a  
    99#define LIB_ZMATRIX_H_
    1010
    11 #include <sstream>
    1211#include <vector>
    1312#include <algorithm>
     
    1817template <class typ> class Matrix{
    1918  int width,height;
    20   std::vector<Vector<typ> > rows;
     19//  std::vector<Vector<typ> > rows;
     20  std::vector<typ> data;
    2121public:
    2222  // rowIterator;
     
    2525  inline int getHeight()const{return height;};
    2626  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_){
    3030    assert(height>=0);
    3131    assert(width>=0);
    32     for(int i=0;i<getHeight();i++)rows[i]=Vector<typ>(width);
    3332  };
    3433  ~Matrix(){
     
    3635  Matrix():width(0),height(0){
    3736  };
    38   Matrix rowVectorMatrix(Vector<typ> const &v)
     37  static Matrix rowVectorMatrix(Vector<typ> const &v)
    3938  {
    4039    Matrix ret(1,v.size());
    41     for(unsigned i=0;i<v.size();i++)ret[0][i]=v[i];
     40    for(int i=0;i<v.size();i++)ret[0][i]=v[i];
    4241    return ret;
    4342  }
     
    4746      assert(i<getWidth());
    4847      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];
    5049      return ret;
    5150    }
     
    5453      Matrix ret(getWidth(),getHeight());
    5554      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];
    5757      return ret;
    5858    }
     
    6060    {
    6161      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);
    6363      return m;
    6464    }
    6565  void append(Matrix const &m)
    6666    {
     67      assert(m.getWidth()==width);
     68          data.resize((height+m.height)*width);
     69          int oldHeight=height;
     70      height+=m.height;
    6771      for(int i=0;i<m.height;i++)
    6872        {
    69           rows.push_back(m[i]);
     73          for(int j=0;j<m.width;j++)
     74                  (*this)[i+oldHeight][j]=m[i][j];
    7075        }
    71       height+=m.height;
    7276    }
    7377  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    }
    8585  void eraseLastRow()
    8686  {
    87     assert(rows.size()>0);
    88     rows.pop_back();
     87    assert(height>0);
     88    data.resize((height-1)*width);
    8989    height--;
    9090  }
     
    111111    {
    112112      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]);
    114116      return p;
    115117    }
     
    147149    for(int i=startRow;i<endRow;i++)
    148150      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
    155317
    156318
     
    164326    for(int i=0;i<getHeight();i++)
    165327      {
    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;
    168330      }
    169331    return false;
     
    180342    if(!a.isZero())
    181343    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);
    184346  }
    185347
     
    189351      {
    190352        if(i)f<<","<<std::endl;
    191         f<<a.rows[i];
     353        f<<a[i].toVector();
    192354      }
    193355    f<<"}"<<std::endl;
     
    197359  std::string toString()const
    198360  {
    199           std::stringstream f;
    200           f<<*this;
    201           return f.str();
     361          std::stringstream f;
     362          f<<*this;
     363          return f.str();
    202364  }
    203365
     
    207369  void swapRows(int i, int j)
    208370  {
    209     std::swap(rows[i],rows[j]);
     371    for(int a=0;a<width;a++)std::swap((*this)[i][a],(*this)[j][a]);
    210372  }
    211373  /**
     
    222384    while(++j<width)
    223385      {
    224         if(!rows[i][j].isZero()) return true;
     386        if(!(*this)[i][j].isZero()) return true;
    225387      }
    226388    return false;
     
    295457    int bestNumberOfNonZero=0;
    296458    for(int i=currentRow;i<height;i++)
    297       if(!rows[i][column].isZero())
     459      if(!(*this)[i][column].isZero())
    298460        {
    299461          int nz=0;
    300462          for(int k=column+1;k<width;k++)
    301             if(!rows[i][k].isZero())nz++;
     463            if(!(*this)[i][k].isZero())nz++;
    302464          if(best==-1)
    303465            {
     
    345507            if(makePivotsOne)
    346508              {//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];
    349511                //                if(!rows[currentRow][i].isOne())
    350512                  {
    351513                    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;
    354516                  }
    355517              }
     
    357519              if(integral)
    358520                {
    359                   if(!rows[j][i].isZero())
     521                  if(!(*this)[j][i].isZero())
    360522                    {
    361523                      typ s;
    362524                      typ t;
    363525
    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;
    367529                        /* We want the (s,t) vector to be as small as possible.
    368530                         * We are allowed to adjust by multiples of (u,v).
     
    378540                        for(int k=0;k<width;k++)
    379541                          {
    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;
    385547                          }
    386548                      }
     
    388550                else
    389551                  {
    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);
    392554                  }
    393555              currentRow++;
     
    412574    while(nextPivot(pivotI,pivotJ))
    413575      {
    414         if(scalePivotsToOne)
    415           rows[pivotI]=rows[pivotI]/rows[pivotI][pivotJ];
     576        if(scalePivotsToOne)
     577          (*this)[pivotI]=(*this)[pivotI].toVector()/(*this)[pivotI][pivotJ];
    416578        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        }
    420582    return ret;
    421583  }
     
    438600      if(!v[pivotJ].isZero())
    439601      {
    440         typ s=-v[pivotJ]/rows[pivotI][pivotJ];
     602        typ s=-v[pivotJ]/(*this)[pivotI][pivotJ];
    441603
    442604        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);
    445607      }
    446608    return v;
     
    472634        while(nextPivot(pivot2I,pivot2J))
    473635          {
    474             ret[k][pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J];
     636            ret[k][pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J];
    475637          }
    476638        ret[k][j]=typ(-1);
     
    508670      while(nextPivot(pivot2I,pivot2J))
    509671        {
    510           diagonalProduct*=rows[pivot2I][pivot2J];
     672          diagonalProduct*=(*this)[pivot2I][pivot2J];
    511673        }
    512674    }
     
    522684      while(nextPivot(pivot2I,pivot2J))
    523685        {
    524           ret[pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J];
     686          ret[pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J];
    525687          lastEntry-=ret[pivot2J]*ret[pivot2J];
    526688        }
     
    546708   * Sort the rows of the matrix.
    547709   */
     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;
    548713  void sortRows()
    549714  {
    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;
    551722  }
    552723  /**
     
    560731    B.appendRow((*this)[0]);
    561732    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());
    563734    *this=B;
    564735  }
     
    573744    assert(bottom.getWidth()==top.getWidth());
    574745    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];
    577748
    578749    return ret;
     
    590761    for(int i=0;i<left.getHeight();i++)
    591762      {
    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];
    594765      }
    595766    return ret;
  • gfanlib/gfanlib_polyhedralfan.cpp

    r6a97c54 r5cea7a  
    7171int PolyhedralFan::getMaxDimension()const
    7272{
    73   if (cones.empty())
    74     return -1;
     73  assert(!cones.empty());
    7574
    7675  return cones.begin()->dimension();
     
    7978int PolyhedralFan::getMinDimension()const
    8079{
    81   if (cones.empty())
    82     return -1;
     80  assert(!cones.empty());
    8381
    8482  return cones.rbegin()->dimension();
     
    355353      bool notAll=false;
    356354      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)
    358356          indices.insert(theCone.indices[j]);
    359357        else
     
    864862int PolyhedralFan::dimensionOfLinealitySpace()const
    865863{
    866   if (cones.empty())
    867     return n;
     864  assert(cones.size());//slow!
    868865  return cones.begin()->dimensionOfLinealitySpace();
    869866}
  • gfanlib/gfanlib_q.h

    r6a97c54 r5cea7a  
    1111#include <string.h>
    1212#include <ostream>
     13#include <assert.h>
    1314#include "gmp.h"
    1415
     
    160161    return mpq_sgn(value);
    161162  }
    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)
    163164  {
    164165/*    mpz_t r;
     
    167168    Integer ret(r);
    168169    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");
    180171    return a;
    181172  }
  • gfanlib/gfanlib_symmetriccomplex.cpp

    r6a97c54 r5cea7a  
    185185int SymmetricComplex::getLinDim()const
    186186{
    187   return linealitySpace.getHeight();
     187        return linealitySpace.getHeight();
    188188}
    189189
     
    595595
    596596    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));
    599598    if(flags&FPF_conesCompressed)polymakeFile.writeStringProperty("CONES_ORBITS",toStringJustCones(getMinDim(),getMaxDim(),false,flags&FPF_group, 0,true,flags&FPF_tPlaneSort));
    600599    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());
    602600
    603601    if(!sym.isTrivial())
     
    637635            stringstream multiplicities;
    638636            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());
    639638  //          log1 fprintf(Stderr,"Done producing list of maximal cones.\n");
    640639          }
  • gfanlib/gfanlib_vector.h

    r6a97c54 r5cea7a  
    7272  // Access
    7373  //--------
    74   typ& operator[](int n)
     74  typ& operator[](int n)__attribute__((always_inline))
    7575    {
    7676      if(!(n>=0 && n<(int)v.size()))outOfRange(n,v.size());
    7777      return (v[n]);
    7878    }
    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]);}
    8182
    8283  //-------------
     
    8586  unsigned int size()const{return v.size();};
    8687  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);}
    8889  void push_back(typ a)
    8990  {
     
    293294  std::string toString()const
    294295  {
    295     std::stringstream f;
    296     f<<*this;
    297     return f.str();
     296          std::stringstream f;
     297          f<<*this;
     298          return f.str();
    298299  }
    299300
    300301  typ gcd()const
    301302  {
     303    typ temp1,temp2;
    302304    typ ret(1);
    303305    for(unsigned i=0;i<size();i++)
    304       ret=typ::gcd(ret,v[i]);
     306      ret=typ::gcd(ret,v[i],temp1,temp2);
    305307    return ret;
    306308  }
     
    309311    assert(!typ::isField());
    310312    return (*this)/gcd();
    311   }
    312 
    313   void debugPrint()const
    314   {
    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;
    325313  }
    326314};
     
    415403
    416404    return ret;
    417 };
    418 
    419 };
     405}
     406
     407}
    420408
    421409#endif /* LIB_ZVECTOR_H_ */
  • gfanlib/gfanlib_z.h

    r6a97c54 r5cea7a  
    1010
    1111#include <string.h>
    12 #include <iostream>
    1312#include <ostream>
    1413
     
    4140    mpz_init_set(value,value_.value);
    4241  }
    43   Integer(mpz_t value_)
     42  Integer(mpz_t const value_)
    4443  {
    4544    mpz_init_set(value,value_);
     
    6059  bool isZero()const{
    6160    return mpz_sgn(value)==0;
    62   }
    63   bool isOne()const{
    64     return mpz_cmp_si(value,1);
    6561  }
    6662  friend std::ostream &operator<<(std::ostream &f, Integer const &a)
     
    7369    return f;
    7470  }
    75   void debugPrint() const
    76   {
    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   }
    8471  Integer& operator+=(const Integer& a)
    8572    {
     
    161148    mpz_init(r);
    162149    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);
    174150    Integer ret(r);
    175151    mpz_clear(r);
     
    354330      return f;
    355331    }
    356     friend void debugPrint(IntegerTemplate const &a)
    357     {
    358       std::cout << a << std::endl;
    359       return;
    360     }
    361332    LimbWord signExtension(LimbWord a)
    362333    {
  • gfanlib/gfanlib_zcone.cpp

    r6a97c54 r5cea7a  
    1212#include <sstream>
    1313
    14 #include <config.h>
    15 #ifdef HAVE_CDD_SETOPER_H
    16 #include <cdd/setoper.h>
    17 #include <cdd/cdd.h>
     14//extern "C"{
     15#ifdef NOCDDPREFIX
     16#include "setoper.h"
     17#include "cdd.h"
    1818#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//}
    2723
    2824namespace 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)
    3469      {
    35         dd_set_global_constants();  /* First, this must be called. */
    36         initialized=true;
    37       }
     70                        dd_set_global_constants();
     71                        initialized=true;
     72      }*/
    3873  }
    3974
     
    85120    int numberOfRows=g.getHeight();
    86121    dd_MatrixPtr A=NULL;
    87     cddinitGmp();
     122    ensureCddInitialisation();
    88123    A=ZMatrix2MatrixGmp(g, err);
    89124    for(int i=numberOfInequalities;i<numberOfRows;i++)
     
    121156    FILE *reading=NULL;
    122157
    123     cddinitGmp();
     158    ensureCddInitialisation();
    124159
    125160    M=ZMatrix2MatrixGmp(g, &err);
     
    128163    // d=M->colsize;
    129164
    130     static dd_Arow temp;
     165    //static dd_Arow temp;
     166    dd_Arow temp;
    131167    dd_InitializeArow(g.getWidth()+1,&temp);
    132168
     
    268304    for(int i=0;i<a.getHeight();i++)
    269305      {
    270         assert(!(a[i].isZero()));
    271         b.insert(a[i].normalized());
     306        assert(!(a[i].toVector().isZero()));
     307        b.insert(a[i].toVector().normalized());
    272308      }
    273309
     
    302338    for(int i=0;i!=original.getHeight();i++)
    303339      for(int j=0;j!=a.getHeight();j++)
    304         if(!dependent(original[i],a[j]))
     340        if(!dependent(original[i].toVector(),a[j].toVector()))
    305341            {
    306342              ZVector const &I=original[i];
     
    337373  void removeRedundantRows(ZMatrix &inequalities, ZMatrix &equations, bool removeInequalityRedundancies)
    338374  {
    339     cddinitGmp();
     375          ensureCddInitialisation();
    340376
    341377    int numberOfEqualities=equations.getHeight();
     
    403439  {
    404440    QVector retUnscaled(inequalities.getWidth());
    405     cddinitGmp();
     441    ensureCddInitialisation();
    406442    int numberOfEqualities=equations.getHeight();
    407443    int numberOfInequalities=inequalities.getHeight();
     
    459495    dd_ErrorType err=dd_NoError;
    460496
    461         cddinitGmp();
     497    ensureCddInitialisation();
    462498
    463499    A=ZMatrix2MatrixGmp(inequalities, equations, &err);
     
    567603    dd_ErrorType err=dd_NoError;
    568604
    569     cddinitGmp();
     605    ensureCddInitialisation();
    570606    A=ZMatrix2MatrixGmp(inequalities, &err);
    571607
     
    748784std::string ZCone::toString()const
    749785{
    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();
    783789}
    784790
     
    983989  }
    984990*/
    985 //          ZCone dual(newGenerators,linealitySpace);
    986           ZCone dual(generators,linealitySpace);
     991//        ZCone dual(newGenerators,linealitySpace);
     992          ZCone dual(generators,linealitySpace);
    987993//  dual.findFacets();
    988994//  dual.canonicalize();
     
    12381244
    12391245  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())
    12411247      {
    1242         ret.appendRow(M[i].subvector(a+b,a+b+n));
     1248        ret.appendRow(M[i].toVector().subvector(a+b,a+b+n));
    12431249      }
    12441250
     
    12521258  assert(temp.getHeight()==1);
    12531259  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)
    12551261      {
    1256         temp[0]=-temp[0];
     1262        temp[0]=-temp[0].toVector();
    12571263        break;
    12581264      }
     
    13381344}
    13391345
    1340 };
     1346}
  • gfanlib/gfanlib_zcone.h

    r6a97c54 r5cea7a  
    1010
    1111#include "gfanlib_matrix.h"
    12 #include <sstream>
    1312
    1413namespace 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();
    1628
    1729/**
     
    358370};
    359371
    360 };
     372}
     373
     374
    361375
    362376
Note: See TracChangeset for help on using the changeset viewer.