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

[Patrick Alken]

On 03/23/2010 09:54 AM, Francesco Abbate wrote:
> > Perhaps a more clear explanation of the Schur form is this: if you take a
> > given 2x2 block on the diagonal of the matrix, call it
> >
> > [ T11 T12 ]
> > [ T21 T22 ]
> >
> > then either T21 = 0, or T21 is nonzero.
> >
> > If T21 is zero, then T11 and T22 are both real eigenvalues of the original
> > matrix.
> >
> > If T21 is nonzero, then T11 = T22, and T11 +/- sqrt(|T21*T12|) are complex
> > conjugate eigenvalues.
> >
> > So, the subdiagonal element (T21) is guaranteed to be 0 in the case of real
> > eigenvalues, and nonzero in the case of complex eigenvalues. I will try to
> > update the documentation to make this a little more clear.
> >
> > The rest of A (below the subdiagonal) is not guaranteed to be 0, as far as I
> > remember.
> >      
> Thank you very much Patrick, now it is completely clear. Well,
> actually this is more or less what I have guessed but I was not sure
> especially about the subdiagonal and the lower part of the matrix.
> Actually I was thinking that for complex eigenvalues the submatrix was
> something like:
>
> [ u  v]
> [ -v u]
>
> à la Cauchy-Riemann but it was a wrong guess... :-)
>
> Otherwise only the problem of the orthogonality of Z remains open...
> maybe it is due to the balancing of the matrix made (maybe) internally
> by the procedure...
>
> In any case, thank you very much for your help!
>
> Francesco
>
>
> _______________________________________________
> Bug-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/bug-gsl
>    

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);

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?