1 /*****************************************
2 * Computer Algebra System SINGULAR *
3 *****************************************/
4 /*
5 * ABSTRACT: eigenvalues of constant square matrices
6 */
7
8
9
10
11 #include "kernel/mod2.h"
12
13 #ifdef HAVE_EIGENVAL
14
15 #include "kernel/structs.h"
16 #include "misc/intvec.h"
17 #include "coeffs/numbers.h"
18 #include "kernel/polys.h"
19 #include "kernel/ideals.h"
20 #include "polys/matpol.h"
21 #include "polys/clapsing.h"
22 #include "kernel/linear_algebra/eigenval.h"
23
24
evSwap(matrix M,int i,int j)25 matrix evSwap(matrix M,int i,int j)
26 {
27 if(i==j)
28 return(M);
29
30 for(int k=1;k<=MATROWS(M);k++)
31 {
32 poly p=MATELEM(M,i,k);
33 MATELEM(M,i,k)=MATELEM(M,j,k);
34 MATELEM(M,j,k)=p;
35 }
36
37 for(int k=1;k<=MATCOLS(M);k++)
38 {
39 poly p=MATELEM(M,k,i);
40 MATELEM(M,k,i)=MATELEM(M,k,j);
41 MATELEM(M,k,j)=p;
42 }
43
44 return(M);
45 }
46
evRowElim(matrix M,int i,int j,int k)47 matrix evRowElim(matrix M,int i,int j,int k)
48 {
49 if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
50 return(M);
51 poly p1=pp_Jet(MATELEM(M,i,k),0,currRing);
52 poly p2=pp_Jet(MATELEM(M,j,k),0,currRing);
53 if ((p1==NULL)||(p2==NULL)) return (M);
54
55 poly p=pNSet(nDiv(pGetCoeff(p1),pGetCoeff(p2)));
56 pNormalize(p);
57
58 for(int l=1;l<=MATCOLS(M);l++)
59 {
60 MATELEM(M,i,l)=pSub(MATELEM(M,i,l),ppMult_qq(p,MATELEM(M,j,l)));
61 pNormalize(MATELEM(M,i,l));
62 }
63 for(int l=1;l<=MATROWS(M);l++)
64 {
65 MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),ppMult_qq(p,MATELEM(M,l,i)));
66 pNormalize(MATELEM(M,l,j));
67 }
68
69 pDelete(&p);
70 pDelete(&p1);
71 pDelete(&p2);
72
73 return(M);
74 }
75
evColElim(matrix M,int i,int j,int k)76 matrix evColElim(matrix M,int i,int j,int k)
77 {
78 if(MATELEM(M,k,i)==0||MATELEM(M,k,j)==0)
79 return(M);
80
81 poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j))));
82 pNormalize(p);
83
84 for(int l=1;l<=MATROWS(M);l++)
85 {
86 MATELEM(M,l,i)=pSub(MATELEM(M,l,i),ppMult_qq(p,MATELEM(M,l,j)));
87 pNormalize(MATELEM(M,l,i));
88 }
89 for(int l=1;l<=MATCOLS(M);l++)
90 {
91 MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),ppMult_qq(p,MATELEM(M,i,l)));
92 pNormalize(MATELEM(M,j,l));
93 }
94
95 pDelete(&p);
96
97 return(M);
98 }
99
evHessenberg(matrix M)100 matrix evHessenberg(matrix M)
101 {
102 int n=MATROWS(M);
103 if(n!=MATCOLS(M))
104 return(M);
105
106 for(int k=1,j=2;k<n-1;k++,j=k+1)
107 {
108 while((j<=n)
109 &&((MATELEM(M,j,k)==NULL)
110 || (p_Totaldegree(MATELEM(M,j,k),currRing)!=0)))
111 j++;
112
113 if(j<=n)
114 {
115 M=evSwap(M,j,k+1);
116
117 for(int i=j+1;i<=n;i++)
118 M=evRowElim(M,i,k+1,k);
119 }
120 }
121
122 return(M);
123 }
124
125 #endif /* HAVE_EIGENVAL */
126