Loading...
Searching...
No Matches
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) .
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.
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) .
pointSet * | minkSumAll (pointSet **pQ, int numq, int dim)
pointSet * | minkSumTwo (pointSet *Q1, pointSet *Q2, int dim)

## Private Attributes

ideal | gls
int | n
int | idelem
int | numSet0
int | msize
intvec * | uRPos
ideal | rmat
simplex * | LP

## 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.

: resMatrixBase(), gls( _gls )
{
pointSet **Qi; // vertices sets of Conv(Supp(f_i)), i=0..idelem
pointSet *E; // all integer lattice points of the minkowski sum of Q0...Qn
int i,k;
int pnt;
int totverts; // total number of exponent vectors in ideal gls
mprfloat shift[MAXVARS+2]; // shiftvector delta, index [1..dim]
if ( (currRing->N) > MAXVARS )
{
WerrorS("resMatrixSparse::resMatrixSparse: Too many variables!");
return;
}
rmat= NULL;
numSet0= 0;
if ( special == SNONE ) linPolyS= 0;
else linPolyS= special;
istate= resMatrixBase::ready;
n= (currRing->N);
idelem= IDELEMS(gls); // should be n+1
// prepare matrix LP->LiPM for Linear Programming
totverts = 0;
for( i=0; i < idelem; i++) totverts += pLength( (gls->m)[i] );
LP = new simplex( idelem+totverts*2+5, totverts+5 ); // rows, cols
// get shift vector
#ifdef mprTEST
shift[0]=0.005; shift[1]=0.003; shift[2]=0.008; shift[3]=0.005; shift[4]=0.002;
shift[5]=0.1; shift[6]=0.3; shift[7]=0.2; shift[8]=0.4; shift[9]=0.2;
#else
randomVector( idelem, shift );
#endif
#ifdef mprDEBUG_PROT
PrintS(" shift vector: ");
for ( i= 1; i <= idelem; i++ ) Print(" %.12f ",(double)shift[i]);
PrintLn();
#endif
// evaluate convex hull for supports of gls
convexHull chnp( LP );
Qi= chnp.newtonPolytopesP( gls );
#ifdef mprMINKSUM
E= minkSumAll( Qi, n+1, n);
#else
// get inner points
mayanPyramidAlg mpa( LP );
E= mpa.getInnerPoints( Qi, shift );
#endif
#ifdef mprDEBUG_PROT
#ifdef mprMINKSUM
PrintS("(MinkSum)");
#endif
PrintS("\n E = (Q_0 + ... + Q_n) \\cap \\N :\n");
for ( pnt= 1; pnt <= E->num; pnt++ )
{
Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
}
PrintLn();
#endif
#ifdef mprTEST
int lift[5][5];
lift[0][1]=3; lift[0][2]=4; lift[0][3]=8; lift[0][4]=2;
lift[1][1]=6; lift[1][2]=1; lift[1][3]=7; lift[1][4]=4;
lift[2][1]=2; lift[2][2]=5; lift[2][3]=9; lift[2][4]=6;
lift[3][1]=2; lift[3][2]=1; lift[3][3]=9; lift[3][4]=5;
lift[4][1]=3; lift[4][2]=7; lift[4][3]=1; lift[4][4]=5;
// now lift everything
for ( i= 0; i <= n; i++ ) Qi[i]->lift( lift[i] );
#else
// now lift everything
for ( i= 0; i <= n; i++ ) Qi[i]->lift();
#endif
E->dim++;
// run Row Content Function for every point in E
for ( pnt= 1; pnt <= E->num; pnt++ )
{
RC( Qi, E, pnt, shift );
}
// remove points not in cells
k= E->num;
for ( pnt= k; pnt > 0; pnt-- )
{
if ( (*E)[pnt]->rcPnt == NULL )
{
E->removePoint(pnt);
mprSTICKYPROT(ST_SPARSE_RCRJ);
}
}
mprSTICKYPROT("\n");
#ifdef mprDEBUG_PROT
PrintS(" points which lie in a cell:\n");
for ( pnt= 1; pnt <= E->num; pnt++ )
{
Print("%d: <",pnt);print_exp( (*E)[pnt], E->dim );PrintS(">\n");
}
PrintLn();
#endif
// unlift to old dimension, sort
for ( i= 0; i <= n; i++ ) Qi[i]->unlift();
E->unlift();
E->sort();
#ifdef mprDEBUG_PROT
Print(" points with a[ij] (%d):\n",E->num);
for ( pnt= 1; pnt <= E->num; pnt++ )
{
Print("Punkt p \\in E[%d]: <",pnt);print_exp( (*E)[pnt], E->dim );
Print(">, RC(p) = (i:%d, j:%d), a[i,j] = <",(*E)[pnt]->rc.set,(*E)[pnt]->rc.pnt);
//print_exp( (Qi[(*E)[pnt]->rc.set])[(*E)[pnt]->rc.pnt], E->dim );PrintS("> = <");
print_exp( (*E)[pnt]->rcPnt, E->dim );PrintS(">\n");
}
#endif
// now create matrix
if (E->num <1)
{
WerrorS("could not handle a degenerate situation: no inner points found");
goto theEnd;
}
if ( createMatrix( E ) != E->num )
{
// this can happen if the shiftvector shift is to large or not generic
istate= resMatrixBase::fatalError;
WerrorS("resMatrixSparse::resMatrixSparse: Error in resMatrixSparse::createMatrix!");
goto theEnd;
}
theEnd:
// clean up
for ( i= 0; i < idelem; i++ )
{
delete Qi[i];
}
omFreeSize( (void *) Qi, idelem * sizeof(pointSet*) );
delete E;
delete LP;
}
convexHull
BOOLEAN convexHull(leftv res, leftv args)
Definition bbcone.cc:1366
i
int i
Definition cfEzgcd.cc:132
k
int k
Definition cfEzgcd.cc:99
resMatrixBase::resMatrixBase
resMatrixBase()
Definition mpr_base.h:28
resMatrixBase::linPolyS
int linPolyS
Definition mpr_base.h:47
resMatrixBase::ready
@ ready
Definition mpr_base.h:26
resMatrixBase::fatalError
@ fatalError
Definition mpr_base.h:26
resMatrixBase::istate
IStateType istate
Definition mpr_base.h:44
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
resMatrixSparse::numSet0
int numSet0
Definition mpr_base.cc:117
resMatrixSparse::rmat
ideal rmat
Definition mpr_base.cc:122
resMatrixSparse::randomVector
void randomVector(const int dim, mprfloat shift[])
Definition mpr_base.cc:1501
resMatrixSparse::LP
simplex * LP
Definition mpr_base.cc:124
resMatrixSparse::gls
ideal gls
Definition mpr_base.cc:114
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
resMatrixSparse::idelem
int idelem
Definition mpr_base.cc:116
resMatrixSparse::n
int n
Definition mpr_base.cc:116
resMatrixSparse::minkSumAll
pointSet * minkSumAll(pointSet **pQ, int numq, int dim)
Definition mpr_base.cc:1549
Print
#define Print
Definition emacs.cc:80
E
REvaluation E(1, terms.length(), IntRandom(25))
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
lift
ideal lift(const ideal J, const ring r, const ideal inI, const ring s)
Definition lift.cc:26
MAXVARS
#define MAXVARS
Definition mpr_base.cc:55
SNONE
#define SNONE
Definition mpr_base.h:14
mprSTICKYPROT
#define mprSTICKYPROT(msg)
Definition mpr_global.h:54
ST_SPARSE_RCRJ
#define ST_SPARSE_RCRJ
Definition mpr_global.h:76
mprfloat
double mprfloat
Definition mpr_global.h:17
omFreeSize
#define omFreeSize(addr, size)
Definition omAllocDecl.h:260
NULL
#define NULL
Definition omList.c:12
pLength
static int pLength(poly a)
Definition p_polys.h:190
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
PrintS
void PrintS(const char *s)
Definition reporter.cc:288
PrintLn
void PrintLn()
Definition reporter.cc:314
IDELEMS
#define IDELEMS(i)
Definition simpleideals.h:23

