resMatrixDense Class Reference

Public Member Functions
	resMatrixDense (const ideal _gls, const int special=SNONE)
	_gls: system of multivariate polynoms special: -1 -> resMatrixDense is a symbolic matrix 0,1, ... -> resMatrixDense ist eine u-Resultante, wobei special das lineare u-Polynom angibt More...
	~resMatrixDense ()
resVector *	getMVector (const int i)
	column vector of matrix, index von 0 ... numVectors-1 More...
ideal	getMatrix ()
	Returns the matrix M in an usable presentation. More...
ideal	getSubMatrix ()
	Returns the submatrix M' of M in an usable presentation. More...
number	getDetAt (const number *evpoint)
	Evaluate the determinant of the matrix M at the point evpoint where the ui's are replaced by the components of evpoint. More...
number	getSubDet ()
	Evaluates the determinant of the submatrix M'. More...
Public Member Functions inherited from resMatrixBase
	resMatrixBase ()
virtual	~resMatrixBase ()
virtual ideal	getMatrix ()
virtual ideal	getSubMatrix ()
virtual poly	getUDet (const number *)
virtual number	getDetAt (const number *)
virtual number	getSubDet ()
virtual long	getDetDeg ()
virtual IStateType	initState () const

Private Member Functions
	resMatrixDense (const resMatrixDense &)
	deactivated copy constructor More...
void	generateBaseData ()
	Generate the "matrix" M. More...
void	generateMonomData (int deg, intvec *polyDegs, intvec *iVO)
	Generates needed set of monoms, split them into sets S0, ... Sn and check if reduced/nonreduced and calculate size of submatrix. More...
void	generateMonoms (poly m, int var, int deg)
	Recursively generate all homogeneous monoms of (currRing->N) variables of degree deg. More...
void	createMatrix ()
	Creates quadratic matrix M of size numVectors for later use. More...

Private Attributes
resVector *	resVectorList
int	veclistmax
int	veclistblock
int	numVectors
int	subSize
matrix	m

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

Constructor & Destructor Documentation

resMatrixDense() [1/2]

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

_gls: system of multivariate polynoms special: -1 -> resMatrixDense is a symbolic matrix 0,1, ... -> resMatrixDense ist eine u-Resultante, wobei special das lineare u-Polynom angibt

Definition at line 2064 of file mpr_base.cc.

  : resMatrixBase()
{
  int i;
 
  sourceRing=currRing;
  gls= idCopy( _gls );
  linPolyS= special;
  m=NULL;
 
  // init all
  generateBaseData();
 
  totDeg= 1;
  for ( i= 0; i < IDELEMS(gls); i++ )
  {
    totDeg*=pTotaldegree( (gls->m)[i] );
  }
 
  mprSTICKYPROT2("  resultant deg: %d\n",totDeg);
 
  istate= resMatrixBase::ready;
}

~resMatrixDense()

resMatrixDense::~resMatrixDense

(

)

Definition at line 2088 of file mpr_base.cc.

{
  int i,j;
  for (i=0; i < numVectors; i++)
  {
    pDelete( &resVectorList[i].mon );
    pDelete( &resVectorList[i].dividedBy );
    for ( j=0; j < resVectorList[i].numColVectorSize; j++ )
    {
        nDelete( resVectorList[i].numColVector+j );
    }
    // OB: ????? (solve_s.tst)
    if (resVectorList[i].numColVector!=NULL)
      omfreeSize( (void *)resVectorList[i].numColVector,
                numVectors * sizeof( number ) );
    if (resVectorList[i].numColParNr!=NULL)
      omfreeSize( (void *)resVectorList[i].numColParNr,
                ((currRing->N)+1) * sizeof(int) );
  }
 
  omFreeSize( (void *)resVectorList, veclistmax*sizeof( resVector ) );
 
  // free matrix m
  if ( m != NULL )
  {
    idDelete((ideal *)&m);
  }
}

