Public Member Functions | Private Member Functions | Private Attributes
resMatrixSparse Class Reference

Public Member Functions

 resMatrixSparse (const ideal _gls, const int special=SNONE)
 
 ~resMatrixSparse ()
 
ideal getMatrix ()
 
number getDetAt (const number *evpoint)
 Fills in resMat[][] with evpoint[] and gets determinant uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) . More...
 
poly getUDet (const number *evpoint)
 
- Public Member Functions inherited from resMatrixBase
 resMatrixBase ()
 
virtual ~resMatrixBase ()
 
virtual ideal getSubMatrix ()
 
virtual number getSubDet ()
 
virtual long getDetDeg ()
 
virtual IStateType initState () const
 

Private Member Functions

 resMatrixSparse (const resMatrixSparse &)
 
void randomVector (const int dim, mprfloat shift[])
 
int RC (pointSet **pQ, pointSet *E, int vert, mprfloat shift[])
 Row Content Function Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j. More...
 
bool remapXiToPoint (const int indx, pointSet **pQ, int *set, int *vtx)
 
int createMatrix (pointSet *E)
 create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) . More...
 
pointSetminkSumAll (pointSet **pQ, int numq, int dim)
 
pointSetminkSumTwo (pointSet *Q1, pointSet *Q2, int dim)
 

Private Attributes

ideal gls
 
int n
 
int idelem
 
int numSet0
 
int msize
 
intvecuRPos
 
ideal rmat
 
simplexLP
 

Additional Inherited Members

- Public Types inherited from resMatrixBase
enum  IStateType {
  none, ready, notInit, fatalError,
  sparseError
}
 
- Protected Attributes inherited from resMatrixBase
IStateType istate
 
ideal gls
 
int linPolyS
 
ring sourceRing
 
int totDeg
 

Detailed Description

Definition at line 69 of file mpr_base.cc.

Constructor & Destructor Documentation

resMatrixSparse::resMatrixSparse ( const ideal  _gls,
const int  special = SNONE 
)

Definition at line 1574 of file mpr_base.cc.