## ~resMatrixSparse()

resMatrixSparse::~resMatrixSparse | ( |  | ) |

Definition at line 1728 of file mpr_base.cc.

{
delete uRPos;
idDelete( &rmat );
}
resMatrixSparse::uRPos
intvec * uRPos
Definition mpr_base.cc:120
idDelete
#define idDelete(H)
delete an ideal
Definition ideals.h:29

## 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.

{
// sparse 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
int i,epos;
int rp,cp;
poly rowp,epp;
poly iterp;
int *epp_mon, *eexp;
epp_mon= (int *)omAlloc( (n+2) * sizeof(int) );
eexp= (int *)omAlloc0(((currRing->N)+1)*sizeof(int));
totDeg= numSet0;
mprSTICKYPROT2(" size of matrix: %d\n", E->num);
mprSTICKYPROT2(" resultant deg: %d\n", numSet0);
uRPos= new intvec( numSet0, pLength((gls->m)[0])+1, 0 );
// sparse Matrix represented as a module where
// each poly is column vector ( pSetComp(p,k) gives the row )
rmat= idInit( E->num, E->num ); // cols, rank= number of rows
msize= E->num;
rp= 1;
rowp= NULL;
epp= pOne();
for ( i= 1; i <= E->num; i++ )
{ // for every row
E->getRowMP( i, epp_mon ); // compute (p-a[ij]), (i,j) = RC(p)
pSetExpV( epp, epp_mon );
//
rowp= ppMult_qq( epp, (gls->m)[(*E)[i]->rc.set] ); // x^(p-a[ij]) * f(i)
cp= 2;
// get column for every monomial in rowp and store it
iterp= rowp;
while ( iterp!=NULL )
{
epos= E->getExpPos( iterp );
if ( epos == 0 )
{
// this can happen, if the shift vector or the lift functions
// are not generically chosen.
Werror("resMatrixSparse::createMatrix: Found exponent not in E, id %d, set [%d, %d]!",
i,(*E)[i]->rc.set,(*E)[i]->rc.pnt);
return i;
}
pSetExpV(iterp,eexp);
pSetComp(iterp, epos );
pSetm(iterp);
if ( (*E)[i]->rc.set == linPolyS )
{ // store coeff positions
IMATELEM(*uRPos,rp,cp)= epos;
cp++;
}
pIter( iterp );
} // while
if ( (*E)[i]->rc.set == linPolyS )
{ // store row
IMATELEM(*uRPos,rp,1)= i-1;
rp++;
}
(rmat->m)[i-1]= rowp;
} // for
pDelete( &epp );
omFreeSize( (void *) epp_mon, (n+2) * sizeof(int) );
omFreeSize( (void *) eexp, ((currRing->N)+1)*sizeof(int));
#ifdef mprDEBUG_ALL
if ( E->num <= 40 )
{
matrix mout= idModule2Matrix( idCopy(rmat) );
print_matrix(mout);
}
for ( i= 1; i <= numSet0; i++ )
{
Print(" row %d contains coeffs of f_%d\n",IMATELEM(*uRPos,i,1),linPolyS);
}
PrintS(" Sparse Matrix done\n");
#endif
return E->num;
}
resMatrixBase::totDeg
int totDeg
Definition mpr_base.h:50
resMatrixSparse::msize
int msize
Definition mpr_base.cc:118
idCopy
ideal idCopy(ideal A)
Definition ideals.h:60
IMATELEM
#define IMATELEM(M, I, J)
Definition intvec.h:86
matrix
ip_smatrix * matrix
Definition matpol.h:43
pIter
#define pIter(p)
Definition monomials.h:37
mprSTICKYPROT2
#define mprSTICKYPROT2(msg, arg)
Definition mpr_global.h:55
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
omAlloc0
#define omAlloc0(size)
Definition omAllocDecl.h:211
pDelete
#define pDelete(p_ptr)
Definition polys.h:187
pSetm
#define pSetm(p)
Definition polys.h:272
ppMult_qq
#define ppMult_qq(p, q)
Definition polys.h:209
pSetExpV
#define pSetExpV(p, e)
Definition polys.h:98
pSetComp
#define pSetComp(p, v)
Definition polys.h:39
pOne
#define pOne()
Definition polys.h:316
Werror
void Werror(const char *fmt,...)
Definition reporter.cc:189
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition simpleideals.cc:35

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

