My Project  debian-1:4.1.2-p1+ds-2
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 66 of file mpr_base.cc.

Constructor & Destructor Documentation

◆ resMatrixSparse() [1/2]

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

Definition at line 1571 of file mpr_base.cc.

1572  : resMatrixBase(), gls( _gls )
1573 {
1574  pointSet **Qi; // vertices sets of Conv(Supp(f_i)), i=0..idelem
1575  pointSet *E; // all integer lattice points of the minkowski sum of Q0...Qn
1576  int i,k;
1577  int pnt;
1578  int totverts; // total number of exponent vectors in ideal gls
1579  mprfloat shift[MAXVARS+2]; // shiftvector delta, index [1..dim]
1580 
1581  if ( (currRing->N) > MAXVARS )
1582  {
1583  WerrorS("resMatrixSparse::resMatrixSparse: Too many variables!");
1584  return;
1585  }
1586 
1587  rmat= NULL;
1588  numSet0= 0;
1589 
1590  if ( special == SNONE ) linPolyS= 0;
1591  else linPolyS= special;
1592 
1594 
1595  n= (currRing->N);
1596  idelem= IDELEMS(gls); // should be n+1
1597 
1598  // prepare matrix LP->LiPM for Linear Programming
1599  totverts = 0;
1600  for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
1601 
1602  LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
1603 
1604  // get shift vector
1605 #ifdef mprTEST
1606  shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
1607  shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
1608 #else
1609  randomVector( idelem, shift );
1610 #endif
1611 #ifdef mprDEBUG_PROT
1612  PrintS(" shift vector: ");
1613  for ( i= 1; i <= idelem; i++ ) Print(" %.12f ",(double)shift[i]);
1614  PrintLn();
1615 #endif
1616 
1617  // evaluate convex hull for supports of gls
1618  convexHull chnp( LP );
1619  Qi= chnp.newtonPolytopesP( gls );
1620 
1621 #ifdef mprMINKSUM
1622  E= minkSumAll( Qi, n+1, n);
1623 #else
1624  // get inner points
1625  mayanPyramidAlg mpa( LP );
1626  E= mpa.getInnerPoints( Qi, shift );
1627 #endif
1628 
1629 #ifdef mprDEBUG_PROT
1630 #ifdef mprMINKSUM
1631  PrintS("(MinkSum)");
1632 #endif
1633  PrintS("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
1634  for ( pnt= 1; pnt <= E->num; pnt++ )
1635  {
1636  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1637  }
1638  PrintLn();
1639 #endif
1640 
1641 #ifdef mprTEST
1642  int lift[5][5];
1643  lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
1644  lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
1645  lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
1646  lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
1647  lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
1648  // now lift everything
1649  for ( i= 0; i <= n; i++ ) Qi[i]->lift( lift[i] );
1650 #else
1651  // now lift everything
1652  for ( i= 0; i <= n; i++ ) Qi[i]->lift();
1653 #endif
1654  E->dim++;
1655 
1656  // run Row Content Function for every point in E
1657  for ( pnt= 1; pnt <= E->num; pnt++ )
1658  {
1659  RC( Qi, E, pnt, shift );
1660  }
1661 
1662  // remove points not in cells
1663  k= E->num;
1664  for ( pnt= k; pnt > 0; pnt-- )
1665  {
1666  if ( (*E)[pnt]->rcPnt == NULL )
1667  {
1668  E->removePoint(pnt);
1670  }
1671  }
1672  mprSTICKYPROT("\n");
1673 
1674 #ifdef mprDEBUG_PROT
1675  PrintS(" points which lie in a cell:\n");
1676  for ( pnt= 1; pnt <= E->num; pnt++ )
1677  {
1678  Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
1679  }
1680  PrintLn();
1681 #endif
1682 
1683  // unlift to old dimension, sort
1684  for ( i= 0; i <= n; i++ ) Qi[i]->unlift();
1685  E->unlift();
1686  E->sort();
1687 
1688 #ifdef mprDEBUG_PROT
1689  Print(" points with a[ij] (%d):\n",E->num);
1690  for ( pnt= 1; pnt <= E->num; pnt++ )
1691  {
1692  Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
1693  Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
1694  //print_exp( (Qi[(*E)[pnt]->rc.set])[(*E)[pnt]->rc.pnt], E->dim );PrintS("> = <");
1695  print_exp( (*E)[pnt]->rcPnt, E->dim );PrintS(">\n");
1696  }
1697 #endif
1698 
1699  // now create matrix
1700  if (E->num <1)
1701  {
1702  WerrorS("could not handle a degenerate situation: no inner points found");
1703  goto theEnd;
1704  }
1705  if ( createMatrix( E ) != E->num )
1706  {
1707  // this can happen if the shiftvector shift is to large or not generic
1709  WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
1710  goto theEnd;
1711  }
1712 
1713  theEnd:
1714  // clean up
1715  for ( i= 0; i < idelem; i++ )
1716  {
1717  delete Qi[i];
1718  }
1719  omFreeSize( (void *) Qi, idelem * sizeof(pointSet*) );
1720 
1721  delete E;
1722 
1723  delete LP;
1724 }

◆ ~resMatrixSparse()

resMatrixSparse::~resMatrixSparse ( )

Definition at line 1728 of file mpr_base.cc.

1729 {
1730  delete uRPos;
1731  idDelete( &rmat );
1732 }

◆ resMatrixSparse() [2/2]

resMatrixSparse::resMatrixSparse ( const resMatrixSparse )
private

Member Function Documentation

◆ createMatrix()

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 1409 of file mpr_base.cc.

1410 {
1411  // sparse matrix
1412  // uRPos[i][1]: row of matrix
1413  // uRPos[i][idelem+1]: col of u(0)
1414  // uRPos[i][2..idelem]: col of u(1) .. u(n)
1415  // i= 1 .. numSet0
1416  int i,epos;
1417  int rp,cp;
1418  poly rowp,epp;
1419  poly iterp;
1420  int *epp_mon, *eexp;
1421 
1422  epp_mon= (int *)omAlloc( (n+2) * sizeof(int) );
1423  eexp= (int *)omAlloc0(((currRing->N)+1)*sizeof(int));
1424 
1425  totDeg= numSet0;
1426 
1427  mprSTICKYPROT2(" size of matrix: %d\n", E->num);
1428  mprSTICKYPROT2(" resultant deg: %d\n", numSet0);
1429 
1430  uRPos= new intvec( numSet0, pLength((gls->m)[0])+1, 0 );
1431 
1432  // sparse Matrix represented as a module where
1433  // each poly is column vector ( pSetComp(p,k) gives the row )
1434  rmat= idInit( E->num, E->num ); // cols, rank= number of rows
1435  msize= E->num;
1436 
1437  rp= 1;
1438  rowp= NULL;
1439  epp= pOne();
1440  for ( i= 1; i <= E->num; i++ )
1441  { // for every row
1442  E->getRowMP( i, epp_mon ); // compute (p-a[ij]), (i,j) = RC(p)
1443  pSetExpV( epp, epp_mon );
1444 
1445  //
1446  rowp= ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] ); // x^(p-a[ij]) * f(i)
1447 
1448  cp= 2;
1449  // get column for every monomial in rowp and store it
1450  iterp= rowp;
1451  while ( iterp!=NULL )
1452  {
1453  epos= E->getExpPos( iterp );
1454  if ( epos == 0 )
1455  {
1456  // this can happen, if the shift vektor or the lift funktions
1457  // are not generically chosen.
1458  Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
1459  i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
1460  return i;
1461  }
1462  pSetExpV(iterp,eexp);
1463  pSetComp(iterp, epos );
1464  pSetm(iterp);
1465  if ( (*E)[i]->rc.set == linPolyS )
1466  { // store coeff positions
1467  IMATELEM(*uRPos,rp,cp)= epos;
1468  cp++;
1469  }
1470  pIter( iterp );
1471  } // while
1472  if ( (*E)[i]->rc.set == linPolyS )
1473  { // store row
1474  IMATELEM(*uRPos,rp,1)= i-1;
1475  rp++;
1476  }
1477  (rmat->m)[i-1]= rowp;
1478  } // for
1479 
1480  pDelete( &epp );
1481  omFreeSize( (void *) epp_mon, (n+2) * sizeof(int) );
1482  omFreeSize( (void *) eexp, ((currRing->N)+1)*sizeof(int));
1483 
1484 #ifdef mprDEBUG_ALL
1485  if ( E->num <= 40 )
1486  {
1487  matrix mout= idModule2Matrix( idCopy(rmat) );
1488  print_matrix(mout);
1489  }
1490  for ( i= 1; i <= numSet0; i++ )
1491  {
1492  Print(" row %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
1493  }
1494  PrintS(" Sparse Matrix done\n");
1495 #endif
1496 
1497  return E->num;
1498 }