1575  : resMatrixBase(), gls( _gls )
1576 {
1577  pointSet **Qi; // vertices sets of Conv(Supp(f_i)), i=0..idelem
1578  pointSet *E; // all integer lattice points of the minkowski sum of Q0...Qn
1579  int i,k;
1580  int pnt;
1581  int totverts; // total number of exponent vectors in ideal gls
1582  mprfloat shift[MAXVARS+2]; // shiftvector delta, index [1..dim]
1583 
1584  if ( (currRing->N) > MAXVARS )
1585  {
1586  WerrorS("resMatrixSparse::resMatrixSparse: Too many variables!");
1587  return;
1588  }
1589 
1590  rmat= NULL;
1591  numSet0= 0;
1592 
1593  if ( special == SNONE ) linPolyS= 0;
1594  else linPolyS= special;
1595 
1597 
1598  n= (currRing->N);
1599  idelem= IDELEMS(gls); // should be n+1
1600 
1601  // prepare matrix LP->LiPM for Linear Programming
1602  totverts = 0;
1603  for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
1604 
1605  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
1606 
1607  // get shift vector
1608 #ifdef mprTEST
1609  shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1610  shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1611 #else
1612  randomVector( idelem, shift );
1613 #endif
1614 #ifdef mprDEBUG_PROT
1615  PrintS(" shift vector: ");
1616  for ( i= 1; i <= idelem; i++ ) Print(" %.12f ",(double)shift[i]);
1617  PrintLn();
1618 #endif
1619 
1620  // evaluate convex hull for supports of gls
1621  convexHull chnp( LP );
1622  Qi= chnp.newtonPolytopesP( gls );
1623 
1624 #ifdef mprMINKSUM
1625  E= minkSumAll( Qi, n+1, n);
1626 #else
1627  // get inner points
1628  mayanPyramidAlg mpa( LP );
1629  E= mpa.getInnerPoints( Qi, shift );
1630 #endif
1631 
1632 #ifdef mprDEBUG_PROT
1633 #ifdef mprMINKSUM
1634  PrintS("(MinkSum)");
1635 #endif
1636  PrintS("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1637  for ( pnt= 1; pnt <= E->num; pnt++ )
1638  {
1639  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1640  }
1641  PrintLn();
1642 #endif
1643 
1644 #ifdef mprTEST
1645  int lift[5][5];
1646  lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1647  lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1648  lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1649  lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1650  lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1651  // now lift everything
1652  for ( i= 0; i <= n; i++ ) Qi[i]->lift( lift[i] );
1653 #else
1654  // now lift everything
1655  for ( i= 0; i <= n; i++ ) Qi[i]->lift();
1656 #endif
1657  E->dim++;
1658 
1659  // run Row Content Function for every point in E
1660  for ( pnt= 1; pnt <= E->num; pnt++ )
1661  {
1662  RC( Qi, E, pnt, shift );
1663  }
1664 
1665  // remove points not in cells
1666  k= E->num;
1667  for ( pnt= k; pnt > 0; pnt-- )
1668  {
1669  if ( (*E)[pnt]->rcPnt == NULL )
1670  {
1671  E->removePoint(pnt);
1673  }
1674  }
1675  mprSTICKYPROT("\n");
1676 
1677 #ifdef mprDEBUG_PROT
1678  PrintS(" points which lie in a cell:\n");
1679  for ( pnt= 1; pnt <= E->num; pnt++ )
1680  {
1681  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1682  }
1683  PrintLn();
1684 #endif
1685 
1686  // unlift to old dimension, sort
1687  for ( i= 0; i <= n; i++ ) Qi[i]->unlift();
1688  E->unlift();
1689  E->sort();
1690 
1691 #ifdef mprDEBUG_PROT
1692  Print(" points with a[ij] (%d):\n",E->num);
1693  for ( pnt= 1; pnt <= E->num; pnt++ )
1694  {
1695  Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
1696  Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1697  //print_exp( (Qi[(*E)[pnt]->rc.set])[(*E)[pnt]->rc.pnt], E->dim );PrintS("> = <");
1698  print_exp( (*E)[pnt]->rcPnt, E->dim );PrintS(">\n");
1699  }
1700 #endif
1701 
1702  // now create matrix
1703  if (E->num <1)
1704  {
1705  WerrorS("could not handle a degenerate situation: no inner points found");
1706  goto theEnd;
1707  }
1708  if ( createMatrix( E ) != E->num )
1709  {
1710  // this can happen if the shiftvector shift is to large or not generic
1712  WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1713  goto theEnd;
1714  }
1715 
1716  theEnd:
1717  // clean up
1718  for ( i= 0; i < idelem; i++ )
1719  {
1720  delete Qi[i];
1721  }
1722  omFreeSize( (void *) Qi, idelem * sizeof(pointSet*) );
1723 
1724  delete E;
1725 
1726  delete LP;
1727 }
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
int RC(pointSet **pQ, pointSet *E, int vert, mprfloat shift[])
Row Content Function Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j...
Definition: mpr_base.cc:1240
void randomVector(const int dim, mprfloat shift[])
Definition: mpr_base.cc:1504
void PrintLn()
Definition: reporter.cc:327
#define Print
Definition: emacs.cc:83
int dim
Definition: mpr_base.cc:172
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
#define SNONE
Definition: mpr_base.h:14
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
double mprfloat
Definition: mpr_global.h:17
static int pLength(poly a)
Definition: p_polys.h:189
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
Definition: lift.cc:24
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
void unlift()
Definition: mpr_base.cc:232
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
#define IDELEMS(i)
Definition: simpleideals.h:24
#define NULL
Definition: omList.c:10
int num
Definition: mpr_base.cc:170
REvaluation E(1, terms.length(), IntRandom(25))
#define MAXVARS
Definition: mpr_base.cc:58
bool removePoint(const int indx)
Definition: mpr_base.cc:500
int linPolyS
Definition: mpr_base.h:47
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
Definition: mpr_base.cc:1552
IStateType istate
Definition: mpr_base.h:44
void sort()
sort lex
Definition: mpr_base.cc:650
#define ST_SPARSE_RCRJ
Definition: mpr_global.h:76
int createMatrix(pointSet *E)
create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) .
Definition: mpr_base.cc:1412
simplex * LP
Definition: mpr_base.cc:127
resMatrixSparse::~resMatrixSparse ( )