resMatrixDense() [2/2]

resMatrixDense::resMatrixDense ( const resMatrixDense &  )

private

deactivated copy constructor

Member Function Documentation

createMatrix()

void resMatrixDense::createMatrix ( )

private

Creates quadratic matrix M of size numVectors for later use.

u0, u1, ...,un are replaced by 0. Entries equal to 0 are not initialized ( == NULL)

Definition at line 2120 of file mpr_base.cc.

{
  int k,i,j;
  resVector *vecp;
 
  m= mpNew( numVectors, numVectors );
 
  for ( i= 1; i <= MATROWS( m ); i++ )
    for ( j= 1; j <= MATCOLS( m ); j++ )
    {
      MATELEM(m,i,j)= pInit();
      pSetCoeff0( MATELEM(m,i,j), nInit(0) );
    }
 
 
  for ( k= 0; k <= numVectors - 1; k++ )
  {
    if ( linPolyS == getMVector(k)->elementOfS )
    {
      mprSTICKYPROT(ST_DENSE_FR);
      for ( i= 0; i < (currRing->N); i++ )
      {
        MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i])= pInit();
      }
    }
    else
    {
      mprSTICKYPROT(ST_DENSE_NR);
      vecp= getMVector(k);
      for ( i= 0; i < numVectors; i++)
      {
        if ( !nIsZero( vecp->getElemNum(i) ) )
        {
          MATELEM(m,numVectors - k,i + 1)= pInit();
          pSetCoeff0( MATELEM(m,numVectors - k,i + 1), nCopy(vecp->getElemNum(i)) );
        }
      }
    }
  } // for
  mprSTICKYPROT("\n");
 
#ifdef mprDEBUG_ALL
  for ( k= numVectors - 1; k >= 0; k-- )
  {
    if ( linPolyS == getMVector(k)->elementOfS )
    {
      for ( i=0; i < (currRing->N); i++ )
      {
        Print(" %d ",(getMVector(k)->numColParNr)[i]);
      }
      PrintLn();
    }
  }
  for (i=1; i <= numVectors; i++)
  {
    for (j=1; j <= numVectors; j++ )
    {
      pWrite0(MATELEM(m,i,j));PrintS("  ");
    }
    PrintLn();
  }
#endif
}

generateBaseData()

void resMatrixDense::generateBaseData ( )

private

Generate the "matrix" M.

Each column is presented by a resVector holding all entries for this column.

Definition at line 2343 of file mpr_base.cc.

