source: git/Singular/linearAlgebra_ip.cc @ b2aa08

spielwiese
Last change on this file since b2aa08 was 762407, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
config.h is for sources files only FIX: config.h should only be used by source (not from inside kernel/mod2.h!) NOTE: each source file should better include mod2.h right after config.h, while headers should better not include mod2.h.
  • Property mode set to 100644
File size: 3.7 KB
Line 
1#include "config.h"
2#include <kernel/mod2.h>
3#include <Singular/lists.h>
4#include <kernel/linearAlgebra.h>
5
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 **/
30lists 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
39lists 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
Note: See TracBrowser for help on using the repository browser.