Definition at line 1731 of file mpr_base.cc.

1732 {
1733  delete uRPos;
1734  idDelete( &rmat );
1735 }
#define idDelete(H)
delete an ideal
Definition: ideals.h:31
intvec * uRPos
Definition: mpr_base.cc:123
resMatrixSparse::resMatrixSparse ( const resMatrixSparse )
private

Member Function Documentation

int resMatrixSparse::createMatrix ( pointSet E)
private

create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) .

. u(n) i= 1 .. numSet0 Returns the dimension of the matrix or -1 in case of an error

Definition at line 1412 of file mpr_base.cc.

1413 {
1414  // sparse matrix
1415  // uRPos[i][1]: row of matrix
1416  // uRPos[i][idelem+1]: col of u(0)
1417  // uRPos[i][2..idelem]: col of u(1) .. u(n)
1418  // i= 1 .. numSet0
1419  int i,epos;
1420  int rp,cp;
1421  poly rowp,epp;
1422  poly iterp;
1423  int *epp_mon, *eexp;
1424 
1425  epp_mon= (int *)omAlloc( (n+2) * sizeof(int) );
1426  eexp= (int *)omAlloc0(((currRing->N)+1)*sizeof(int));
1427 
1428  totDeg= numSet0;
1429 
1430  mprSTICKYPROT2(" size of matrix: %d\n", E->num);
1431  mprSTICKYPROT2(" resultant deg: %d\n", numSet0);
1432 
1433  uRPos= new intvec( numSet0, pLength((gls->m)[0])+1, 0 );
1434 
1435  // sparse Matrix represented as a module where
1436  // each poly is column vector ( pSetComp(p,k) gives the row )
1437  rmat= idInit( E->num, E->num ); // cols, rank= number of rows
1438  msize= E->num;
1439 
1440  rp= 1;
1441  rowp= NULL;
1442  epp= pOne();
1443  for ( i= 1; i <= E->num; i++ )
1444  { // for every row
1445  E->getRowMP( i, epp_mon ); // compute (p-a[ij]), (i,j) = RC(p)
1446  pSetExpV( epp, epp_mon );
1447 
1448  //
1449  rowp= ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] ); // x^(p-a[ij]) * f(i)
1450 
1451  cp= 2;
1452  // get column for every monomial in rowp and store it
1453  iterp= rowp;
1454  while ( iterp!=NULL )
1455  {
1456  epos= E->getExpPos( iterp );
1457  if ( epos == 0 )
1458  {
1459  // this can happen, if the shift vektor or the lift funktions
1460  // are not generically chosen.
1461  Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1462  i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1463  return i;
1464  }
1465  pSetExpV(iterp,eexp);
1466  pSetComp(iterp, epos );
1467  pSetm(iterp);
1468  if ( (*E)[i]->rc.set == linPolyS )
1469  { // store coeff positions
1470  IMATELEM(*uRPos,rp,cp)= epos;
1471  cp++;
1472  }
1473  pIter( iterp );
1474  } // while
1475  if ( (*E)[i]->rc.set == linPolyS )
1476  { // store row
1477  IMATELEM(*uRPos,rp,1)= i-1;
1478  rp++;
1479  }
1480  (rmat->m)[i-1]= rowp;
1481  } // for
1482 
1483  pDelete( &epp );
1484  omFreeSize( (void *) epp_mon, (n+2) * sizeof(int) );
1485  omFreeSize( (void *) eexp, ((currRing->N)+1)*sizeof(int));
1486 
1487 #ifdef mprDEBUG_ALL
1488  if ( E->num <= 40 )
1489  {
1490  matrix mout= idModule2Matrix( idCopy(rmat) );
1491  print_matrix(mout);
1492  }
1493  for ( i= 1; i <= numSet0; i++ )
1494  {
1495  Print(" row %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
1496  }
1497  PrintS(" Sparse Matrix done\n");
1498 #endif
1499 
1500  return E->num;
1501 }
#define ppMult_qq(p, q)
Definition: polys.h:179
#define pSetm(p)
Definition: polys.h:241
#define Print
Definition: emacs.cc:83
#define mprSTICKYPROT2(msg, arg)
Definition: mpr_global.h:55
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
void getRowMP(const int indx, int *vert)
Definition: mpr_base.cc:602
#define omAlloc(size)
Definition: omAllocDecl.h:210
static int pLength(poly a)
Definition: p_polys.h:189
#define pIter(p)
Definition: monomials.h:44
int getExpPos(const poly p)
Definition: mpr_base.cc:581
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
Definition: intvec.h:14
#define pSetExpV(p, e)
Definition: polys.h:97
#define pSetComp(p, v)
Definition: polys.h:38
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
#define pOne()
Definition: polys.h:286
ideal idCopy(ideal A)
Definition: ideals.h:73
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
#define NULL
Definition: omList.c:10
int num
Definition: mpr_base.cc:170
int linPolyS
Definition: mpr_base.h:47
#define pDelete(p_ptr)
Definition: polys.h:157
polyrec * poly
Definition: hilb.h:10
#define IMATELEM(M, I, J)
Definition: intvec.h:77
void Werror(const char *fmt,...)
Definition: reporter.cc:199
#define omAlloc0(size)
Definition: omAllocDecl.h:211
intvec * uRPos
Definition: mpr_base.cc:123
number resMatrixSparse::getDetAt ( const number *  evpoint)
virtual

