# source:git/kernel/eigenval.cc@bc669c

fieker-DuValspielwiese
Last change on this file since bc669c was 341696, checked in by Hans Schönemann <hannes@…>, 15 years ago
• Property mode set to `100644`
File size: 2.1 KB
Line
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* \$Id\$ */
5/*
6* ABSTRACT: eigenvalues of constant square matrices
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_EIGENVAL
12
13#include "febase.h"
14#include "structs.h"
15//#include "ipid.h"
16#include "intvec.h"
17#include "numbers.h"
18#include "polys.h"
19#include "ideals.h"
20//#include "lists.h"
21#include "matpol.h"
22#include "clapsing.h"
23#include "eigenval.h"
24
25
26matrix evSwap(matrix M,int i,int j)
27{
28  if(i==j)
29    return(M);
30
31  for(int k=1;k<=MATROWS(M);k++)
32  {
33    poly p=MATELEM(M,i,k);
34    MATELEM(M,i,k)=MATELEM(M,j,k);
35    MATELEM(M,j,k)=p;
36  }
37
38  for(int k=1;k<=MATCOLS(M);k++)
39  {
40    poly p=MATELEM(M,k,i);
41    MATELEM(M,k,i)=MATELEM(M,k,j);
42    MATELEM(M,k,j)=p;
43  }
44
45  return(M);
46}
47
48matrix evRowElim(matrix M,int i,int j,int k)
49{
50  if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
51    return(M);
52
53  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k))));
54  pNormalize(p);
55
56  for(int l=1;l<=MATCOLS(M);l++)
57  {
58    MATELEM(M,i,l)=pSub(MATELEM(M,i,l),ppMult_qq(p,MATELEM(M,j,l)));
59    pNormalize(MATELEM(M,i,l));
60  }
61  for(int l=1;l<=MATROWS(M);l++)
62  {
64    pNormalize(MATELEM(M,l,j));
65  }
66
67  pDelete(&p);
68
69  return(M);
70}
71
72matrix evColElim(matrix M,int i,int j,int k)
73{
74  if(MATELEM(M,k,i)==0||MATELEM(M,k,j)==0)
75    return(M);
76
77  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j))));
78  pNormalize(p);
79
80  for(int l=1;l<=MATROWS(M);l++)
81  {
82    MATELEM(M,l,i)=pSub(MATELEM(M,l,i),ppMult_qq(p,MATELEM(M,l,j)));
83    pNormalize(MATELEM(M,l,i));
84  }
85  for(int l=1;l<=MATCOLS(M);l++)
86  {
88    pNormalize(MATELEM(M,j,l));
89  }
90
91  pDelete(&p);
92
93  return(M);
94}
95
96matrix evHessenberg(matrix M)
97{
98  int n=MATROWS(M);
99  if(n!=MATCOLS(M))
100    return(M);
101
102  for(int k=1,j=2;k<n-1;k++,j=k+1)
103  {
104    while(j<=n&&MATELEM(M,j,k)==0)
105      j++;
106
107    if(j<=n)
108    {
109      M=evSwap(M,j,k+1);
110
111      for(int i=j+1;i<=n;i++)
112        M=evRowElim(M,i,k+1,k);
113    }
114  }
115
116  return(M);
117}
118
119#endif /* HAVE_EIGENVAL */
Note: See TracBrowser for help on using the repository browser.