◆ getDetAt()

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 1795 of file mpr_base.cc.

1796 {
1797  int i,cp;
1798  poly pp,phelp,piter;
1799 
1800  mprPROTnl("smCallDet");
1801 
1802  for ( i= 1; i <= numSet0; i++ )
1803  {
1804  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1805  pDelete( &pp );
1806  pp= NULL;
1807  phelp= pp;
1808  piter= NULL;
1809  // u_1,..,u_n
1810  for ( cp= 2; cp <= idelem; cp++ )
1811  {
1812  if ( !nIsZero(evpoint[cp-1]) )
1813  {
1814  phelp= pOne();
1815  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1816  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1817  pSetmComp( phelp );
1818  if ( piter )
1819  {
1820  pNext(piter)= phelp;
1821  piter= phelp;
1822  }
1823  else
1824  {
1825  pp= phelp;
1826  piter= phelp;
1827  }
1828  }
1829  }
1830  // u0
1831  phelp= pOne();
1832  pSetCoeff( phelp, nCopy(evpoint[0]) );
1833  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1834  pSetmComp( phelp );
1835  pNext(piter)= phelp;
1836  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1837  }
1838 
1839  mprSTICKYPROT(ST__DET); // 1
1840 
1841  poly pres= sm_CallDet( rmat, currRing );
1842  number numres= nCopy( pGetCoeff( pres ) );
1843  pDelete( &pres );
1844 
1845  mprSTICKYPROT(ST__DET); // 2
1846 
1847  return ( numres );
1848 }