{
int i,cp;
poly pp,phelp,piter;
mprPROTnl("smCallDet");
for ( i= 1; i <= numSet0; i++ )
{
pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
pDelete( &pp );
pp= NULL;
phelp= pp;
piter= NULL;
// u_1,..,u_n
for ( cp= 2; cp <= idelem; cp++ )
{
if ( !nIsZero(evpoint[cp-1]) )
{
phelp= pOne();
pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
pSetmComp( phelp );
if ( piter )
{
pNext(piter)= phelp;
piter= phelp;
}
else
{
pp= phelp;
piter= phelp;
}
}
}
// u0
phelp= pOne();
pSetCoeff( phelp, nCopy(evpoint[0]) );
pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
pSetmComp( phelp );
pNext(piter)= phelp;
(rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
}
mprSTICKYPROT(ST__DET); // 1
poly pres= sm_CallDet( rmat, currRing );
number numres= nCopy( pGetCoeff( pres ) );
pDelete( &pres );
mprSTICKYPROT(ST__DET); // 2
return ( numres );
}
pp
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f ).
Definition cf_gcd.cc:676
phelp
#define phelp
Definition libparse.cc:1242
pNext
#define pNext(p)
Definition monomials.h:36
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:44
ST__DET
#define ST__DET
Definition mpr_global.h:78
mprPROTnl
#define mprPROTnl(msg)
Definition mpr_global.h:42
nIsZero
#define nIsZero(n)
Definition numbers.h:19
nCopy
#define nCopy(n)
Definition numbers.h:15
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition polys.h:32
pSetmComp
#define pSetmComp(p)
TODO:
Definition polys.h:274
sm_CallDet
poly sm_CallDet(ideal I, const ring R)
Definition sparsmat.cc:302