Fills in resMat[][] with evpoint[] and gets determinant uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2..idelem]: col of u(1) .

. u(n) i= 1 .. numSet0

Reimplemented from resMatrixBase.

Definition at line 1798 of file mpr_base.cc.

1799 {
1800  int i,cp;
1801  poly pp,phelp,piter;
1802 
1803  mprPROTnl("smCallDet");
1804 
1805  for ( i= 1; i <= numSet0; i++ )
1806  {
1807  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1808  pDelete( &pp );
1809  pp= NULL;
1810  phelp= pp;
1811  piter= NULL;
1812  // u_1,..,u_n
1813  for ( cp= 2; cp <= idelem; cp++ )
1814  {
1815  if ( !nIsZero(evpoint[cp-1]) )
1816  {
1817  phelp= pOne();
1818  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1819  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1820  pSetmComp( phelp );
1821  if ( piter )
1822  {
1823  pNext(piter)= phelp;
1824  piter= phelp;
1825  }
1826  else
1827  {
1828  pp= phelp;
1829  piter= phelp;
1830  }
1831  }
1832  }
1833  // u0
1834  phelp= pOne();
1835  pSetCoeff( phelp, nCopy(evpoint[0]) );
1836  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1837  pSetmComp( phelp );
1838  pNext(piter)= phelp;
1839  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1840  }
1841 
1842  mprSTICKYPROT(ST__DET); // 1
1843 
1844  poly pres= sm_CallDet( rmat, currRing );
1845  number numres= nCopy( pGetCoeff( pres ) );
1846  pDelete( &pres );
1847 
1848  mprSTICKYPROT(ST__DET); // 2
1849 
1850  return ( numres );
1851 }
#define pSetmComp(p)
TODO:
Definition: polys.h:243
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:359
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
#define mprPROTnl(msg)
Definition: mpr_global.h:42
poly pp
Definition: myNF.cc:296
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
#define pSetComp(p, v)
Definition: polys.h:38
int i
Definition: cfEzgcd.cc:123
#define pOne()
Definition: polys.h:286
#define ST__DET
Definition: mpr_global.h:78
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:10
#define pDelete(p_ptr)
Definition: polys.h:157
#define nCopy(n)
Definition: numbers.h:15
#define pNext(p)
Definition: monomials.h:43
#define phelp
Definition: libparse.cc:1240
polyrec * poly
Definition: hilb.h:10
#define IMATELEM(M, I, J)
Definition: intvec.h:77
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
intvec * uRPos
Definition: mpr_base.cc:123
ideal resMatrixSparse::getMatrix ( )
virtual

Reimplemented from resMatrixBase.

Definition at line 1737 of file mpr_base.cc.

