Re: [Bug-gsl] Re: [Help-gsl] error in Schur decomposition of a non-symmetric matrix ?

[Francesco Abbate]

2010/3/23 Patrick Alken <address@hidden>:
> On 03/23/2010 09:54 AM, Francesco Abbate wrote:
> After looking into this a bit more, I found that the code is performing
> balancing of matrices for the case of computing eigenvectors, which does not
> preserve orthogonality of the Schur vectors. I probably did it this way
> because its the LAPACK default. This issue should only affect the function
> gsl_eigen_nonsymmv_Z(). The gsl_eigen_nonsymm_Z() function should give you
> an orthogonal matrix Z, provided you haven't turned on balancing.
>
> The only way to give the user control over balancing in the eigenvector case
> would be to make another function like gsl_eigen_nonsymmv_params()...which
> may be best solution.
>
> Otherwise for now you can just use gsl_eigen_nonsymm_Z() to get the matrix
> Z, and a separate call to gsl_eigen_nonsymmv() to get the eigenvectors. Or
> if you really need something efficient right away, go to line 98 of
> eigen/nonsymmv.c and change it to:
>
> gsl_eigen_nonsymm_params(1, 0, w->nonsymm_workspace_p);

Ok, Patrick, that looks ok. Probably a small fix should be included in
the release of GSL as you suggests by adding a parameter that you can
configure also for the 'nonsymmv' function.

> I'll think a little bit more about whether it would make sense to give the
> user the option of balancing for the nonsymmv functions. Out of curiosity,
> why do you need both the Schur vectors and the eigenvectors?

I'm developing a software, GSL shell, that offer an high level
interface to GSL routines. I was implementing the routines for solving
eigensystems and for completeness I've also provided a "shur" function
that return the Schur decomposition of the matrix. If someone want to
try out, the Subverson repository at Savannah is up to date.

Here a sample session:
--------------------------------------------------------------------------------------
GSL Shell, Copyright (C) 2009 Francesco Abbate
GNU Scientific Library, Copyright (C) The GSL Team
Lua 5.1.4, Copyright (C) 1994-2008 Lua.org, PUC-Rio (complex double int32)
> dofile('examples/matrix-algebra.lua')
> A = vandermonde {-1, -2, 3, 4}
> =A
[ -1  1 -1  1 ]
[ -8  4 -2  1 ]
[ 27  9  3  1 ]
[ 64 16  4  1 ]
> e, v = eignonsymmv(A)
> =e
[         -6.41391 ]
[ 5.54555+3.08545i ]
[ 5.54555-3.08545i ]
[           2.3228 ]
> = cmul(cinverse(v), c(A), v)
[         -6.41391                0                0                0 ]
[                0 5.54555+3.08545i                0                0 ]
[                0                0 5.54555-3.08545i                0 ]
[                0                0                0           2.3228 ]
> T, Z = schur(A)
> = mul(Z, T, inverse(Z))
[ -1  1 -1  1 ]
[ -8  4 -2  1 ]
[ 27  9  3  1 ]
[ 64 16  4  1 ]
----------------------------------------------------------

Francesco