## getMatrix()

ideal resMatrixSparse::getMatrix | ( |  | ) |

virtual

Reimplemented from resMatrixBase.

Definition at line 1734 of file mpr_base.cc.

{
int i,/*j,*/cp;
poly pp,phelp,piter,pgls;
// copy original sparse res matrix
if (rmat==NULL) return NULL; //in case of error before
ideal rmat_out= idCopy(rmat);
// now fill in coeffs of f0
for ( i= 1; i <= numSet0; i++ )
{
pgls= (gls->m)[0]; // f0
// get matrix row and delete it
pp= (rmat_out->m)[IMATELEM(*uRPos,i,1)];
pDelete( &pp );
pp= NULL;
phelp= pp;
piter= NULL;
// u_1,..,u_k
cp=2;
while ( pNext(pgls)!=NULL )
{
phelp= pOne();
pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
pSetmComp( phelp );
if ( piter!=NULL )
{
pNext(piter)= phelp;
piter= phelp;
}
else
{
pp= phelp;
piter= phelp;
}
cp++;
pIter( pgls );
}
// u0, now pgls points to last monom
phelp= pOne();
pSetCoeff( phelp, nCopy(pGetCoeff(pgls)) );
//pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
pSetComp( phelp, IMATELEM(*uRPos,i,pLength((gls->m)[0])+1) );
pSetmComp( phelp );
if (piter!=NULL) pNext(piter)= phelp;
else pp= phelp;
(rmat_out->m)[IMATELEM(*uRPos,i,1)]= pp;
}
return rmat_out;
}