◆ getMatrix()

ideal resMatrixSparse::getMatrix ( )
virtual

Reimplemented from resMatrixBase.

Definition at line 1734 of file mpr_base.cc.

1735 {
1736  int i,/*j,*/cp;
1737  poly pp,phelp,piter,pgls;
1738 
1739  // copy original sparse res matrix
1740  ideal rmat_out= idCopy(rmat);
1741 
1742  // now fill in coeffs of f0
1743  for ( i= 1; i <= numSet0; i++ )
1744  {
1745 
1746  pgls= (gls->m)[0]; // f0
1747 
1748  // get matrix row and delete it
1749  pp= (rmat_out->m)[IMATELEM(*uRPos,i,1)];
1750  pDelete( &pp );
1751  pp= NULL;
1752  phelp= pp;
1753  piter= NULL;
1754 
1755  // u_1,..,u_k
1756  cp=2;
1757  while ( pNext(pgls)!=NULL )
1758  {
1759  phelp= pOne();
1760  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1761  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1762  pSetmComp( phelp );
1763  if ( piter!=NULL )
1764  {
1765  pNext(piter)= phelp;
1766  piter= phelp;
1767  }
1768  else
1769  {
1770  pp= phelp;
1771  piter= phelp;
1772  }
1773  cp++;
1774  pIter( pgls );
1775  }
1776  // u0, now pgls points to last monom
1777  phelp= pOne();
1778  pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
1779  //pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1780  pSetComp( phelp, IMATELEM(*uRPos,i,pLength((gls->m)[0])+1) );
1781  pSetmComp( phelp );
1782  if (piter!=NULL) pNext(piter)= phelp;
1783  else pp= phelp;
1784  (rmat_out->m)[IMATELEM(*uRPos,i,1)]= pp;
1785  }
1786 
1787  return rmat_out;
1788 }

◆ getUDet()

poly resMatrixSparse::getUDet ( const number *  evpoint)
virtual

Reimplemented from resMatrixBase.

Definition at line 1855 of file mpr_base.cc.

