My Project  debian-1:4.1.1-p2+ds-4build2
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 68 of file mpr_base.cc.

Constructor & Destructor Documentation

◆ resMatrixSparse() [1/2]

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

Definition at line 1573 of file mpr_base.cc.

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

◆ ~resMatrixSparse()

resMatrixSparse::~resMatrixSparse ( )

Definition at line 1730 of file mpr_base.cc.

1731 {
1732  delete uRPos;
1733  idDelete( &rmat );
1734 }

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

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

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

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

◆ getMatrix()

ideal resMatrixSparse::getMatrix ( )
virtual

Reimplemented from resMatrixBase.

Definition at line 1736 of file mpr_base.cc.

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

◆ getUDet()

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

Reimplemented from resMatrixBase.

Definition at line 1857 of file mpr_base.cc.

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

◆ minkSumAll()

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

Definition at line 1551 of file mpr_base.cc.

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

◆ minkSumTwo()

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

Definition at line 1523 of file mpr_base.cc.

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

◆ randomVector()

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

Definition at line 1503 of file mpr_base.cc.

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

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

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

◆ remapXiToPoint()

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

Definition at line 1220 of file mpr_base.cc.

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

Field Documentation

◆ gls

ideal resMatrixSparse::gls
private

Definition at line 116 of file mpr_base.cc.

◆ idelem

int resMatrixSparse::idelem
private

Definition at line 118 of file mpr_base.cc.

◆ LP

simplex* resMatrixSparse::LP
private

Definition at line 126 of file mpr_base.cc.

◆ msize

int resMatrixSparse::msize
private

Definition at line 120 of file mpr_base.cc.

◆ n

int resMatrixSparse::n
private

Definition at line 118 of file mpr_base.cc.

◆ numSet0

int resMatrixSparse::numSet0
private

Definition at line 119 of file mpr_base.cc.

◆ rmat

ideal resMatrixSparse::rmat
private

Definition at line 124 of file mpr_base.cc.

◆ uRPos

intvec* resMatrixSparse::uRPos
private