{
  int k,j,i;
  number matEntry;
  poly pmatchPos;
  poly pi,factor,pmp;
 
  // holds the degrees of F0, F1, ..., Fn
  intvec polyDegs( IDELEMS(gls) );
  for ( k= 0; k < IDELEMS(gls); k++ )
    polyDegs[k]= pTotaldegree( (gls->m)[k] );
 
  // the internal Variable Ordering
  // make sure that the homogenization variable goes last!
  intvec iVO( (currRing->N) );
  if ( linPolyS != SNONE )
  {
    iVO[(currRing->N) - 1]= linPolyS;
    int p=0;
    for ( k= (currRing->N) - 1; k >= 0; k-- )
    {
      if ( k != linPolyS )
      {
        iVO[p]= k;
        p++;
      }
    }
  }
  else
  {
    linPolyS= 0;
    for ( k= 0; k < (currRing->N); k++ )
      iVO[k]= (currRing->N) - k - 1;
  }
 
  // the critical degree d= sum( deg(Fi) ) - n
  int sumDeg= 0;
  for ( k= 0; k < polyDegs.rows(); k++ )
    sumDeg+= polyDegs[k];
  sumDeg-= polyDegs.rows() - 1;
 
  // generate the base data
  generateMonomData( sumDeg, &polyDegs, &iVO );
 
  // generate "matrix"
  for ( k= numVectors - 1; k >= 0; k-- )
  {
    if ( resVectorList[k].elementOfS != linPolyS )
    {
      // column k is a normal column with numerical or symbolic entries
      // init stuff
      resVectorList[k].numColParNr= NULL;
      resVectorList[k].numColVectorSize= numVectors;
      resVectorList[k].numColVector= (number *)omAlloc( numVectors*sizeof( number ) );
      for ( i= 0; i < numVectors; i++ ) resVectorList[k].numColVector[i]= nInit(0);
 
      // compute row poly
      poly pi= ppMult_qq( (gls->m)[ resVectorList[k].elementOfS ] , resVectorList[k].mon );
      pi= pp_DivideM( pi, resVectorList[k].dividedBy, currRing  );
 
      // fill in "matrix"
      while ( pi != NULL )
      {
        matEntry= nCopy(pGetCoeff(pi));
        pmatchPos= pLmInit( pi );
        pSetCoeff0( pmatchPos, nInit(1) );
 
        for ( i= 0; i < numVectors; i++)
          if ( pLmEqual( pmatchPos, resVectorList[i].mon ) )
            break;
 
        resVectorList[k].numColVector[numVectors - i - 1] = nCopy(matEntry);
 
        pDelete( &pmatchPos );
        nDelete( &matEntry );
 
        pIter( pi );
      }
      pDelete( &pi );
    }
    else
    {
      // column is a special column, i.e. is generated by S0 and F0
      // safe only the positions of the ui's in the column
      //mprPROTInl(" setup of numColParNr ",k);
      resVectorList[k].numColVectorSize= 0;
      resVectorList[k].numColVector= NULL;
      resVectorList[k].numColParNr= (int *)omAlloc0( ((currRing->N)+1) * sizeof(int) );
 
      pi= (gls->m)[ resVectorList[k].elementOfS ];
      factor= pp_DivideM( resVectorList[k].mon, resVectorList[k].dividedBy, currRing );
 
      j=0;
      while ( pi  != NULL )
      { // fill in "matrix"
        pmp= pMult( pCopy( factor ), pHead( pi ) );
        pTest( pmp );
 
        for ( i= 0; i < numVectors; i++)
          if ( pLmEqual( pmp, resVectorList[i].mon ) )
            break;
 
        resVectorList[k].numColParNr[j]= i;
        pDelete( &pmp );
        pIter( pi );
        j++;
      }
      pDelete( &pi );
      pDelete( &factor );
    }
  } // for ( k= numVectors - 1; k >= 0; k-- )
 
  mprSTICKYPROT2(" size of matrix:    %d\n",numVectors);
  mprSTICKYPROT2(" size of submatrix: %d\n",subSize);
 
  // create the matrix M
  createMatrix();
 
}

generateMonomData()

void resMatrixDense::generateMonomData ( int  deg, intvec *  polyDegs, intvec *  iVO  )

private

Generates needed set of monoms, split them into sets S0, ... Sn and check if reduced/nonreduced and calculate size of submatrix.

Definition at line 2227 of file mpr_base.cc.