1856 {
1857  int i,cp;
1858  poly pp,phelp/*,piter*/;
1859 
1860  mprPROTnl("smCallDet");
1861 
1862  for ( i= 1; i <= numSet0; i++ )
1863  {
1864  pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
1865  pDelete( &pp );
1866  phelp= NULL;
1867  // piter= NULL;
1868  for ( cp= 2; cp <= idelem; cp++ )
1869  { // u1 .. un
1870  if ( !nIsZero(evpoint[cp-1]) )
1871  {
1872  phelp= pOne();
1873  pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
1874  pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
1875  //pSetmComp( phelp );
1876  pSetm( phelp );
1877  //Print("comp %d\n",IMATELEM(*uRPos,i,cp));
1878  #if 0
1879  if ( piter!=NULL )
1880  {
1881  pNext(piter)= phelp;
1882  piter= phelp;
1883  }
1884  else
1885  {
1886  pp= phelp;
1887  piter= phelp;
1888  }
1889  #else
1890  pp=pAdd(pp,phelp);
1891  #endif
1892  }
1893  }
1894  // u0
1895  phelp= pOne();
1896  pSetExp(phelp,1,1);
1897  pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
1898  // Print("comp %d\n",IMATELEM(*uRPos,i,idelem+1));
1899  pSetm( phelp );
1900  #if 0
1901  pNext(piter)= phelp;
1902  #else
1903  pp=pAdd(pp,phelp);
1904  #endif
1905  pTest(pp);
1906  (rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
1907  }
1908 
1909  mprSTICKYPROT(ST__DET); // 1
1910 
1911  poly pres= sm_CallDet( rmat, currRing );
1912 
1913  mprSTICKYPROT(ST__DET); // 2
1914 
1915  return ( pres );
1916 }

◆ minkSumAll()

pointSet * resMatrixSparse::minkSumAll ( pointSet **  pQ,
int  numq,
int  dim 
)
private

Definition at line 1549 of file mpr_base.cc.

1550 {
1551  pointSet *vs,*vs_old;
1552  int j;
1553 
1554  vs= new pointSet( dim );
1555 
1556  for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
1557 
1558  for ( j= 1; j < numq; j++ )
1559  {
1560  vs_old= vs;
1561  vs= minkSumTwo( vs_old, pQ[j], dim );
1562 
1563  delete vs_old;
1564  }
1565 
1566  return vs;
1567 }

◆ minkSumTwo()

pointSet * resMatrixSparse::minkSumTwo ( pointSet Q1,
pointSet Q2,
int  dim 
)
private

Definition at line 1521 of file mpr_base.cc.

1522 {
1523  pointSet *vs;
1524  onePoint vert;
1525  int j,k,l;
1526 
1527  vert.point=(Coord_t*)omAlloc( ((currRing->N)+2) * sizeof(Coord_t) );
1528 
1529  vs= new pointSet( dim );
1530 
1531  for ( j= 1; j <= Q1->num; j++ )
1532  {
1533  for ( k= 1; k <= Q2->num; k++ )
1534  {
1535  for ( l= 1; l <= dim; l++ )
1536  {
1537  vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
1538  }
1539  vs->mergeWithExp( &vert );
1540  //vs->addPoint( &vert );
1541  }
1542  }
1543 
1544  omFreeSize( (void *) vert.point, ((currRing->N)+2) * sizeof(Coord_t) );
1545 
1546  return vs;
1547 }

◆ randomVector()

void resMatrixSparse::randomVector ( const int  dim,
mprfloat  shift[] 
)
private

Definition at line 1501 of file mpr_base.cc.

1502 {
1503  int i,j;
1504  i= 1;
1505 
1506  while ( i <= dim )
1507  {
1508  shift[i]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
1509  i++;
1510  for ( j= 1; j < i-1; j++ )
1511  {
1512  if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
1513  {
1514  i--;
1515  break;
1516  }
1517  }
1518  }
1519 }

◆ RC()

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 1237 of file mpr_base.cc.

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

◆ remapXiToPoint()

bool resMatrixSparse::remapXiToPoint ( const int  indx,
pointSet **  pQ,
int *  set,
int *  vtx 
)
private

Definition at line 1218 of file mpr_base.cc.

1219 {
1220  int i,nn= (currRing->N);
1221  int loffset= 0;
1222  for ( i= 0; i <= nn; i++ )
1223  {
1224  if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) )
1225  {
1226  *set= i;
1227  *pnt= indx-loffset;
1228  return true;
1229  }
1230  else loffset+= pQ[i]->num;
1231  }
1232  return false;
1233 }