1738 {
1739  int i,/*j,*/cp;
1740  poly pp,phelp,piter,pgls;
1741 
1742  // copy original sparse res matrix
1743  ideal rmat_out= idCopy(rmat);
1744 
1745  // now fill in coeffs of f0
1746  for ( i= 1; i <= numSet0; i++ )
1747  {
1748 
1749  pgls= (gls->m)[0]; // f0
1750 
1751  // get matrix row and delete it
1752  pp= (rmat_out->m)[IMATELEM(*uRPos,i,1)];
1753  pDelete( &pp );
1754  pp= NULL;
1755  phelp= pp;
1756  piter= NULL;
1757 
1758  // u_1,..,u_k
1759  cp=2;
1760  while ( pNext(pgls)!=NULL )
1761  {
1762  phelp= pOne();
1763  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1764  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1765  pSetmComp( phelp );
1766  if ( piter!=NULL )
1767  {
1768  pNext(piter)= phelp;
1769  piter= phelp;
1770  }
1771  else
1772  {
1773  pp= phelp;
1774  piter= phelp;
1775  }
1776  cp++;
1777  pIter( pgls );
1778  }
1779  // u0, now pgls points to last monom
1780  phelp= pOne();
1781  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1782  //pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1783  pSetComp( phelp, IMATELEM(*uRPos,i,pLength((gls->m)[0])+1) );
1784  pSetmComp( phelp );
1785  if (piter!=NULL) pNext(piter)= phelp;
1786  else pp= phelp;
1787  (rmat_out->m)[IMATELEM(*uRPos,i,1)]= pp;
1788  }
1789 
1790  return rmat_out;
1791 }
#define pSetmComp(p)
TODO:
Definition: polys.h:243
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
static int pLength(poly a)
Definition: p_polys.h:189
poly pp
Definition: myNF.cc:296
#define pIter(p)
Definition: monomials.h:44
#define pSetComp(p, v)
Definition: polys.h:38
int i
Definition: cfEzgcd.cc:123
#define pOne()
Definition: polys.h:286
ideal idCopy(ideal A)
Definition: ideals.h:73
#define NULL
Definition: omList.c:10
#define pDelete(p_ptr)
Definition: polys.h:157
#define nCopy(n)
Definition: numbers.h:15
#define pNext(p)
Definition: monomials.h:43
#define phelp
Definition: libparse.cc:1240
polyrec * poly
Definition: hilb.h:10
#define IMATELEM(M, I, J)
Definition: intvec.h:77
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
intvec * uRPos
Definition: mpr_base.cc:123
poly resMatrixSparse::getUDet ( const number *  evpoint)
virtual

Reimplemented from resMatrixBase.

Definition at line 1858 of file mpr_base.cc.

1859 {
1860  int i,cp;
1861  poly pp,phelp/*,piter*/;
1862 
1863  mprPROTnl("smCallDet");
1864 
1865  for ( i= 1; i <= numSet0; i++ )
1866  {
1867  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1868  pDelete( &pp );
1869  phelp= NULL;
1870  // piter= NULL;
1871  for ( cp= 2; cp <= idelem; cp++ )
1872  { // u1 .. un
1873  if ( !nIsZero(evpoint[cp-1]) )
1874  {
1875  phelp= pOne();
1876  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1877  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1878  //pSetmComp( phelp );
1879  pSetm( phelp );
1880  //Print("comp %d\n",IMATELEM(*uRPos,i,cp));
1881  #if 0
1882  if ( piter!=NULL )
1883  {
1884  pNext(piter)= phelp;
1885  piter= phelp;
1886  }
1887  else
1888  {
1889  pp= phelp;
1890  piter= phelp;
1891  }
1892  #else
1893  pp=pAdd(pp,phelp);
1894  #endif
1895  }
1896  }
1897  // u0
1898  phelp= pOne();
1899  pSetExp(phelp,1,1);
1900  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1901  // Print("comp %d\n",IMATELEM(*uRPos,i,idelem+1));
1902  pSetm( phelp );
1903  #if 0
1904  pNext(piter)= phelp;
1905  #else
1906  pp=pAdd(pp,phelp);
1907  #endif
1908  pTest(pp);
1909  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1910  }
1911 
1912  mprSTICKYPROT(ST__DET); // 1
1913 
1914  poly pres= sm_CallDet( rmat, currRing );
1915 
1916  mprSTICKYPROT(ST__DET); // 2
1917 
1918  return ( pres );
1919 }
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
#define pSetm(p)
Definition: polys.h:241
#define pAdd(p, q)
Definition: polys.h:174
#define pSetExp(p, i, v)
Definition: polys.h:42
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:359
#define pTest(p)
Definition: polys.h:387
#define mprPROTnl(msg)
Definition: mpr_global.h:42
poly pp
Definition: myNF.cc:296
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
#define pSetComp(p, v)
Definition: polys.h:38
int i
Definition: cfEzgcd.cc:123
#define pOne()
Definition: polys.h:286
#define ST__DET
Definition: mpr_global.h:78
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:10
#define pDelete(p_ptr)
Definition: polys.h:157
#define nCopy(n)
Definition: numbers.h:15
#define pNext(p)
Definition: monomials.h:43
#define phelp
Definition: libparse.cc:1240
polyrec * poly
Definition: hilb.h:10
#define IMATELEM(M, I, J)
Definition: intvec.h:77
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
intvec * uRPos
Definition: mpr_base.cc:123
pointSet * resMatrixSparse::minkSumAll ( pointSet **  pQ,
int  numq,
int  dim 
)
private

