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