Field Documentation

◆ gls

ideal resMatrixSparse::gls
private

Definition at line 114 of file mpr_base.cc.

◆ idelem

int resMatrixSparse::idelem
private

Definition at line 116 of file mpr_base.cc.

◆ LP

simplex* resMatrixSparse::LP
private

Definition at line 124 of file mpr_base.cc.

◆ msize

int resMatrixSparse::msize
private

Definition at line 118 of file mpr_base.cc.

◆ n

int resMatrixSparse::n
private

Definition at line 116 of file mpr_base.cc.

◆ numSet0

int resMatrixSparse::numSet0
private

Definition at line 117 of file mpr_base.cc.

◆ rmat

ideal resMatrixSparse::rmat
private

Definition at line 122 of file mpr_base.cc.

◆ uRPos

intvec* resMatrixSparse::uRPos
private

Definition at line 120 of file mpr_base.cc.


The documentation for this class was generated from the following file:
dim
int dim(ideal I, ring r)
Definition: tropicalStrategy.cc:22
idCopy
ideal idCopy(ideal A)
Definition: ideals.h:59
resMatrixBase::resMatrixBase
resMatrixBase()
Definition: mpr_base.h:27
ip_smatrix
Definition: matpol.h:13
simplex::m
int m
Definition: mpr_numeric.h:197
j
int j
Definition: facHensel.cc:105
resMatrixSparse::n
int n
Definition: mpr_base.cc:116
k
int k
Definition: cfEzgcd.cc:92
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:28
resMatrixSparse::createMatrix
int createMatrix(pointSet *E)
create coeff matrix uRPos[i][1]: row of matrix uRPos[i][idelem+1]: col of u(0) uRPos[i][2....
Definition: mpr_base.cc:1409
onePoint
Definition: mpr_base.cc:139
ST_SPARSE_RCRJ
#define ST_SPARSE_RCRJ
Definition: mpr_global.h:75
lift
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
Definition: lift.cc:25
simplex::LiPM
mprfloat ** LiPM
Definition: mpr_numeric.h:204
setID
Definition: mpr_base.cc:133
resMatrixSparse::minkSumAll
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
Definition: mpr_base.cc:1549
E
REvaluation E(1, terms.length(), IntRandom(25))
ST_SPARSE_RC
#define ST_SPARSE_RC
Definition: mpr_global.h:74
pDelete
#define pDelete(p_ptr)
Definition: polys.h:175
pointSet
Definition: mpr_base.cc:160
pSetComp
#define pSetComp(p, v)
Definition: polys.h:37
mprPROTnl
#define mprPROTnl(msg)
Definition: mpr_global.h:41
ppMult_qq
#define ppMult_qq(p, q)
Definition: polys.h:196
mprSTICKYPROT
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:53
resMatrixSparse::randomVector
void randomVector(const int dim, mprfloat shift[])
Definition: mpr_base.cc:1501
found
bool found
Definition: facFactorize.cc:56
mprSTICKYPROT2
#define mprSTICKYPROT2(msg, arg)
Definition: mpr_global.h:54
resMatrixSparse::idelem
int idelem
Definition: mpr_base.cc:116
pLength
static unsigned pLength(poly a)
Definition: p_polys.h:182
for
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:65
i
int i
Definition: cfEzgcd.cc:125
cd
CanonicalForm cd(bCommonDen(FF))
Definition: cfModGcd.cc:4030
PrintS
void PrintS(const char *s)
Definition: reporter.cc:283
omFreeSize
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:258
pTest
#define pTest(p)
Definition: polys.h:398
resMatrixBase::totDeg
int totDeg
Definition: mpr_base.h:49
simplex::icase
int icase
Definition: mpr_numeric.h:200
SNONE
#define SNONE
Definition: mpr_base.h:13
pointSet::addPoint
bool addPoint(const onePointP vert)
Adds a point to pointSet, copy vert[0,...,dim] ot point[num+1][0,...,dim].
Definition: mpr_base.cc:464
resMatrixBase::istate
IStateType istate
Definition: mpr_base.h:43
convexHull
Definition: mpr_base.cc:249
resMatrixSparse::gls
ideal gls
Definition: mpr_base.cc:114
resMatrixSparse::msize
int msize
Definition: mpr_base.cc:118
size
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
pOne
#define pOne()
Definition: polys.h:299
intvec
Definition: intvec.h:18
pIter
#define pIter(p)
Definition: monomials.h:34
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:208
pointSet::dim
int dim
Definition: mpr_base.cc:169
ST__DET
#define ST__DET
Definition: mpr_global.h:77
pp
CanonicalForm pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:248
Coord_t
unsigned int Coord_t
Definition: mpr_base.cc:131
SCALEDOWN
#define SCALEDOWN
Definition: mpr_base.cc:51
MAXRVVAL
#define MAXRVVAL
Definition: mpr_base.cc:54
pAdd
#define pAdd(p, q)
Definition: polys.h:191
resMatrixSparse::numSet0
int numSet0
Definition: mpr_base.cc:117
MAXVARS
#define MAXVARS
Definition: mpr_base.cc:55
simplex
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:193
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:30
setID::set
int set
Definition: mpr_base.cc:135
nIsZero
#define nIsZero(n)
Definition: numbers.h:18
IMATELEM
#define IMATELEM(M, I, J)
Definition: intvec.h:85
resMatrixSparse::uRPos
intvec * uRPos
Definition: mpr_base.cc:120
pointSet::mergeWithExp
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet \cap point = \emptyset.
Definition: mpr_base.cc:512
resMatrixSparse::remapXiToPoint
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
Definition: mpr_base.cc:1218
resMatrixSparse::LP
simplex * LP
Definition: mpr_base.cc:124
simplex::m3
int m3
Definition: mpr_numeric.h:199
resMatrixBase::linPolyS
int linPolyS
Definition: mpr_base.h:46
RVMULT
#define RVMULT
Definition: mpr_base.cc:53
mprfloat
double mprfloat
Definition: mpr_global.h:16
pSetmComp
#define pSetmComp(p)
TODO:
Definition: polys.h:258
phelp
#define phelp
Definition: libparse.cc:1241
resMatrixBase::ready
Definition: mpr_base.h:25
Print
#define Print
Definition: emacs.cc:79
mayanPyramidAlg
Definition: mpr_base.cc:277
resMatrixBase::fatalError
Definition: mpr_base.h:25
Werror
void Werror(const char *fmt,...)
Definition: reporter.cc:188
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:34
resMatrixSparse::RC
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:1237
WerrorS
void WerrorS(const char *s)
Definition: feFopen.cc:24
onePoint::point
Coord_t * point
Definition: mpr_base.cc:141
NULL
#define NULL
Definition: omList.c:11
sm_CallDet
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:301
pSetm
#define pSetm(p)
Definition: polys.h:256
pSetExpV
#define pSetExpV(p, e)
Definition: polys.h:94
SIMPLEX_EPS
#define SIMPLEX_EPS
Definition: mpr_numeric.h:180
l
int l
Definition: cfEzgcd.cc:93
simplex::compute
void compute()
Definition: mpr_numeric.cc:1093
pointSet::num
int num
Definition: mpr_base.cc:167
pSetExp
#define pSetExp(p, i, v)
Definition: polys.h:41
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
resMatrixSparse::rmat
ideal rmat
Definition: mpr_base.cc:122
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:23
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:41
PrintLn
void PrintLn()
Definition: reporter.cc:309
simplex::n
int n
Definition: mpr_numeric.h:198
siRand
int siRand()
Definition: sirandom.c:42
resMatrixSparse::minkSumTwo
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
Definition: mpr_base.cc:1521
nCopy
#define nCopy(n)
Definition: numbers.h:14
setID::pnt
int pnt
Definition: mpr_base.cc:136
pNext
#define pNext(p)
Definition: monomials.h:33
omAlloc0
#define omAlloc0(size)
Definition: omAllocDecl.h:209
simplex::iposv
int * iposv
Definition: mpr_numeric.h:202