Definition at line 1552 of file mpr_base.cc.

1553 {
1554  pointSet *vs,*vs_old;
1555  int j;
1556 
1557  vs= new pointSet( dim );
1558 
1559  for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
1560 
1561  for ( j= 1; j < numq; j++ )
1562  {
1563  vs_old= vs;
1564  vs= minkSumTwo( vs_old, pQ[j], dim );
1565 
1566  delete vs_old;
1567  }
1568 
1569  return vs;
1570 }
int j
Definition: myNF.cc:70
int dim(ideal I, ring r)
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
Definition: mpr_base.cc:1524
int num
Definition: mpr_base.cc:170
bool addPoint(const onePointP vert)
Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
Definition: mpr_base.cc:467
pointSet * resMatrixSparse::minkSumTwo ( pointSet Q1,
pointSet Q2,
int  dim 
)
private

Definition at line 1524 of file mpr_base.cc.

1525 {
1526  pointSet *vs;
1527  onePoint vert;
1528  int j,k,l;
1529 
1530  vert.point=(Coord_t*)omAlloc( ((currRing->N)+2) * sizeof(Coord_t) );
1531 
1532  vs= new pointSet( dim );
1533 
1534  for ( j= 1; j <= Q1->num; j++ )
1535  {
1536  for ( k= 1; k <= Q2->num; k++ )
1537  {
1538  for ( l= 1; l <= dim; l++ )
1539  {
1540  vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
1541  }
1542  vs->mergeWithExp( &vert );
1543  //vs->addPoint( &vert );
1544  }
1545  }
1546 
1547  omFreeSize( (void *) vert.point, ((currRing->N)+2) * sizeof(Coord_t) );
1548 
1549  return vs;
1550 }
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int k
Definition: cfEzgcd.cc:93
#define omAlloc(size)
Definition: omAllocDecl.h:210
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
int j
Definition: myNF.cc:70
int dim(ideal I, ring r)
Coord_t * point
Definition: mpr_base.cc:144
int num
Definition: mpr_base.cc:170
unsigned int Coord_t
Definition: mpr_base.cc:134
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet point = .
Definition: mpr_base.cc:515
int l
Definition: cfEzgcd.cc:94
void resMatrixSparse::randomVector ( const int  dim,
mprfloat  shift[] 
)
private

Definition at line 1504 of file mpr_base.cc.

1505 {
1506  int i,j;
1507  i= 1;
1508 
1509  while ( i <= dim )
1510  {
1511  shift[i]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
1512  i++;
1513  for ( j= 1; j < i-1; j++ )
1514  {
1515  if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
1516  {
1517  i--;
1518  break;
1519  }
1520  }
1521  }
1522 }
double mprfloat
Definition: mpr_global.h:17
int j
Definition: myNF.cc:70
#define MAXRVVAL
Definition: mpr_base.cc:57
int dim(ideal I, ring r)
int i
Definition: cfEzgcd.cc:123
#define SIMPLEX_EPS
Definition: mpr_numeric.h:181
#define RVMULT
Definition: mpr_base.cc:56
int siRand()
Definition: sirandom.c:41
int resMatrixSparse::RC ( pointSet **  pQ,
pointSet E,
int  vert,
mprfloat  shift[] 
)
private

