1 | /** |
---|
2 | * Computes all eigenvalues of a given real quadratic matrix with |
---|
3 | * multiplicites. |
---|
4 | * |
---|
5 | * The method assumes that the current ground field is the complex numbers. |
---|
6 | * Computations are based on the QR double shift algorithm involving |
---|
7 | * Hessenberg form and householder transformations. |
---|
8 | * If the algorithm works, then it returns a list with two entries which |
---|
9 | * are again lists of the same size: |
---|
10 | * _[1][i] is the i-th mutually distinct eigenvalue that was found, |
---|
11 | * _[2][i] is the (int) multiplicity of _[1][i]. |
---|
12 | * If the algorithm does not work (due to an ill-posed matrix), a list with |
---|
13 | * the single entry (int)0 is returned. |
---|
14 | * 'tol1' is used for detection of deflation in the actual qr double shift |
---|
15 | * algorithm. |
---|
16 | * 'tol2' is used for ending Heron's iteration whenever square roots |
---|
17 | * are being computed. |
---|
18 | * 'tol3' is used to distinguish between distinct eigenvalues: When |
---|
19 | * the Euclidean distance between two computed eigenvalues is less then |
---|
20 | * tol3, then they will be regarded equal, resulting in a higher |
---|
21 | * multiplicity of the corresponding eigenvalue. |
---|
22 | * |
---|
23 | * @return a list with one entry (int)0, or two entries which are again lists |
---|
24 | **/ |
---|
25 | lists qrDoubleShift( |
---|
26 | const matrix A, /**< [in] the quadratic matrix */ |
---|
27 | const number tol1, /**< [in] tolerance for deflation */ |
---|
28 | const number tol2, /**< [in] tolerance for square roots */ |
---|
29 | const number tol3, /**< [in] tolerance for distinguishing |
---|
30 | eigenvalues */ |
---|
31 | const ring r= currRing |
---|
32 | ); |
---|
33 | |
---|
34 | lists qrDoubleShift(const matrix A, const number tol1, const number tol2, |
---|
35 | const number tol3, const ring R) |
---|
36 | { |
---|
37 | int n = MATROWS(A); |
---|
38 | matrix* queue = new matrix[n]; |
---|
39 | queue[0] = mp_Copy(A,R); int queueL = 1; |
---|
40 | number* eigenVs = new number[n]; int eigenL = 0; |
---|
41 | /* here comes the main call: */ |
---|
42 | bool worked = qrDS(n, queue, queueL, eigenVs, eigenL, tol1, tol2,R); |
---|
43 | lists result = (lists)omAlloc(sizeof(slists)); |
---|
44 | if (!worked) |
---|
45 | { |
---|
46 | for (int i = 0; i < eigenL; i++) |
---|
47 | nDelete(&eigenVs[i]); |
---|
48 | delete [] eigenVs; |
---|
49 | for (int i = 0; i < queueL; i++) |
---|
50 | idDelete((ideal*)&queue[i]); |
---|
51 | delete [] queue; |
---|
52 | result->Init(1); |
---|
53 | result->m[0].rtyp = INT_CMD; |
---|
54 | result->m[0].data = (void*)0; /* a list with a single entry |
---|
55 | which is the int zero */ |
---|
56 | } |
---|
57 | else |
---|
58 | { |
---|
59 | /* now eigenVs[0..eigenL-1] contain all eigenvalues; among them, there |
---|
60 | may be equal entries */ |
---|
61 | number* distinctEVs = new number[n]; int distinctC = 0; |
---|
62 | int* mults = new int[n]; |
---|
63 | for (int i = 0; i < eigenL; i++) |
---|
64 | { |
---|
65 | int index = similar(distinctEVs, distinctC, eigenVs[i], tol3); |
---|
66 | if (index == -1) /* a new eigenvalue */ |
---|
67 | { |
---|
68 | distinctEVs[distinctC] = nCopy(eigenVs[i]); |
---|
69 | mults[distinctC++] = 1; |
---|
70 | } |
---|
71 | else mults[index]++; |
---|
72 | nDelete(&eigenVs[i]); |
---|
73 | } |
---|
74 | delete [] eigenVs; |
---|
75 | |
---|
76 | lists eigenvalues = (lists)omAlloc(sizeof(slists)); |
---|
77 | eigenvalues->Init(distinctC); |
---|
78 | lists multiplicities = (lists)omAlloc(sizeof(slists)); |
---|
79 | multiplicities->Init(distinctC); |
---|
80 | for (int i = 0; i < distinctC; i++) |
---|
81 | { |
---|
82 | eigenvalues->m[i].rtyp = NUMBER_CMD; |
---|
83 | eigenvalues->m[i].data = (void*)nCopy(distinctEVs[i]); |
---|
84 | multiplicities->m[i].rtyp = INT_CMD; |
---|
85 | multiplicities->m[i].data = (void*)mults[i]; |
---|
86 | nDelete(&distinctEVs[i]); |
---|
87 | } |
---|
88 | delete [] distinctEVs; delete [] mults; |
---|
89 | |
---|
90 | result->Init(2); |
---|
91 | result->m[0].rtyp = LIST_CMD; |
---|
92 | result->m[0].data = (char*)eigenvalues; |
---|
93 | result->m[1].rtyp = LIST_CMD; |
---|
94 | result->m[1].data = (char*)multiplicities; |
---|
95 | } |
---|
96 | return result; |
---|
97 | } |
---|
98 | |
---|