Definition at line 122 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:23
idCopy
ideal idCopy(ideal A)
Definition: ideals.h:60
resMatrixBase::resMatrixBase
resMatrixBase()
Definition: mpr_base.h:28
ip_smatrix
Definition: matpol.h:14
simplex::m
int m
Definition: mpr_numeric.h:198
j
int j
Definition: facHensel.cc:105
resMatrixSparse::n
int n
Definition: mpr_base.cc:118
k
int k
Definition: cfEzgcd.cc:92
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
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:1411
onePoint
Definition: mpr_base.cc:141
ST_SPARSE_RCRJ
#define ST_SPARSE_RCRJ
Definition: mpr_global.h:76
lift
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
Definition: lift.cc:26
simplex::LiPM
mprfloat ** LiPM
Definition: mpr_numeric.h:205
setID
Definition: mpr_base.cc:135
resMatrixSparse::minkSumAll
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
Definition: mpr_base.cc:1551
E
REvaluation E(1, terms.length(), IntRandom(25))
ST_SPARSE_RC
#define ST_SPARSE_RC
Definition: mpr_global.h:75
pDelete
#define pDelete(p_ptr)
Definition: polys.h:173
pointSet
Definition: mpr_base.cc:162
pSetComp
#define pSetComp(p, v)
Definition: polys.h:38
mprPROTnl
#define mprPROTnl(msg)
Definition: mpr_global.h:42
ppMult_qq
#define ppMult_qq(p, q)
Definition: polys.h:195
mprSTICKYPROT
#define mprSTICKYPROT(msg)
Definition: mpr_global.h:54
resMatrixSparse::randomVector
void randomVector(const int dim, mprfloat shift[])
Definition: mpr_base.cc:1503
found
bool found
Definition: facFactorize.cc:56
mprSTICKYPROT2
#define mprSTICKYPROT2(msg, arg)
Definition: mpr_global.h:55
resMatrixSparse::idelem
int idelem
Definition: mpr_base.cc:118
pLength
static unsigned pLength(poly a)
Definition: p_polys.h:192
for
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:65
currRing
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
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:284
omFreeSize
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
pTest
#define pTest(p)
Definition: polys.h:401
resMatrixBase::totDeg
int totDeg
Definition: mpr_base.h:50
simplex::icase
int icase
Definition: mpr_numeric.h:201
SNONE
#define SNONE
Definition: mpr_base.h:14
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:466
resMatrixBase::istate
IStateType istate
Definition: mpr_base.h:44
convexHull
Definition: mpr_base.cc:251
resMatrixSparse::gls
ideal gls
Definition: mpr_base.cc:116
resMatrixSparse::msize
int msize
Definition: mpr_base.cc:120
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:301
intvec
Definition: intvec.h:17
pIter
#define pIter(p)
Definition: monomials.h:44
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:210
pointSet::dim
int dim
Definition: mpr_base.cc:171
ST__DET
#define ST__DET
Definition: mpr_global.h:78
pp
CanonicalForm pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:253
Coord_t
unsigned int Coord_t
Definition: mpr_base.cc:133
SCALEDOWN
#define SCALEDOWN
Definition: mpr_base.cc:53
MAXRVVAL
#define MAXRVVAL
Definition: mpr_base.cc:56
pAdd
#define pAdd(p, q)
Definition: polys.h:190
resMatrixSparse::numSet0
int numSet0
Definition: mpr_base.cc:119
MAXVARS
#define MAXVARS
Definition: mpr_base.cc:57
simplex
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:31
setID::set
int set
Definition: mpr_base.cc:137
nIsZero
#define nIsZero(n)
Definition: numbers.h:20
IMATELEM
#define IMATELEM(M, I, J)
Definition: intvec.h:85
resMatrixSparse::uRPos
intvec * uRPos
Definition: mpr_base.cc:122
pointSet::mergeWithExp
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet \cap point = \emptyset.
Definition: mpr_base.cc:514
resMatrixSparse::remapXiToPoint
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
Definition: mpr_base.cc:1220
resMatrixSparse::LP
simplex * LP
Definition: mpr_base.cc:126
simplex::m3
int m3
Definition: mpr_numeric.h:200
resMatrixBase::linPolyS
int linPolyS
Definition: mpr_base.h:47
RVMULT
#define RVMULT
Definition: mpr_base.cc:55
mprfloat
double mprfloat
Definition: mpr_global.h:17
pSetmComp
#define pSetmComp(p)
TODO:
Definition: polys.h:259
phelp
#define phelp
Definition: libparse.cc:1240
resMatrixBase::ready
@ ready
Definition: mpr_base.h:26
Print
#define Print
Definition: emacs.cc:80
mayanPyramidAlg
Definition: mpr_base.cc:279
resMatrixBase::fatalError
@ fatalError
Definition: mpr_base.h:26
Werror
void Werror(const char *fmt,...)
Definition: reporter.cc:189
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:37
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:1239
WerrorS
void WerrorS(const char *s)
Definition: feFopen.cc:24
onePoint::point
Coord_t * point
Definition: mpr_base.cc:143
NULL
#define NULL
Definition: omList.c:10
sm_CallDet
poly sm_CallDet(ideal I, const ring R)
Definition: sparsmat.cc:357
pSetm
#define pSetm(p)
Definition: polys.h:257
pSetExpV
#define pSetExpV(p, e)
Definition: polys.h:97
SIMPLEX_EPS
#define SIMPLEX_EPS
Definition: mpr_numeric.h:181
l
int l
Definition: cfEzgcd.cc:93
simplex::compute
void compute()
Definition: mpr_numeric.cc:1099
pointSet::num
int num
Definition: mpr_base.cc:169
pSetExp
#define pSetExp(p, i, v)
Definition: polys.h:42
resMatrixSparse::rmat
ideal rmat
Definition: mpr_base.cc:124
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:26
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:51
PrintLn
void PrintLn()
Definition: reporter.cc:310
simplex::n
int n
Definition: mpr_numeric.h:199
siRand
int siRand()
Definition: sirandom.c:41
resMatrixSparse::minkSumTwo
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
Definition: mpr_base.cc:1523
nCopy
#define nCopy(n)
Definition: numbers.h:16
setID::pnt
int pnt
Definition: mpr_base.cc:138
pNext
#define pNext(p)
Definition: monomials.h:43
omAlloc0
#define omAlloc0(size)
Definition: omAllocDecl.h:211
simplex::iposv
int * iposv
Definition: mpr_numeric.h:203