Re: [ESA-SoC] Generalized eigenvalue

[eem2314]

Rafael Gonzalez wrote
> Hello
> 
> I tried to address this problem a few years ago but I couldn't continue in
> the GSoC program. Later, someone else also tried to take the problem but I
> stop hearing from him. I'm copying you the mail I wrote to him where I
> explain what is what I think is the issue.
> 
> -----
> 
> I'm answering between lines:
> 
>>> I can't see any FORTRAN calls in eig.cc file.
> 
> The code in "eig.cc" file contains the logic to process the input
> arguments
> and determine if single or double precision is used. This file includes
> other header files such as "EIG.h" and "fEIG.h" (single precision
> computation). You can find them in the following route in the source code
> folder: "octave-3.XX/liboctave/numeric". The callings to the FORTRAN
> routines are In the source files "EIG.cc" and "fEIG.cc".
> 
>>> And current eig(A,B) function is giving compatible results with matlab.
> 
> Yes, the eig function is working. But there are some particular cases
> where
> Octave might give some roundoff errors. Matlab, on the other side,
> improves
> this problem by making a balancing procedure, which can also be done in
> Octave but it requires calling explicitly the "balance" function. You can
> see an example of this in the Matlab documentation of the "eig" function:
> Examples
> 
> The matrix
> 
> B = [ 3     -2      -.9    2*eps
>      -2      4       1    -eps
>      -eps/4  eps/2  -1     0
>      -.5    -.5      .1    1   ];
> 
> has elements on the order of roundoff error. It is an example for which
> the
> nobalance option is necessary to compute the eigenvectors correctly. Try
> the statements
> 
> [VB,DB] = eig(B)
> B*VB - VB*DB
> [VN,DN] = eig(B,'nobalance')
> B*VN - VN*DN
> 
> ---
> 
> In Octave, the first case will give a big roundoff error, which won't be
> the same as Matlab. The second case will yield an error, since such option
> is invalid. You may decide whether to use or not the balancing process,
> but
> an "eig" function that does it automatically should be very useful.
> 
> Also, if you compare the documentation of both functions you'll see that
> Matlab gives more options to choose from. I think it would be nice for
> Octave to implement the same, or at least a compatible form to call the
> eig
> functions, for compatibility. Even the same can be done with the "eigs"
> function, but that might not be desire for the Octave project, alias you
> can include it in your proposal.
> 
> The following is a table of the FORTRAN routines that are called with the
> different configuration of parameters:
> 
> 
>    *Case* *Routine*  Real symmetric A DSYEV  Real nonsymmetric A:
>  With preliminary balance step DGEEV (with SCLFAC =
> 2 instead of 8 in DGEBAL)  d = eig(A,'nobalance') DGEHRD, DHSEQR  [V,D] =
> eig(A,'nobalance') DGEHRD, DORGHR, DHSEQR, DTREVC  Hermitian A ZHEEV
> Non-Hermitian A:
>  With preliminary balance step ZGEEV (with SCLFAC =
> 2 instead of 8 in ZGEBAL)  d = eig(A,'nobalance') ZGEHRD, ZHSEQR  [V,D] =
> eig(A,'nobalance') ZGEHRD, ZUNGHR, ZHSEQR, ZTREVC   Real symmetric A,
> symmetric positive definite B. DSYGV      Special case:
>     eig(A,B,'qz') for real A, B
>     (same as real nonsymmetric A, real
>     general B) DGGEV  Real nonsymmetric A, real general B DGGEV  Complex
> Hermitian A,
> Hermitian positive definite B. ZHEGV      Special case:
>     eig(A,B,'qz') for complex A or B
>     (same as complex non-Hermitian A,
>     complex B) ZGGEV  Complex non-Hermitian A, complex B ZGGEV
> 
> One thing that you'll notice in the EIG.cc code is that there is no ZGEBAL
> nor DGEBAL, which are the FORTRAN balancing routines. They should be in
> the
> definition of the "balance" function but it should be desirable to perform
> the balancing as the table describes.
> 
> So, first things first, you must get familiar with the FORTRAN routines
> and
> its parameters. I don't think you even will need to code in FORTRAN, but
> calling the LAPACK libraries from C is a must. Then, look how to implement
> the proper calls in the Octave code as in the table, in order to give
> Octave the same calling scheme and functionality. And finally the testing
> and validation of your implementation including both Eigenvalue and
> Generalized Eigenvalue problems in similar cases as in the example I
> included. Perhaps, you'll also need to include some unit tests to assess
> the compliance of your implementation during compilation.
> 
> 
> 2015-04-15 7:55 GMT-05:00 András Mihálykó <

> mihalykoandras+levl@

> >:
> 
>> Hi,
>> I saw the project ideas list of ESA Summer of Code, and saw the
>> generalized
>> eigenvalue problem, however I didn't actually get, what should be done
>> with
>> it. I read this through:
>>
>> http://octave.1599824.n4.nabble.com/General-eigenvalue-problem-proposal-td4651990.html
>> but so far I saw, the generalized eigenvector problem is done, just
>> balancing isn't in the eig() function. Am I right or maybe I've
>> misunderstood something? Please, could you explain me, is this the main
>> focus of this project, or what else should be done with eig()?
>> Regards,
>> András Mihálykó
>>
>>
>>
>> --
>> View this message in context:
>> http://octave.1599824.n4.nabble.com/ESA-SoC-Generalized-eigenvalue-tp4669870.html
>> Sent from the Octave - Maintainers mailing list archive at Nabble.com.
>>
>>
> 
> 
> -- 
> Rafael González García

The Lapack folks have done some work on the balancing problem that the ML
example shows:

http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4270&p=10278&hilit=geev#p10278

This improved balancing function is in lapack 3.5.
They have not addressed the issue with *geev I reported to them (same link)
so that is
something to be aware of. Also if you do balancing before calling the
eigenvalue routine
instead of letting the eigenvalue routine handle it, you need to scale the
eigenvectors.
It seems to me the best option is to just use the *x eigenvalue routines
which have optional
balancing.



--
View this message in context: 
http://octave.1599824.n4.nabble.com/ESA-SoC-Generalized-eigenvalue-tp4669870p4669892.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.