Loading...
Searching...
No Matches
linearAlgebra_ip.cc
Go to the documentation of this file.
#include "kernel/mod2.h"
#include "Singular/lists.h"
#include "kernel/linear_algebra/linearAlgebra.h"
/**
* 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.
*
* @return a list with one entry (int)0, or two entries which are again lists
**/
lists qrDoubleShift(
const matrix A, /**< [in] the quadratic matrix */
const number tol1, /**< [in] tolerance for deflation */
const number tol2, /**< [in] tolerance for square roots */
const number tol3, /**< [in] tolerance for distinguishing
eigenvalues */
const ring r= currRing
);
lists qrDoubleShift(const matrix A, const number tol1, const number tol2,
const number tol3, const ring R)
{
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
linearAlgebra.h
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.
Definition linearAlgebra_ip.cc:41
lists.h
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
mod2.h
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
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
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