1 | /***************************************** |
---|
2 | * Computer Algebra System SINGULAR * |
---|
3 | *****************************************/ |
---|
4 | /* |
---|
5 | * ABSTRACT: eigenvalues of constant square matrices |
---|
6 | */ |
---|
7 | |
---|
8 | #include "config.h" |
---|
9 | #include <kernel/mod2.h> |
---|
10 | |
---|
11 | #ifdef HAVE_EIGENVAL |
---|
12 | |
---|
13 | #include <kernel/febase.h> |
---|
14 | #include <kernel/structs.h> |
---|
15 | //#include "ipid.h" |
---|
16 | #include <misc/intvec.h> |
---|
17 | #include <coeffs/numbers.h> |
---|
18 | #include <kernel/polys.h> |
---|
19 | #include <kernel/ideals.h> |
---|
20 | //#include "lists.h" |
---|
21 | #include <polys/matpol.h> |
---|
22 | #include <polys/clapsing.h> |
---|
23 | #include <kernel/eigenval.h> |
---|
24 | |
---|
25 | |
---|
26 | matrix 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 | |
---|
48 | matrix 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 | { |
---|
63 | MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),ppMult_qq(p,MATELEM(M,l,i))); |
---|
64 | pNormalize(MATELEM(M,l,j)); |
---|
65 | } |
---|
66 | |
---|
67 | pDelete(&p); |
---|
68 | |
---|
69 | return(M); |
---|
70 | } |
---|
71 | |
---|
72 | matrix 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 | { |
---|
87 | MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),ppMult_qq(p,MATELEM(M,i,l))); |
---|
88 | pNormalize(MATELEM(M,j,l)); |
---|
89 | } |
---|
90 | |
---|
91 | pDelete(&p); |
---|
92 | |
---|
93 | return(M); |
---|
94 | } |
---|
95 | |
---|
96 | matrix 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 */ |
---|