{
  int i,j,k;
 
  // init monomData
  veclistblock= 512;
  veclistmax= veclistblock;
  resVectorList= (resVector *)omAlloc( veclistmax*sizeof( resVector ) );
 
  // Init resVector()s
  for ( j= veclistmax - 1; j >= 0; j-- ) resVectorList[j].init();
  numVectors= 0;
 
  // Generate all monoms of degree deg
  poly start= pOne();
  generateMonoms( start, 1, deg );
  pDelete( & start );
 
  mprSTICKYPROT("\n");
 
  // Check for reduced monoms
  // First generate polyDegs.rows() monoms
  //  x(k)^(polyDegs[k]),  0 <= k < polyDegs.rows()
  ideal pDegDiv= idInit( polyDegs->rows(), 1 );
  for ( k= 0; k < polyDegs->rows(); k++ )
  {
    poly p= pOne();
    pSetExp( p, k + 1, (*polyDegs)[k] );
    pSetm( p );
    (pDegDiv->m)[k]= p;
  }
 
  // Now check each monom if it is reduced.
  // A monom monom is called reduced if there exists
  // exactly one x(k)^(polyDegs[k]) that divides the monom.
  int divCount;
  for ( j= numVectors - 1; j >= 0; j-- )
  {
    divCount= 0;
    for ( k= 0; k < IDELEMS(pDegDiv); k++ )
      if ( pLmDivisibleByNoComp( (pDegDiv->m)[k], resVectorList[j].mon ) )
        divCount++;
    resVectorList[j].isReduced= (divCount == 1);
  }
 
  // create the sets S(k)s
  // a monom x(i)^deg, deg given, is element of the set S(i)
  // if all x(0)^(polyDegs[0]) ... x(i-1)^(polyDegs[i-1]) DONT divide
  // x(i)^deg and only x(i)^(polyDegs[i]) divides x(i)^deg
  bool doInsert;
  for ( k= 0; k < iVO->rows(); k++)
  {
    //mprPROTInl(" ------------ var:",(*iVO)[k]);
    for ( j= numVectors - 1; j >= 0; j-- )
    {
      //mprPROTPnl("testing monom",resVectorList[j].mon);
      if ( resVectorList[j].elementOfS == SFREE )
      {
        //mprPROTnl("\tfree");
        if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[k] ], resVectorList[j].mon ) )
        {
          //mprPROTPnl("\tdivisible by ",(pDegDiv->m)[ (*iVO)[k] ]);
          doInsert=TRUE;
          for ( i= 0; i < k; i++ )
          {
            //mprPROTPnl("\tchecking db ",(pDegDiv->m)[ (*iVO)[i] ]);
            if ( pLmDivisibleByNoComp( (pDegDiv->m)[ (*iVO)[i] ], resVectorList[j].mon ) )
            {
              //mprPROTPnl("\t and divisible by",(pDegDiv->m)[ (*iVO)[i] ]);
              doInsert=FALSE;
              break;
            }
          }
          if ( doInsert )
          {
            //mprPROTInl("\t------------------> S ",(*iVO)[k]);
            resVectorList[j].elementOfS= (*iVO)[k];
            resVectorList[j].dividedBy= pCopy( (pDegDiv->m)[ (*iVO)[i] ] );
          }
        }
      }
    }
  }
 
  // size of submatrix M', equal to number of nonreduced monoms
  // (size of matrix M is equal to number of monoms=numVectors)
  subSize= 0;
  int sub;
  for ( i= 0; i < polyDegs->rows(); i++ )
  {
    sub= 1;
    for ( k= 0; k < polyDegs->rows(); k++ )
      if ( i != k ) sub*= (*polyDegs)[k];
    subSize+= sub;
  }
  subSize= numVectors - subSize;
 
  // pDegDiv wieder freigeben!
  idDelete( &pDegDiv );
 
#ifdef mprDEBUG_ALL
  // Print a list of monoms and their properties
  PrintS("// \n");
  for ( j= numVectors - 1; j >= 0; j-- )
  {
    Print("// %s, S(%d),  db ",
          resVectorList[j].isReduced?"reduced":"nonreduced",
          resVectorList[j].elementOfS);
    pWrite0(resVectorList[j].dividedBy);
    PrintS("  monom ");
    pWrite(resVectorList[j].mon);
  }
  Print("// size: %d, subSize: %d\n",numVectors,subSize);
#endif
}

generateMonoms()

void resMatrixDense::generateMonoms ( poly  m, int  var, int  deg  )

private

Recursively generate all homogeneous monoms of (currRing->N) variables of degree deg.

Definition at line 2187 of file mpr_base.cc.

{
  if ( deg == 0 )
  {
    poly mon = pCopy( mm );
 
    if ( numVectors == veclistmax )
    {
      resVectorList= (resVector * )omReallocSize( resVectorList,
                                            (veclistmax) * sizeof( resVector ),
                                            (veclistmax + veclistblock) * sizeof( resVector ) );
      int k;
      for ( k= veclistmax; k < (veclistmax + veclistblock); k++ )
        resVectorList[k].init();
      veclistmax+= veclistblock;
      mprSTICKYPROT(ST_DENSE_MEM);
 
    }
    resVectorList[numVectors].init( mon );
    numVectors++;
    mprSTICKYPROT(ST_DENSE_NMON);
    return;
  }
  else
  {
    if ( var == (currRing->N)+1 ) return;
    poly newm = pCopy( mm );
    while ( deg >= 0 )
    {
      generateMonoms( newm, var+1, deg );
      pIncrExp( newm, var );
      pSetm( newm );
      deg--;
    }
    pDelete( & newm );
  }
 
  return;
}

getDetAt()

number resMatrixDense::getDetAt ( const number *  evpoint )

virtual

Evaluate the determinant of the matrix M at the point evpoint where the ui's are replaced by the components of evpoint.

Uses singclap_det from factory.

Reimplemented from resMatrixBase.

Definition at line 2550 of file mpr_base.cc.

{
  int k,i;
 
  // copy evaluation point into matrix
  // p0, p1, ..., pn replace u0, u1, ..., un
  for ( k= numVectors - 1; k >= 0; k-- )
  {
    if ( linPolyS == getMVector(k)->elementOfS )
    {
      for ( i= 0; i < (currRing->N); i++ )
      {
        number np=pGetCoeff(MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]));
        if (np!=NULL) nDelete(&np);
        pSetCoeff0( MATELEM(m,numVectors-k,numVectors-(getMVector(k)->numColParNr)[i]),
                   nCopy(evpoint[i]) );
      }
    }
  }
 
  mprSTICKYPROT(ST__DET);
 
  // evaluate determinant of matrix m using factory singclap_det
  poly res= singclap_det( m, currRing );
 
  // avoid errors for det==0
  number numres;
  if ( (res!=NULL)  && (!nIsZero(pGetCoeff( res ))) )
  {
    numres= nCopy( pGetCoeff( res ) );
  }
  else
  {
    numres= nInit(0);
    mprPROT("0");
  }
  pDelete( &res );
 
  mprSTICKYPROT(ST__DET);
 
  return( numres );
}

getMatrix()

ideal resMatrixDense::getMatrix ( )

virtual

Returns the matrix M in an usable presentation.

Reimplemented from resMatrixBase.

Definition at line 2469 of file mpr_base.cc.

{
  int i,j;
 
  // copy matrix
  matrix resmat= mpNew(numVectors,numVectors);
  poly p;
  for (i=1; i <= numVectors; i++)
  {
    for (j=1; j <= numVectors; j++ )
    {
      p=MATELEM(m,i,j);
      if (( p!=NULL)
      && (!nIsZero(pGetCoeff(p)))
      && (pGetCoeff(p)!=NULL)
      )
      {
        MATELEM(resmat,i,j)= pCopy( p );
      }
    }
  }
  for (i=0; i < numVectors; i++)
  {
    if ( resVectorList[i].elementOfS == linPolyS )
    {
      for (j=1; j <= (currRing->N); j++ )
      {
        if ( MATELEM(resmat,numVectors-i,
                     numVectors-resVectorList[i].numColParNr[j-1])!=NULL )
          pDelete( &MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]) );
        MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1])= pOne();
        // FIX ME
        if ( FALSE )
        {
          pSetCoeff( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), n_Param(j,currRing) );
        }
        else
        {
          pSetExp( MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]), j, 1 );
          pSetm(MATELEM(resmat,numVectors-i,numVectors-resVectorList[i].numColParNr[j-1]));
        }
      }
    }
  }
 
  // obachman: idMatrix2Module frees resmat !!
  ideal resmod= id_Matrix2Module(resmat,currRing);
  return resmod;
}

getMVector()

resVector * resMatrixDense::getMVector

(

const int

i

)

column vector of matrix, index von 0 ... numVectors-1

Definition at line 2463 of file mpr_base.cc.

{
  assume( i >= 0 && i < numVectors );
  return &resVectorList[i];
}

getSubDet()

number resMatrixDense::getSubDet ( )

virtual

Evaluates the determinant of the submatrix M'.

Since the matrix is numerically, no evaluation point is needed. Uses singclap_det from factory.

Reimplemented from resMatrixBase.

Definition at line 2593 of file mpr_base.cc.

{
  int k,i,j,l;
  resVector *vecp;
 
  // generate quadratic matrix mat of size subSize
  matrix mat= mpNew( subSize, subSize );
 
  for ( i= 1; i <= MATROWS( mat ); i++ )
  {
    for ( j= 1; j <= MATCOLS( mat ); j++ )
    {
      MATELEM(mat,i,j)= pInit();
      pSetCoeff0( MATELEM(mat,i,j), nInit(0) );
    }
  }
  j=1;
  for ( k= numVectors - 1; k >= 0; k-- )
  {
    vecp= getMVector(k);
    if ( vecp->isReduced ) continue;
    l=1;
    for ( i= numVectors - 1; i >= 0; i-- )
    {
      if ( getMVector(i)->isReduced ) continue;
      if ( vecp->getElemNum(numVectors - i - 1) && !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
      {
        pSetCoeff(MATELEM(mat, j , l ), nCopy(vecp->getElemNum(numVectors - i - 1)));
      }
      /* else
      {
           MATELEM(mat, j , l )= pOne();
           pSetCoeff(MATELEM(mat, j , l ), nInit(0) );
      }
      */
      l++;
    }
    j++;
  }
 
  poly res= singclap_det( mat, currRing );
 
  number numres;
  if ((res != NULL) && (!nIsZero(pGetCoeff( res ))) )
  {
    numres= nCopy(pGetCoeff( res ));
  }
  else
  {
    numres= nInit(0);
  }
  pDelete( &res );
  return numres;
}

getSubMatrix()

ideal resMatrixDense::getSubMatrix ( )

virtual

Returns the submatrix M' of M in an usable presentation.

Reimplemented from resMatrixBase.

Definition at line 2519 of file mpr_base.cc.

{
  int k,i,j,l;
  resVector *vecp;
 
  // generate quadratic matrix resmat of size subSize
  matrix resmat= mpNew( subSize, subSize );
 
  j=1;
  for ( k= numVectors - 1; k >= 0; k-- )
  {
    vecp= getMVector(k);
    if ( vecp->isReduced ) continue;
    l=1;
    for ( i= numVectors - 1; i >= 0; i-- )
    {
      if ( getMVector(i)->isReduced ) continue;
      if ( !nIsZero(vecp->getElemNum(numVectors - i - 1)) )
      {
        MATELEM(resmat,j,l)= pCopy( vecp->getElem(numVectors-i-1) );
      }
      l++;
    }
    j++;
  }
 
  // obachman: idMatrix2Module frees resmat !!
  ideal resmod= id_Matrix2Module(resmat,currRing);
  return resmod;
}

Field Documentation

m

matrix resMatrixDense::m

private

Definition at line 1995 of file mpr_base.cc.

numVectors

int resMatrixDense::numVectors

private

Definition at line 1992 of file mpr_base.cc.

resVectorList

resVector* resMatrixDense::resVectorList

private

Definition at line 1988 of file mpr_base.cc.

subSize

int resMatrixDense::subSize

private

Definition at line 1993 of file mpr_base.cc.

veclistblock

int resMatrixDense::veclistblock

private

Definition at line 1991 of file mpr_base.cc.

veclistmax

int resMatrixDense::veclistmax

private

Definition at line 1990 of file mpr_base.cc.

---

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

• kernel/numeric/mpr_base.cc