/*****************************************
*  Computer Algebra System SINGULAR      *
*****************************************/
/* $Id: eigenval.cc,v 1.3 2004-03-04 10:08:24 Singular Exp $ */
/*
* ABSTRACT: eigenvalues of constant square matrices
*/

#include "mod2.h"

#ifdef HAVE_EIGENVAL

#include "febase.h"
#include "structs.h"
//#include "ipid.h"
#include "intvec.h"
#include "numbers.h"
#include "polys.h"
#include "ideals.h"
//#include "lists.h"
#include "matpol.h"
#include "clapsing.h"
#include "eigenval.h"


matrix evSwap(matrix M,int i,int j)
{
  if(i==j)
    return(M);

  for(int k=1;k<=MATROWS(M);k++)
  {
    poly p=MATELEM(M,i,k);
    MATELEM(M,i,k)=MATELEM(M,j,k);
    MATELEM(M,j,k)=p;
  }

  for(int k=1;k<=MATCOLS(M);k++)
  {
    poly p=MATELEM(M,k,i);
    MATELEM(M,k,i)=MATELEM(M,k,j);
    MATELEM(M,k,j)=p;
  }

  return(M);
}

matrix evRowElim(matrix M,int i,int j,int k)
{
  if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
    return(M);

  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k))));
  pNormalize(p);

  for(int l=1;l<=MATCOLS(M);l++)
  {
    MATELEM(M,i,l)=pSub(MATELEM(M,i,l),ppMult_qq(p,MATELEM(M,j,l)));
    pNormalize(MATELEM(M,i,l));
  }
  for(int l=1;l<=MATROWS(M);l++)
  {
    MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),ppMult_qq(p,MATELEM(M,l,i)));
    pNormalize(MATELEM(M,l,j));
  }

  pDelete(&p);

  return(M);
}

matrix evColElim(matrix M,int i,int j,int k)
{
  if(MATELEM(M,k,i)==0||MATELEM(M,k,j)==0)
    return(M);

  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j))));
  pNormalize(p);

  for(int l=1;l<=MATROWS(M);l++)
  {
    MATELEM(M,l,i)=pSub(MATELEM(M,l,i),ppMult_qq(p,MATELEM(M,l,j)));
    pNormalize(MATELEM(M,l,i));
  }
  for(int l=1;l<=MATCOLS(M);l++)
  {
    MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),ppMult_qq(p,MATELEM(M,i,l)));
    pNormalize(MATELEM(M,j,l));
  }

  pDelete(&p);

  return(M);
}

matrix evHessenberg(matrix M)
{
  int n=MATROWS(M);
  if(n!=MATCOLS(M))
    return(M);

  for(int k=1,j=2;k<n-1;k++,j=k+1)
  {
    while(j<=n&&MATELEM(M,j,k)==0)
      j++;

    if(j<=n)
    {
      M=evSwap(M,j,k+1);

      for(int i=j+1;i<=n;i++)
        M=evRowElim(M,i,k+1,k);
    }
  }

  return(M);
}

#endif /* HAVE_EIGENVAL */