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