My Project  debian-1:4.1.1-p2+ds-4build3
Functions
linearAlgebra_ip.cc File Reference
#include "kernel/mod2.h"
#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]A[in] the quadratic matrix
[in]tol1[in] tolerance for deflation
[in]tol2[in] tolerance for square roots
[in]tol3[in] tolerance 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 }
ip_smatrix
Definition: matpol.h:15
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
NUMBER_CMD
@ NUMBER_CMD
Definition: grammar.cc:287
result
return result
Definition: facAbsBiFact.cc:76
LIST_CMD
@ LIST_CMD
Definition: tok.h:118
qrDS
bool qrDS(const int, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
Definition: linearAlgebra.cc:1090
similar
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].
Definition: linearAlgebra.cc:1188
i
int i
Definition: cfEzgcd.cc:125
INT_CMD
@ INT_CMD
Definition: tok.h:96
sleftv::data
void * data
Definition: subexpr.h:88
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:210
slists::m
sleftv * m
Definition: lists.h:45
slists
Definition: lists.h:23
sleftv::rtyp
int rtyp
Definition: subexpr.h:91
lists
slists * lists
Definition: mpr_numeric.h:146
mults
static int mults
Definition: fast_mult.cc:13
nDelete
#define nDelete(n)
Definition: numbers.h:17
R
#define R
Definition: sirandom.c:26
slists::Init
INLINE_THIS void Init(int l=0)
mp_Copy
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:63
MATROWS
#define MATROWS(i)
Definition: matpol.h:28
index
static int index(p_Length length, p_Ord ord)
Definition: p_Procs_Impl.h:592
A
#define A
Definition: sirandom.c:23
nCopy
#define nCopy(n)
Definition: numbers.h:16