Row Content Function Finds the largest i such that F[i] is a point, F[i]= a[ij] in A[i] for some j.

Returns -1 iff the point vert does not lie in a cell

Definition at line 1240 of file mpr_base.cc.

1241 {
1242  int i, j, k,c ;
1243  int size;
1244  bool found= true;
1245  mprfloat cd;
1246  int onum;
1247  int bucket[MAXVARS+2];
1248  setID *optSum;
1249 
1250  LP->n = 1;
1251  LP->m = n + n + 1; // number of constrains
1252 
1253  // fill in LP matrix
1254  for ( i= 0; i <= n; i++ )
1255  {
1256  size= pQ[i]->num;
1257  for ( k= 1; k <= size; k++ )
1258  {
1259  LP->n++;
1260 
1261  // objective funtion, minimize
1262  LP->LiPM[1][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
1263 
1264  // lambdas sum up to 1
1265  for ( j = 0; j <= n; j++ )
1266  {
1267  if ( i==j )
1268  LP->LiPM[j+2][LP->n] = -1.0;
1269  else
1270  LP->LiPM[j+2][LP->n] = 0.0;
1271  }
1272 
1273  // the points
1274  for ( j = 1; j <= n; j++ )
1275  {
1276  LP->LiPM[j+n+2][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[j] );
1277  }
1278  }
1279  }
1280 
1281  for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
1282  for ( j= 1; j <= n; j++ )
1283  {
1284  LP->LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
1285  }
1286  LP->n--;
1287 
1288  LP->LiPM[1][1] = 0.0;
1289 
1290 #ifdef mprDEBUG_ALL
1291  PrintLn();
1292  Print(" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
1293  print_mat(LP->LiPM, LP->m+1, LP->n+1);
1294 #endif
1295 
1296  LP->m3= LP->m;
1297 
1298  LP->compute();
1299 
1300  if ( LP->icase < 0 )
1301  {
1302  // infeasibility: the point does not lie in a cell -> remove it
1303  return -1;
1304  }
1305 
1306  // store result
1307  (*E)[vert]->point[E->dim]= (int)(-LP->LiPM[1][1] * SCALEDOWN);
1308 
1309 #ifdef mprDEBUG_ALL
1310  Print(" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
1311  //print_bmat(LP->LiPM, NumCons + 1, LP->n+1-NumCons, LP->n+1, LP->iposv); // ( rows= M+1, cols= N+1-m3 )
1312  //print_mat(LP->LiPM, NumCons+1, LP->n);
1313 #endif
1314 
1315 #if 1
1316  // sort LP results
1317  while (found)
1318  {
1319  found=false;
1320  for ( i= 1; i < LP->m; i++ )
1321  {
1322  if ( LP->iposv[i] > LP->iposv[i+1] )
1323  {
1324 
1325  c= LP->iposv[i];
1326  LP->iposv[i]=LP->iposv[i+1];
1327  LP->iposv[i+1]=c;
1328 
1329  cd=LP->LiPM[i+1][1];
1330  LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
1331  LP->LiPM[i+2][1]=cd;
1332 
1333  found= true;
1334  }
1335  }
1336  }
1337 #endif
1338 
1339 #ifdef mprDEBUG_ALL
1340  print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
1341  PrintS(" now split into sets\n");
1342 #endif
1343 
1344 
1345  // init bucket
1346  for ( i= 0; i <= E->dim; i++ ) bucket[i]= 0;
1347  // remap results of LP to sets Qi
1348  c=0;
1349  optSum= (setID*)omAlloc( (LP->m) * sizeof(struct setID) );
1350  for ( i= 0; i < LP->m; i++ )
1351  {
1352  //Print("% .15f\n",LP->LiPM[i+2][1]);
1353  if ( LP->LiPM[i+2][1] > 1e-12 )
1354  {
1355  if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
1356  {
1357  Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
1358  WerrorS(" resMatrixSparse::RC: remapXiToPoint failed!");
1359  return -1;
1360  }
1361  bucket[optSum[c].set]++;
1362  c++;
1363  }
1364  }
1365 
1366  onum= c;
1367  // find last min in bucket[]: maximum i such that Fi is a point
1368  c= 0;
1369  for ( i= 1; i < E->dim; i++ )
1370  {
1371  if ( bucket[c] >= bucket[i] )
1372  {
1373  c= i;
1374  }
1375  }
1376  // find matching point set
1377  for ( i= onum - 1; i >= 0; i-- )
1378  {
1379  if ( optSum[i].set == c )
1380  break;
1381  }
1382  // store
1383  (*E)[vert]->rc.set= c;
1384  (*E)[vert]->rc.pnt= optSum[i].pnt;
1385  (*E)[vert]->rcPnt= (*pQ[c])[optSum[i].pnt];
1386  // count
1387  if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
1388 
1389 #ifdef mprDEBUG_PROT
1390  Print("\n Point E[%d] was <",vert);print_exp((*E)[vert],E->dim-1);Print(">, bucket={");
1391  for ( j= 0; j < E->dim; j++ )
1392  {
1393  Print(" %d",bucket[j]);
1394  }
1395  PrintS(" }\n optimal Sum: Qi ");
1396  for ( j= 0; j < LP->m; j++ )
1397  {
1398  Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
1399  }
1400  Print(" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].pnt);
1401 #endif
1402 
1403  // clean up
1404  omFreeSize( (void *) optSum, (LP->m) * sizeof(struct setID) );
1405 
1407 
1408  return (int)(-LP->LiPM[1][1] * SCALEDOWN);
1409 }
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
void PrintLn()
Definition: reporter.cc:327
void compute()
#define Print
Definition: emacs.cc:83
int dim
Definition: mpr_base.cc:172
CanonicalForm cd(bCommonDen(FF))
Definition: cfModGcd.cc:4030
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
#define omAlloc(size)
Definition: omAllocDecl.h:210
double mprfloat
Definition: mpr_global.h:17
bool found
Definition: facFactorize.cc:56
#define ST_SPARSE_RC
Definition: mpr_global.h:75
for(int i=0;i< R->ExpL_Size;i++) Print("%09lx "
Definition: cfEzgcd.cc:66
int j
Definition: myNF.cc:70
#define SCALEDOWN
Definition: mpr_base.cc:54
int * iposv
Definition: mpr_numeric.h:203
P bucket
Definition: myNF.cc:79
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
mprfloat ** LiPM
Definition: mpr_numeric.h:205
int num
Definition: mpr_base.cc:170
#define MAXVARS
Definition: mpr_base.cc:58
int linPolyS
Definition: mpr_base.h:47
int set
Definition: mpr_base.cc:138
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
Definition: mpr_base.cc:1221
int icase
Definition: mpr_numeric.h:201
int pnt
Definition: mpr_base.cc:139
void Werror(const char *fmt,...)
Definition: reporter.cc:199
simplex * LP
Definition: mpr_base.cc:127
bool resMatrixSparse::remapXiToPoint ( const int  indx,
pointSet **  pQ,
int *  set,
int *  vtx 
)
private

Definition at line 1221 of file mpr_base.cc.

1222 {
1223  int i,nn= (currRing->N);
1224  int loffset= 0;
1225  for ( i= 0; i <= nn; i++ )
1226  {
1227  if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) )
1228  {
1229  *set= i;
1230  *pnt= indx-loffset;
1231  return true;
1232  }
1233  else loffset+= pQ[i]->num;
1234  }
1235  return false;
1236 }
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
int i
Definition: cfEzgcd.cc:123
int num
Definition: mpr_base.cc:170
int pnt
Definition: mpr_base.cc:139

Field Documentation

ideal resMatrixSparse::gls
private

Definition at line 117 of file mpr_base.cc.

int resMatrixSparse::idelem
private

Definition at line 119 of file mpr_base.cc.

simplex* resMatrixSparse::LP
private

Definition at line 127 of file mpr_base.cc.

int resMatrixSparse::msize
private

Definition at line 121 of file mpr_base.cc.

int resMatrixSparse::n
private

Definition at line 119 of file mpr_base.cc.

int resMatrixSparse::numSet0
private

Definition at line 120 of file mpr_base.cc.

ideal resMatrixSparse::rmat
private

Definition at line 125 of file mpr_base.cc.

intvec* resMatrixSparse::uRPos
private

Definition at line 123 of file mpr_base.cc.


The documentation for this class was generated from the following file: