eigenval_ip.cc File Reference

#include "kernel/mod2.h"
#include "Singular/tok.h"
#include "Singular/ipid.h"
#include "misc/intvec.h"
#include "coeffs/numbers.h"
#include "kernel/polys.h"
#include "kernel/ideals.h"
#include "Singular/lists.h"
#include "polys/matpol.h"
#include "polys/clapsing.h"
#include "kernel/linear_algebra/eigenval.h"
#include "Singular/ipshell.h"
#include "Singular/eigenval_ip.h"

Go to the source code of this file.

Functions
BOOLEAN	evSwap (leftv res, leftv h)
BOOLEAN	evRowElim (leftv res, leftv h)
BOOLEAN	evHessenberg (leftv res, leftv h)
lists	evEigenvals (matrix M)
BOOLEAN	evEigenvals (leftv res, leftv h)

Function Documentation

evEigenvals() [1/2]

BOOLEAN evEigenvals	(	leftv	res,
		leftv	h
	)

Definition at line 269 of file eigenval_ip.cc.

{
  if(currRing)
  {
    if(h&&h->Typ()==MATRIX_CMD)
    {
      matrix M=(matrix)h->CopyD();
      res->rtyp=LIST_CMD;
      res->data=(void *)evEigenvals(M);
      return FALSE;
    }
    WerrorS("<matrix> expected");
    return TRUE;
  }
  WerrorS("no ring active");
  return TRUE;
}

evEigenvals() [2/2]

lists evEigenvals

(

matrix

M

)

Definition at line 120 of file eigenval_ip.cc.

{
  lists l=(lists)omAllocBin(slists_bin);
  if(MATROWS(M)!=MATCOLS(M))
  {
    l->Init(0);
    return(l);
  }
 
  M=evHessenberg(M);
 
  int n=MATCOLS(M);
  ideal e=idInit(n,1);
  intvec *m=new intvec(n);
 
  poly t=pOne();
  pSetExp(t,1,1);
  pSetm(t);
 
  for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
  {
    while(j<=n&&MATELEM(M,j,j-1)!=NULL)
      j++;
    if(j==j0+1)
    {
      e->m[k]=pHead(MATELEM(M,j0,j0));
      (*m)[k]=1;
      k++;
    }
    else
    {
      int n0=j-j0;
      matrix M0=mpNew(n0,n0);
 
      j0--;
      for(int i=1;i<=n0;i++)
        for(int j=1;j<=n0;j++)
          MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
      for(int i=1;i<=n0;i++)
        MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
 
      intvec *m0;
      ideal e0=singclap_factorize(mp_DetBareiss(M0,currRing),&m0,2, currRing);
      if (e0==NULL)
      {
        l->Init(0);
        return(l);
      }
 
      for(int i=0;i<IDELEMS(e0);i++)
      {
        if(pNext(e0->m[i])==NULL)
        {
          (*m)[k]=(*m0)[i];
          k++;
        }
        else
        if(pGetExp(e0->m[i],1)<2&&pGetExp(pNext(e0->m[i]),1)<2&&
           pNext(pNext(e0->m[i]))==NULL)
        {
          number e1=nCopy(pGetCoeff(e0->m[i]));
          e1=nInpNeg(e1);
          if(pGetExp(pNext(e0->m[i]),1)==0)
            e->m[k]=pNSet(nDiv(pGetCoeff(pNext(e0->m[i])),e1));
          else
            e->m[k]=pNSet(nDiv(e1,pGetCoeff(pNext(e0->m[i]))));
          nDelete(&e1);
          pNormalize(e->m[k]);
          (*m)[k]=(*m0)[i];
          k++;
        }
        else
        {
          e->m[k]=e0->m[i];
          pNormalize(e->m[k]);
          e0->m[i]=NULL;
          (*m)[k]=(*m0)[i];
          k++;
        }
      }
 
      delete(m0);
      idDelete(&e0);
    }
  }
 
  pDelete(&t);
  idDelete((ideal *)&M);
 
  for(int i0=0;i0<n-1;i0++)
  {
    for(int i1=i0+1;i1<n;i1++)
    {
      if(pEqualPolys(e->m[i0],e->m[i1]))
      {
        (*m)[i0]+=(*m)[i1];
        (*m)[i1]=0;
      }
      else
      {
        if((e->m[i0]==NULL&&!nGreaterZero(pGetCoeff(e->m[i1])))||
           (e->m[i1]==NULL&&
            (nGreaterZero(pGetCoeff(e->m[i0]))||pNext(e->m[i0])!=NULL))||
           e->m[i0]!=NULL&&e->m[i1]!=NULL&&
            ((pNext(e->m[i0])!=NULL&&pNext(e->m[i1])==NULL)||
             (pNext(e->m[i0])==NULL&&pNext(e->m[i1])==NULL&&
             nGreater(pGetCoeff(e->m[i0]),pGetCoeff(e->m[i1])))))
        {
          poly e1=e->m[i0];
          e->m[i0]=e->m[i1];
          e->m[i1]=e1;
          int m1=(*m)[i0];
          (*m)[i0]=(*m)[i1];
          (*m)[i1]=m1;
        }
      }
    }
  }
 
  int n0=0;
  for(int i=0;i<n;i++)
    if((*m)[i]>0)
      n0++;
 
  ideal e0=idInit(n0,1);
  intvec *m0=new intvec(n0);
 
  for(int i=0,i0=0;i<n;i++)
    if((*m)[i]>0)
    {
      e0->m[i0]=e->m[i];
      e->m[i]=NULL;
      (*m0)[i0]=(*m)[i];
      i0++;
    }
 
  idDelete(&e);
  delete(m);
 
  l->Init(2);
  l->m[0].rtyp=IDEAL_CMD;
  l->m[0].data=e0;
  l->m[1].rtyp=INTVEC_CMD;
  l->m[1].data=m0;
 
  return(l);
}

evHessenberg()

BOOLEAN evHessenberg	(	leftv	res,
		leftv	h
	)

Definition at line 101 of file eigenval_ip.cc.

{
  if(currRing)
  {
    if(h&&h->Typ()==MATRIX_CMD)
    {
      matrix M=(matrix)h->Data();
      res->rtyp=MATRIX_CMD;
      res->data=(void *)evHessenberg(mp_Copy(M, currRing));
      return FALSE;
    }
    WerrorS("<matrix> expected");
    return TRUE;
  }
  WerrorS("no ring active");
  return TRUE;
}

evRowElim()

BOOLEAN evRowElim	(	leftv	res,
		leftv	h
	)

Definition at line 51 of file eigenval_ip.cc.

{
  if(currRing)
  {
    const short t[]={4,MATRIX_CMD,INT_CMD,INT_CMD,INT_CMD};
    if (iiCheckTypes(h,t,1))
    {
      matrix M=(matrix)h->CopyD();
      h=h->next;
      int i=(int)(long)h->Data();
      h=h->next;
      int j=(int)(long)h->Data();
      h=h->next;
      int k=(int)(long)h->Data();
      res->rtyp=MATRIX_CMD;
      res->data=(void *)evRowElim(M,i,j,k);
      return FALSE;
    }
    return TRUE;
  }
  WerrorS("no ring active");
  return TRUE;
}

evSwap()

BOOLEAN evSwap	(	leftv	res,
		leftv	h
	)

Definition at line 29 of file eigenval_ip.cc.

{
  if(currRing)
  {
    const short t[]={3,MATRIX_CMD,INT_CMD,INT_CMD};
    if (iiCheckTypes(h,t,1))
    {
      matrix M=(matrix)h->Data();
      h=h->next;
      int i=(int)(long)h->Data();
      h=h->next;
      int j=(int)(long)h->Data();
      res->rtyp=MATRIX_CMD;
      res->data=(void *)evSwap(mp_Copy(M, currRing),i,j);
      return FALSE;
    }
    return TRUE;
  }
  WerrorS("no ring active");
  return TRUE;
}