# source:git/Singular/linearAlgebra_ip.cc@a2174a

spielwiese
Last change on this file since a2174a was a2174a, checked in by Hans Schoenemann <hannes@…>, 12 years ago
fix linearAlgebra_ip.cc: includes
• Property mode set to `100644`
File size: 3.7 KB
Line
1#include <kernel/mod2.h>
2#include <Singular/lists.h>
3#include <kernel/liniearAlgebra.h>
4
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 * _[i] is the i-th mutually distinct eigenvalue that was found,
15 * _[i] is the (int) multiplicity of _[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 **/
29lists 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
38lists 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 = 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.rtyp = INT_CMD;
58    result->m.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.rtyp = LIST_CMD;
96    result->m.data = (char*)eigenvalues;
97    result->m.rtyp = LIST_CMD;
98    result->m.data = (char*)multiplicities;
99  }
100  return result;
101}
102
Note: See TracBrowser for help on using the repository browser.