My Project
Loading...
Searching...
No Matches
Functions
linearAlgebra_ip.h File Reference
#include "Singular/lists.h"
#include "kernel/linear_algebra/linearAlgebra.h"

Go to the source code of this file.

Functions

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. More...
 

Function Documentation

◆ qrDoubleShift()

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.

The method assumes that the current ground field is the complex numbers. Computations are based on the QR double shift algorithm involving Hessenberg form and householder transformations. If the algorithm works, then it returns a list with two entries which are again lists of the same size: _[1][i] is the i-th mutually distinct eigenvalue that was found, _[2][i] is the (int) multiplicity of _[1][i]. If the algorithm does not work (due to an ill-posed matrix), a list with the single entry (int)0 is returned. 'tol1' is used for detection of deflation in the actual qr double shift algorithm. 'tol2' is used for ending Heron's iteration whenever square roots are being computed. 'tol3' is used to distinguish between distinct eigenvalues: When the Euclidean distance between two computed eigenvalues is less then tol3, then they will be regarded equal, resulting in a higher multiplicity of the corresponding eigenvalue.

Returns
a list with one entry (int)0, or two entries which are again lists
Parameters
[in]Athe quadratic matrix
[in]tol1tolerance for deflation
[in]tol2tolerance for square roots
[in]tol3tolerance for distinguishing eigenvalues

Definition at line 41 of file linearAlgebra_ip.cc.

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}
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].
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
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
@ LIST_CMD
Definition: tok.h:118
@ INT_CMD
Definition: tok.h:96