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