## getUDet()

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

virtual

Reimplemented from resMatrixBase.

Definition at line 1856 of file mpr_base.cc.

{
int i,cp;
poly pp,phelp/*,piter*/;
mprPROTnl("smCallDet");
for ( i= 1; i <= numSet0; i++ )
{
pp= (rmat->m)[IMATELEM(*uRPos,i,1)];
pDelete( &pp );
phelp= NULL;
// piter= NULL;
for ( cp= 2; cp <= idelem; cp++ )
{ // u1 .. un
if ( !nIsZero(evpoint[cp-1]) )
{
phelp= pOne();
pSetCoeff( phelp, nCopy(evpoint[cp-1]) );
pSetComp( phelp, IMATELEM(*uRPos,i,cp) );
//pSetmComp( phelp );
pSetm( phelp );
//Print("comp %d\n",IMATELEM(*uRPos,i,cp));
#if 0
if ( piter!=NULL )
{
pNext(piter)= phelp;
piter= phelp;
}
else
{
pp= phelp;
piter= phelp;
}
#else
pp=pAdd(pp,phelp);
#endif
}
}
// u0
phelp= pOne();
pSetExp(phelp,1,1);
pSetComp( phelp, IMATELEM(*uRPos,i,idelem+1) );
// Print("comp %d\n",IMATELEM(*uRPos,i,idelem+1));
pSetm( phelp );
#if 0
pNext(piter)= phelp;
#else
pp=pAdd(pp,phelp);
#endif
pTest(pp);
(rmat->m)[IMATELEM(*uRPos,i,1)]= pp;
}
mprSTICKYPROT(ST__DET); // 1
poly pres= sm_CallDet( rmat, currRing );
mprSTICKYPROT(ST__DET); // 2
return ( pres );
}
pAdd
#define pAdd(p, q)
Definition polys.h:204
pTest
#define pTest(p)
Definition polys.h:415
pSetExp
#define pSetExp(p, i, v)
Definition polys.h:43

## minkSumAll()

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

private

Definition at line 1549 of file mpr_base.cc.

{
pointSet *vs,*vs_old;
int j;
vs= new pointSet( dim );
for ( j= 1; j <= pQ[0]->num; j++ ) vs->addPoint( (*pQ[0])[j] );
for ( j= 1; j < numq; j++ )
{
vs_old= vs;
vs= minkSumTwo( vs_old, pQ[j], dim );
delete vs_old;
}
return vs;
}
pointSet::addPoint
bool addPoint(const onePointP vert)
Adds a point to pointSet, copy vert[0,...,dim] to point[num+1][0,...,dim].
Definition mpr_base.cc:464
pointSet::num
int num
Definition mpr_base.cc:167
resMatrixSparse::minkSumTwo
pointSet * minkSumTwo(pointSet *Q1, pointSet *Q2, int dim)
Definition mpr_base.cc:1521
j
int j
Definition facHensel.cc:110
dim
int dim(ideal I, ring r)
Definition tropicalStrategy.cc:23

## minkSumTwo()

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

private

Definition at line 1521 of file mpr_base.cc.

{
pointSet *vs;
onePoint vert;
int j,k,l;
vert.point=(Coord_t*)omAlloc( ((currRing->N)+2) * sizeof(Coord_t) );
vs= new pointSet( dim );
for ( j= 1; j <= Q1->num; j++ )
{
for ( k= 1; k <= Q2->num; k++ )
{
for ( l= 1; l <= dim; l++ )
{
vert.point[l]= (*Q1)[j]->point[l] + (*Q2)[k]->point[l];
}
vs->mergeWithExp( &vert );
//vs->addPoint( &vert );
}
}
omFreeSize( (void *) vert.point, ((currRing->N)+2) * sizeof(Coord_t) );
return vs;
}
l
int l
Definition cfEzgcd.cc:100
pointSet::mergeWithExp
bool mergeWithExp(const onePointP vert)
Adds point to pointSet, iff pointSet \cap point = \emptyset.
Definition mpr_base.cc:512
Coord_t
unsigned int Coord_t
Definition mpr_base.cc:131
onePoint::point
Coord_t * point
Definition mpr_base.cc:141

