# source:git/Singular/linearAlgebra.cut@6ed8c4

spielwiese
Last change on this file since 6ed8c4 was 6ed8c4, checked in by mlee <martinlee84@…>, 12 years ago
added const ring to some functions in linearAlgebra.* (TODO complete this) switched on compatibility layers in numbers.h and polys.h added nGreater macro to numbers.h commented out lists in ring.h moved lists dependend functions from linearAlgebra.* to linearAlgebra.cut
• Property mode set to `100644`
File size: 3.6 KB
Line
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 * _[i] is the i-th mutually distinct eigenvalue that was found,
11 * _[i] is the (int) multiplicity of _[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 **/
25lists 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
34lists 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 = 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.rtyp = INT_CMD;
54    result->m.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.rtyp = LIST_CMD;
92    result->m.data = (char*)eigenvalues;
93    result->m.rtyp = LIST_CMD;
94    result->m.data = (char*)multiplicities;
95  }
96  return result;
97}
98
Note: See TracBrowser for help on using the repository browser.