Re: [ESA-SoC] Generalized eigenvalue

[Rafael Gonzalez]

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ó <address@hidden>:

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