## randomVector()

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

private

Definition at line 1501 of file mpr_base.cc.

{
int i,j;
i= 1;
while ( i <= dim )
{
shift[i]= (mprfloat) (RVMULT*(siRand()%MAXRVVAL)/(mprfloat)MAXRVVAL);
i++;
for ( j= 1; j < i-1; j++ )
{
if ( (shift[j] < shift[i-1] + SIMPLEX_EPS) && (shift[j] > shift[i-1] - SIMPLEX_EPS) )
{
i--;
break;
}
}
}
}
RVMULT
#define RVMULT
Definition mpr_base.cc:53
MAXRVVAL
#define MAXRVVAL
Definition mpr_base.cc:54
SIMPLEX_EPS
#define SIMPLEX_EPS
Definition mpr_numeric.h:181
siRand
int siRand()
Definition sirandom.c:42

## 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.

{
int i, j, k,c ;
int size;
bool found= true;
mprfloat cd;
int onum;
int bucket[MAXVARS+2];
setID *optSum;
LP->n = 1;
LP->m = n + n + 1; // number of constrains
// fill in LP matrix
for ( i= 0; i <= n; i++ )
{
size= pQ[i]->num;
for ( k= 1; k <= size; k++ )
{
LP->n++;
// objective function, minimize
LP->LiPM[1][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[pQ[i]->dim] / SCALEDOWN );
// lambdas sum up to 1
for ( j = 0; j <= n; j++ )
{
if ( i==j )
LP->LiPM[j+2][LP->n] = -1.0;
else
LP->LiPM[j+2][LP->n] = 0.0;
}
// the points
for ( j = 1; j <= n; j++ )
{
LP->LiPM[j+n+2][LP->n] = - ( (mprfloat) (*pQ[i])[k]->point[j] );
}
}
}
for ( j = 0; j <= n; j++ ) LP->LiPM[j+2][1] = 1.0;
for ( j= 1; j <= n; j++ )
{
LP->LiPM[j+n+2][1]= (mprfloat)(*E)[vert]->point[j] - shift[j];
}
LP->n--;
LP->LiPM[1][1] = 0.0;
#ifdef mprDEBUG_ALL
PrintLn();
Print(" n= %d, LP->m=M= %d, LP->n=N= %d\n",n,LP->m,LP->n);
print_mat(LP->LiPM, LP->m+1, LP->n+1);
#endif
LP->m3= LP->m;
LP->compute();
if ( LP->icase < 0 )
{
// infeasibility: the point does not lie in a cell -> remove it
return -1;
}
// store result
(*E)[vert]->point[E->dim]= (int)(-LP->LiPM[1][1] * SCALEDOWN);
#ifdef mprDEBUG_ALL
Print(" simplx returned %d, Objective value = %f\n", LP->icase, LP->LiPM[1][1]);
//print_bmat(LP->LiPM, NumCons + 1, LP->n+1-NumCons, LP->n+1, LP->iposv); // ( rows= M+1, cols= N+1-m3 )
//print_mat(LP->LiPM, NumCons+1, LP->n);
#endif
#if 1
// sort LP results
while (found)
{
found=false;
for ( i= 1; i < LP->m; i++ )
{
if ( LP->iposv[i] > LP->iposv[i+1] )
{
c= LP->iposv[i];
LP->iposv[i]=LP->iposv[i+1];
LP->iposv[i+1]=c;
cd=LP->LiPM[i+1][1];
LP->LiPM[i+1][1]=LP->LiPM[i+2][1];
LP->LiPM[i+2][1]=cd;
found= true;
}
}
}
#endif
#ifdef mprDEBUG_ALL
print_bmat(LP->LiPM, LP->m + 1, LP->n+1-LP->m, LP->n+1, LP->iposv);
PrintS(" now split into sets\n");
#endif
// init bucket
for ( i= 0; i <= E->dim; i++ ) bucket[i]= 0;
// remap results of LP to sets Qi
c=0;
optSum= (setID*)omAlloc( (LP->m) * sizeof(struct setID) );
for ( i= 0; i < LP->m; i++ )
{
//Print("% .15f\n",LP->LiPM[i+2][1]);
if ( LP->LiPM[i+2][1] > 1e-12 )
{
if ( !remapXiToPoint( LP->iposv[i+1], pQ, &(optSum[c].set), &(optSum[c].pnt) ) )
{
Werror(" resMatrixSparse::RC: Found bad solution in LP: %d!",LP->iposv[i+1]);
WerrorS(" resMatrixSparse::RC: remapXiToPoint failed!");
return -1;
}
bucket[optSum[c].set]++;
c++;
}
}
onum= c;
// find last min in bucket[]: maximum i such that Fi is a point
c= 0;
for ( i= 1; i < E->dim; i++ )
{
if ( bucket[c] >= bucket[i] )
{
c= i;
}
}
// find matching point set
for ( i= onum - 1; i >= 0; i-- )
{
if ( optSum[i].set == c )
break;
}
// store
(*E)[vert]->rc.set= c;
(*E)[vert]->rc.pnt= optSum[i].pnt;
(*E)[vert]->rcPnt= (*pQ[c])[optSum[i].pnt];
// count
if ( (*E)[vert]->rc.set == linPolyS ) numSet0++;
#ifdef mprDEBUG_PROT
Print("\n Point E[%d] was <",vert);print_exp((*E)[vert],E->dim-1);Print(">, bucket={");
for ( j= 0; j < E->dim; j++ )
{
Print(" %d",bucket[j]);
}
PrintS(" }\n optimal Sum: Qi ");
for ( j= 0; j < LP->m; j++ )
{
Print(" [ %d, %d ]",optSum[j].set,optSum[j].pnt);
}
Print(" -> i= %d, j = %d\n",(*E)[vert]->rc.set,optSum[i].pnt);
#endif
// clean up
omFreeSize( (void *) optSum, (LP->m) * sizeof(struct setID) );
mprSTICKYPROT(ST_SPARSE_RC);
return (int)(-LP->LiPM[1][1] * SCALEDOWN);
}
size
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition cf_ops.cc:600
cd
CanonicalForm cd(bCommonDen(FF))
Definition cfModGcd.cc:4097
pointSet::dim
int dim
Definition mpr_base.cc:169
resMatrixSparse::remapXiToPoint
bool remapXiToPoint(const int indx, pointSet **pQ, int *set, int *vtx)
Definition mpr_base.cc:1218
found
bool found
Definition facFactorize.cc:55
SCALEDOWN
#define SCALEDOWN
Definition mpr_base.cc:51
setID::set
int set
Definition mpr_base.cc:135
setID::pnt
int pnt
Definition mpr_base.cc:136
ST_SPARSE_RC
#define ST_SPARSE_RC
Definition mpr_global.h:75

## remapXiToPoint()

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

private

Definition at line 1218 of file mpr_base.cc.

{
int i,nn= (currRing->N);
int loffset= 0;
for ( i= 0; i <= nn; i++ )
{
if ( (loffset < indx) && (indx <= pQ[i]->num + loffset) )
{
*set= i;
*pnt= indx-loffset;
return true;
}
else loffset+= pQ[i]->num;
}
return false;
}
num
CanonicalForm num(const CanonicalForm &f)
Definition canonicalform.h:337

## 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:

- kernel/numeric/mpr_base.cc