My Project
Loading...
Searching...
No Matches
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
41lists 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:132
int rtyp
Definition: subexpr.h:91
void * data
Definition: subexpr.h:88
Definition: lists.h:24
sleftv * m
Definition: lists.h:46
INLINE_THIS void Init(int l=0)
return result
Definition: facAbsBiFact.cc:75
STATIC_VAR int mults
Definition: fast_mult.cc:13
@ NUMBER_CMD
Definition: grammar.cc:288
#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:57
#define MATROWS(i)
Definition: matpol.h:26
slists * lists
Definition: mpr_numeric.h:146
#define nDelete(n)
Definition: numbers.h:16
#define nCopy(n)
Definition: numbers.h:15
#define omAlloc(size)
Definition: omAllocDecl.h:210
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
VAR 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:27
#define A
Definition: sirandom.c:24
@ LIST_CMD
Definition: tok.h:118
@ INT_CMD
Definition: tok.h:96