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