Loading...
Searching...
No Matches
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.

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

Computes all eigenvalues of a given real quadratic matrix with multiplicities.

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 | the quadratic matrix
  [in] | tol1 | tolerance for deflation
  [in] | tol2 | tolerance for square roots
  [in] | tol3 | tolerance for distinguishing eigenvalues

Definition at line 41 of file linearAlgebra_ip.cc.

{
int n = MATROWS(A);
matrix* queue = new matrix[n];
queue[0] = mp_Copy(A,R); int queueL = 1;
number* eigenVs = new number[n]; int eigenL = 0;
/* here comes the main call: */
bool worked = qrDS(n, queue, queueL, eigenVs, eigenL, tol1, tol2,R);
lists result = (lists)omAlloc(sizeof(slists));
if (!worked)
{
for (int i = 0; i < eigenL; i++)
nDelete(&eigenVs[i]);
delete [] eigenVs;
for (int i = 0; i < queueL; i++)
idDelete((ideal*)&queue[i]);
delete [] queue;
result->Init(1);
result->m[0].rtyp = INT_CMD;
result->m[0].data = (void*)0; /* a list with a single entry
which is the int zero */
}
else
{
/* now eigenVs[0..eigenL-1] contain all eigenvalues; among them, there
may be equal entries */
number* distinctEVs = new number[n]; int distinctC = 0;
int* mults = new int[n];
for (int i = 0; i < eigenL; i++)
{
int index = similar(distinctEVs, distinctC, eigenVs[i], tol3);
if (index == -1) /* a new eigenvalue */
{
distinctEVs[distinctC] = nCopy(eigenVs[i]);
mults[distinctC++] = 1;
}
else mults[index]++;
nDelete(&eigenVs[i]);
}
delete [] eigenVs;
lists eigenvalues = (lists)omAlloc(sizeof(slists));
eigenvalues->Init(distinctC);
lists multiplicities = (lists)omAlloc(sizeof(slists));
multiplicities->Init(distinctC);
for (int i = 0; i < distinctC; i++)
{
eigenvalues->m[i].rtyp = NUMBER_CMD;
eigenvalues->m[i].data = (void*)nCopy(distinctEVs[i]);
multiplicities->m[i].rtyp = INT_CMD;
multiplicities->m[i].data = (void*)(long)mults[i];
nDelete(&distinctEVs[i]);
}
delete [] distinctEVs; delete [] mults;
result->Init(2);
result->m[0].rtyp = LIST_CMD;
result->m[0].data = (char*)eigenvalues;
result->m[1].rtyp = LIST_CMD;
result->m[1].data = (char*)multiplicities;
}
return result;
}
i
int i
Definition cfEzgcd.cc:132
sleftv::rtyp
int rtyp
Definition subexpr.h:91
sleftv::data
void * data
Definition subexpr.h:88
slists
Definition lists.h:24
slists::m
sleftv * m
Definition lists.h:46
slists::Init
INLINE_THIS void Init(int l=0)
result
return result
Definition facAbsBiFact.cc:76
mults
STATIC_VAR int mults
Definition fast_mult.cc:13
NUMBER_CMD
@ NUMBER_CMD
Definition grammar.cc:289
idDelete
#define idDelete(H)
delete an ideal
Definition ideals.h:29
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:1099
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:1197
mp_Copy
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition matpol.cc:57
matrix
ip_smatrix * matrix
Definition matpol.h:43
MATROWS
#define MATROWS(i)
Definition matpol.h:26
lists
slists * lists
Definition mpr_numeric.h:146
nDelete
#define nDelete(n)
Definition numbers.h:16
nCopy
#define nCopy(n)
Definition numbers.h:15
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
index
static int index(p_Length length, p_Ord ord)
Definition p_Procs_Impl.h:592
R
#define R
Definition sirandom.c:27
A
#define A
Definition sirandom.c:24
LIST_CMD
@ LIST_CMD
Definition tok.h:118
INT_CMD
@ INT_CMD
Definition tok.h:96