bugGNU Octave - Bugs: bug #42097, 64bit Slackware,eig returning...

 
 

bug #42097: 64bit Slackware,eig returning incorrect results

Submitter:  None
Submitted:  Fri 11 Apr 2014 02:25:36 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Invalid / Not an Octave Bug Assigned to:  None
Originator Name:  Russell Carden Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.8.0
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 18 Apr 2014 05:24:44 PM UTC, comment #32: 

Looking more closely, that must have been what happenned.  It was the plain arapack that was causing the difficulty.  Perhaps, I compiled things out of order: compiled arpack while either lapack and/or blas were not installed.  The standard arpack has options to include and it checks for blas and lapack.  When removing things to find out what was wrong, I had overlooked arpack, otherwise I would have found my problem sooner and not have bother you all.  Thanks for your help.  I'm glad to have my confidence octave's eig's restored.

Russell Carden <russell_l_carden>
Thu 17 Apr 2014 02:20:35 PM UTC, comment #31: 

arpack-ng removed the duplicate and apparently buggy BLAS libraries a long time:

http://forge.scilab.org/index.php/p/arpack-ng/source/commit/c9d9d92d5108fea1ba617684d55cd8faff1fd146/

Are you sure you're getting a modified dlahqr from arpack-ng and not from the old arpack?

At any rate, can we conclude at this point that this is not a bug in Octave, but in the libraries linked in to Octave? I think so. I'm closing this bug as invalid. If you disagree and think the bug is in Octave, we can reopen it.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Thu 17 Apr 2014 01:35:43 PM UTC, comment #30: 

Is ARPACK providing a BLAS function?  I don't think it should be doing that.  It should just be using the BLAS library.


John W. Eaton <jwe>
Group administrator
Thu 17 Apr 2014 01:21:01 PM UTC, comment #29: 

This is weird. ARPACK should not be called at all for eig, only for eigs.

What is the bug in arpack-ng? Can you please report it there?

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Wed 16 Apr 2014 09:23:27 PM UTC, comment #28: 

So both the test program and octave produce the same results up
to line 466 of dhseqr.f.
466  CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILO,
467          $                   IHI, Z, LDZ, INFO )
At this point I am able to step into DLAHQR in my test program, but in gdb for octave I step right over this function call.  After this function call the results from octave are incorrect. 
After some googling, I found that I should tell gdb
set step-mode on
I then found that that there is a dlahqr_ that is used in arpack
and this is the dlahqr_ that octave is calling.
Note, the version of arpack that I installed is arpack-ng.  Supposedly arpack-ng is actively maintained while the regular arpack is not.  When I saw this, I removed arpack-ng and installed the regular arpack, and now I have the correct results.
Apparently the regular arpack build script has an important flag/option that arpack-ng is missing.

Russell Carden <russell_l_carden>
Wed 16 Apr 2014 12:47:09 PM UTC, comment #27: 

I have blas, lapack and octave compiled for debugging.  Everything looks good up until the schur factorization is computed.  That is the balancing and the transformation to upper Hessenberg form all produce the same results.  Things diverge after the schur factorization.  I would have more info but for some reason when I reached the call for the schur factorization, I decided I should setup blas for debugging.  Hope to have more infor this afternoon.


Russell Carden <russell_l_carden>
Tue 15 Apr 2014 06:51:06 PM UTC, comment #26: 

The next thing I would try is to run both your small test program and Octave with gdb and then step in to the dgeev call and ask for a stack trace.  Then I think you should be able to verify that they are both actually using the same LAPACK library.  You could also set a breakpoint for one of the BLAS library functions that DGEEV calls and then also verify that they are both using the same BLAS library.

Beyond that, if you compiled the BLAS and LAPACK libraries yourself, then you can step through the two programs side by side and find where the computation goes wrong in one but not the other.  That might give some clues about what is going wrong.

John W. Eaton <jwe>
Group administrator
Tue 15 Apr 2014 06:31:25 PM UTC, comment #25: 

I removed -fPIC and -O2 from the compiler options in the build script.  No change, eig still returns incorect results for the test matrix, A=[0 1 0 ; 0 0 1 ; 1 0 0 ].

Russell Carden <russell_l_carden>
Tue 15 Apr 2014 02:04:18 PM UTC, comment #24: 

Sometimes we've seen compiler options have an effect.  I would take out the '-fPIC' flag option from the build script for the x86_64 architecture (yours).  I let configure/make figure out which files want -fPIC.

Rik <rik5>
Group administrator
Tue 15 Apr 2014 02:21:05 AM UTC, comment #23: 

I tried one last thing.  I ripped out the F77 macros and
called the dgeev in the same manner that I call it in the test program, and I still get the wrong results from octave.
Is there a compiler option/flag that could be causing the problem?

Russell Carden <russell_l_carden>
Mon 14 Apr 2014 08:10:58 PM UTC, comment #22: 

Yes, the lapack functions for general matrices work just fine for symmetric matrices.  This is true of lapack.  However, on my machine for whatever reason, when octave calls lapack it doesn't work for matrices that are not symmetric and not upper hessenberg.
Entries below the first diagonal seem to lead to perturbations in other parts of the matrix. 

Note, I have observed for a couple test matrices that
if I do the following
[V D]=eig(A*2^j)
D=D/2^j
for j greater than number of bits in double precision, 53
then V D form an accurate to machine precision eigenvector eigeinvalue factorization of A.  Another, words whatever is wacky in the call that octave makes to lapack,  the perturbation does not depend on the norm of A.


I have compared the inputs between of the fortran call to dgeev
and octave's call to dgeev.  The inputs are the same.  I even redid my test program in C.
Octave's call to dgeev is garbled with F77_* function calls.
Does anyone know how these work, so I might be able to determine exactly how dgeev is called

  F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
                           F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
                           n, tmp_data, n, pwr, pwi, dummy,
                           idummy, pvr, n, &dummy_work, lwork, info
                           F77_CHAR_ARG_LEN (1)
                           F77_CHAR_ARG_LEN (1)));

I compiled octave with -g and no optimization.  beyond the F77_XFCN, gdb only picks up a call in misc/cquit.c to
octave_restore_current_context.  Beyond that call I don't know what is going on. 


Input Matrix  and WORK after DGEEV in octave
Note, I gave the entries of WORK the values 1:sizeof(WORK),
so I could see exactly which entries were altered.
Note the indexing differes between fortran and C/octave.

after A(0,0)==-0.341164
after A(1,0)==-0.862241
after A(2,0)==0.000000
after A(0,1)==1.564735
after A(1,1)==-0.341164
after A(2,1)==0.000000
after A(0,2)==-0.083729
after A(1,2)==0.320588
after A(2,2)==0.682328
after Work(0)==102.000000
after Work(1)==1.000000
after Work(2)==1.000000
after Work(3)==3.000000
after Work(4)==0.000000
after Work(5)==6.000000
after Work(6)==1.000000
after Work(7)==0.000000
after Work(8)==0.000000
after Work(9)==10.000000
after Work(10)==11.000000
after Work(11)==12.000000
after Work(12)==13.000000
after Work(13)==14.000000
after Work(14)==15.000000
after Work(15)==16.000000
after Work(16)==17.000000
after Work(17)==18.000000
after Work(18)==19.000000
after Work(19)==20.000000

Matrix and WORK after dgeev in test program
A after computing eigenvalues
A(1,1)=  -0.500
A(2,1)=  -0.866
A(3,1)=   0.000
A(1,2)=   0.866
A(2,2)=  -0.500
A(3,2)=   0.000
A(1,3)=   0.000
A(2,3)=   0.000
A(3,3)=   1.000
Work after computing eigenvalues
WORK(  1)= 102.000
WORK(  2)=   1.000
WORK(  3)=   1.000
WORK(  4)=   3.000
WORK(  5)=   0.000
WORK(  6)=   6.000
WORK(  7)=   1.000
WORK(  8)=   1.000  <--- 0 in octave results
WORK(  9)=   9.000  <--- 0 in octave results
WORK( 10)=  10.000
WORK( 11)=  11.000
WORK( 12)=  12.000
WORK( 13)=  13.000
WORK( 14)=  14.000
WORK( 15)=  15.000
WORK( 16)=  16.000
WORK( 17)=  17.000
WORK( 18)=  18.000
WORK( 19)=  19.000
WORK( 20)=  20.000







Russell Carden <russell_l_carden>
Mon 14 Apr 2014 03:32:33 PM UTC, comment #21: 

I would say that Octave could call a different function for symmetric matrices because it might be more efficient to do so, but the lapack function for general matrices should also work if the matrix happens to be symmetric, shouldn't it?

John W. Eaton <jwe>
Group administrator
Mon 14 Apr 2014 02:50:35 PM UTC, comment #20: 

Ignore my remark about the 4x4 example for which eig(A^2) returned the correct results.  A^2 happens to be symmetric, in which case octave should call a different function.

Russell Carden <russell_l_carden>
Mon 14 Apr 2014 02:47:12 PM UTC, comment #19: 

I almost have all the code in place to do what you are suggesting.
I will do that later this afternoon.
However, let me note The results are incorrect immediately after the call to DGEEV.  The eigenvalues and the eigenvectors are incorrect.  If you compute the perturbation to A, as I did in a previous comment,  you find that the presenece of the A(3,1) entry results in a perturbation to the A(2,1) entry of the matrix.  This perturbation is not a linear function of the A(3,1) entry.
I just did a test on a similarly structured matrix
A=diag(ones(3,1),1), A(4,1)=1.  eig(A) returns incorrect results,
However when I apply eig(A^2) which is not an upper Hessenberg matrix, I get correct results.  Octave always seems to give the right answer for upper Hessenberg matrices.


Russell Carden <russell_l_carden>
Mon 14 Apr 2014 02:23:55 PM UTC, comment #18: 

I would debug this problem by finding out precisely what arguments Octave is passing to DGEEV and what values it returns before Octave post-processes the results to return them in the form of diagonal matrices, etc.  Then I would create a separate program that passes exactly the same arguments to DGEEV but that does not depend on Octave in any way.  I would expect the separate program to show the same problem.  If it does, the problem is in the lapack or blas libraries.  If it does not, then I'd have to determine why there is a difference.  Maybe the problem is in the way Octave is post-processing the results returned from DGEEV.

For the record, I'm not able to reproduce the results you shown in comment #5.  I get the same results as shown in comment #6.

John W. Eaton <jwe>
Group administrator
Mon 14 Apr 2014 01:55:05 PM UTC, comment #17: 

Actually src/octave-cli appears to be a script but the script mentions src/.lib/octave-cli.
Applying ldd to this file returns along with many other lines

liblapack.so.3 => /usr/lib64/../lib64/liblapack.so.3
(0x00007f95d94d7000)
        libblas.so.3 => /usr/lib64/../lib64/libblas.so.3 (0x00007f95d9281000)
        libreadline.so.5 => /usr/lib64/../lib64/libreadline.so.5 (0x00007f95d9042000)
        libncurses.so.5 => /usr/lib64/../lib64/libncurses.so.5 (0x00007f95d8dec000)
        libpcre.so.1 => /usr/lib64/../lib64/libpcre.so.1 (0x00007f95d8baa000)
        libdl.so.2 => /lib64/libdl.so.2 (0x00007f95d89a5000)
        libgfortran.so.3 => /usr/lib64/../lib64/libgfortran.so.3 (0x00007f95d868e000)

So indicates that octave is linking to the appropriate files, right?

Russell Carden <russell_l_carden>
Mon 14 Apr 2014 01:38:27 PM UTC, comment #16: 

The file src/octave is a simple wrapper program that execs the real Octave process.  Try running ldd on src/octave-cli instead.

John W. Eaton <jwe>
Group administrator
Mon 14 Apr 2014 01:23:42 PM UTC, comment #15: 

First of all ZGEEV is not the LAPACK routine I should be using for my test matrix, since my test matrix is real.
So I updated my test program to use DGEEV.  I call DGEEV
twice, once to compute the optimal size for the variable work
and then again to do the actual eigenvalue computation.
The same operations are performed in octave
Again the test program returns the correct results.
Octave does not.  I added printf statements to EIG.cc
to see the inputs in the calls to DGEEV that octave makes
I have looked over the important inputs, and they are correct and they are the same as the ones in the test program.  I will look over the other variables sometime this afternoon.

The function call in EIG.cc is

      F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
                               F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
                               n, btmp_data, n, pwr, pwi, dummy,
                               idummy, pvr, n, pwork, lwork, info
                               F77_CHAR_ARG_LEN (1)
                               F77_CHAR_ARG_LEN (1)));

I see the mention of F77 in the octave code.  Note I compiled the test program as a .f90 not a .f77.  Might that have some effect.  (I am not familiar with FOTRAN, I just know that gfortran gave me strange warning until I gave my test program the right extension)

See the ldd outputs below, there is no menition of lapack or blas from the octave binary.  Am I using ldd correctly.
I am applying it to the right file.  When I apply to liboctave.so, I see mention of the correct lapack and blas files.
Is something going wrong with the linking.  If so, where do I need to look.  If not, why does it work at all?


Output from ldd for test program
ldd ./a.out
        linux-vdso.so.1 (0x00007fffbe379000)
        liblapack.so.3 => /usr/lib64/liblapack.so.3 (0x00007f03c49f2000)
        libblas.so.3 => /usr/lib64/libblas.so.3 (0x00007f03c479c000)
        libgfortran.so.3 => /usr/lib64/libgfortran.so.3 (0x00007f03c4484000)
        libm.so.6 => /lib64/libm.so.6 (0x00007f03c4182000)
        libgcc_s.so.1 => /usr/lib64/libgcc_s.so.1 (0x00007f03c3f6c000)
        libquadmath.so.0 => /usr/lib64/libquadmath.so.0 (0x00007f03c3d30000)
        libc.so.6 => /lib64/libc.so.6 (0x00007f03c3967000)
        /lib64/ld-linux-x86-64.so.2 (0x00007f03c527f000)
Output of ldd for ./src/octave
ldd ./src/octave                   
        linux-vdso.so.1 (0x00007fff527ff000)
        libX11.so.6 => /usr/lib64/libX11.so.6 (0x00007fd28c800000)
        libxcb.so.1 => /usr/lib64/libxcb.so.1 (0x00007fd28c5e3000)
        libXau.so.6 => /usr/lib64/libXau.so.6 (0x00007fd28c3df000)
        libXdmcp.so.6 => /usr/lib64/libXdmcp.so.6 (0x00007fd28c1da000)
        libdl.so.2 => /lib64/libdl.so.2 (0x00007fd28bfd6000)
        libutil.so.1 => /lib64/libutil.so.1 (0x00007fd28bdd2000)
        libstdc++.so.6 => /usr/lib64/libstdc++.so.6 (0x00007fd28bad0000)
        libm.so.6 => /lib64/libm.so.6 (0x00007fd28b7ce000)
        libgcc_s.so.1 => /usr/lib64/libgcc_s.so.1 (0x00007fd28b5b7000)
        libpthread.so.0 => /lib64/libpthread.so.0 (0x00007fd28b39a000)
        libc.so.6 => /lib64/libc.so.6 (0x00007fd28afd1000)
        /lib64/ld-linux-x86-64.so.2 (0x00007fd28cb7b000)


Russell Carden <russell_l_carden>
Sat 12 Apr 2014 07:59:08 PM UTC, comment #14: 

Could you grep -i Octave's source code for ZGEEV and add some printf-debugging commands in there to see if you're getting the back same values as what you got from your standalone program? Or if you're knowledgeable enough, can you put a breakpoint at that location and run Octave under gdb? Also check with ldd that your standalone program and the Octave binary are linked to the same libraries.

I'm sorry, but you really are going to have to do your own debugging here. I am not yet convinced that the problem does not lie with the libraries linked to Octave in your build.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Sat 12 Apr 2014 07:43:44 PM UTC, comment #13: 

I have not altered the build scripts from slackbuild.org.
I use them as is.  As I said in a previous comment, I tested out LAPACK's ZGEEV, and it returns the right answer. 

uname -m      
returns                                   
x86_64                                                                                                                      
gcc --version                                      
returns
gcc (GCC) 4.8.2                                                                                                             
Copyright (C) 2013 Free Software Foundation, Inc.                                                                           
This is free software; see the source for copying conditions.  There is NO                                                  
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                                                 

Russell Carden <russell_l_carden>
Sat 12 Apr 2014 03:54:22 PM UTC, comment #12: 

I took a look at the Octave build script from buildscripts.org and it looks pretty simple.  What does 'uname -m' return for you?  Are you running the script as is?  Did you modify it to throw any other options to configure or change the values for CFLAGS, CXXFLAGS, FFLAGS?  What does 'gcc --version' return?

On a working Ubuntu distribution the BLAS package is libblas3gf (1.2.20110419-2ubuntu1).  I don't know what patches they may have added, but it appears they at least started with the 20110419 version.  The LAPACK version is 3.3.1-1.

Getting the right BLAS/LAPACK combination has been an issue in the past and I'm guessing that is still where the problem lies.

Rik <rik5>
Group administrator
Sat 12 Apr 2014 02:05:23 PM UTC, comment #11: 

I have the version of BLAS from netlib.org.
Slackbuilds.org has time stamp on it of 20110419.
This timestamp is what is listed as the version on slackbuilds.org.
I have used reference blas in the past on other systems( not 64bit)
without any problems.

No, I did not have two different versions of lapack installed at the same time.  I removed the version of lapack for atlas before installing the version of lapack from netlib.

Russell Carden <russell_l_carden>
Sat 12 Apr 2014 01:33:00 PM UTC, comment #10: 


> I used the lapack and blas versions from slackbuilds.org
> LAPACK version 3.4.1
> BLAS version 3.2.2


What does this mean? BLAS is an API that has many implementations, e.g. ATLAS, OpenBLAS, Intel MKL, Apple's Accelerate framework, as well as the old Fortran reference BLAS, which has been frozen since the 70's.

I can't find out which BLAS is 3.2.2. The closest I can guess is an old version of LAPACK that had that version number.

Did you really install two different LAPACK versions at the same time?

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Sat 12 Apr 2014 01:44:39 AM UTC, comment #9: 

I used the lapack and blas versions from slackbuilds.org
LAPACK version 3.4.1
BLAS version 3.2.2

I tested the same matrix with LAPACK in FORTRAN,
and I get the correct results.  So the problem has to be
somewhere else. 

Eigenvalues
 (-0.50000000000000011     , 0.86602540378443904     )
 (-0.50000000000000011     ,-0.86602540378443893     )
 (  1.0000000000000000     ,  0.0000000000000000     )
First Eigenvector
 (-0.28867513459481287     ,-0.49999999999999989     )
 ( 0.57735026918962606     ,  0.0000000000000000     )
 (-0.28867513459481331     , 0.49999999999999978     )
Second Eigenvector
 (-0.28867513459481270     , 0.49999999999999989     )
 ( 0.57735026918962584     ,  0.0000000000000000     )
 (-0.28867513459481320     ,-0.49999999999999967     )
Third Eigenvector 
 ( 0.57735026918962562     , -1.6177811532852776E-032)
 ( 0.57735026918962573     , -1.2325951644078307E-032)
 ( 0.57735026918962584     ,  0.0000000000000000     )

Russell Carden <russell_l_carden>
Fri 11 Apr 2014 10:17:16 PM UTC, comment #8: 

Which LAPACK version? Which BLAS implementation?

The bug is almost certainly in those two. Octave doesn't actually perform any significant part of the eigenvalue computations itself.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Fri 11 Apr 2014 07:09:30 PM UTC, comment #7: 

I compiled octave, blas, and lapack using the scripts on slackbuild.org for version 14.1 of slackware.  I had been using using atlas until I encountered this problem with eig.  Since encountering this problem with eig, I removed atlas and the packages affected by its removal and then reinstalled.

Everything seems to work for upper Hessenberg matrices.  Scaling the matrix, computing the eigenvalues and then scaling back produces incorrect results.

For the matrix from my previous comment, using the V and D returned from eig

octave:45> V*D/V
ans =

 Columns 1 and 2:

   0.000000000000001 - 0.000000000000000i   1.000000000000000 + 0.000000000000000i
  -1.000000000000001 + 0.000000000000000i  -0.000000000000001 - 0.000000000000000i
   1.000000000000000 - 0.000000000000000i   0.000000000000000 + 0.000000000000000i

 Column 3:

  -0.000000000000000 + 0.000000000000000i
   1.000000000000000 - 0.000000000000000i
  -0.000000000000000 - 0.000000000000000i

or V*D/V=
  0 1 0
 -1 0 0
  1 0 0

The entry (3,1) entry seems to result in a perturbation of the (2,1) entry.  The perturbation is not a linear function of A(3,1).  See the following

octave:115> for i=[1:38], A(3,1)=i; [V D]=eig(A);Ahat=V*D/V; disp([i Ahat(2,1), Ahat(3,1)]), end
   1.0000 + 0.0000i  -1.0000 + 0.0000i   1.0000 + 0.0000i
   2.0000 + 0.0000i  -1.0000 - 0.0000i   2.0000 + 0.0000i
   3.0000 + 0.0000i  -2.0000 + 0.0000i   3.0000 + 0.0000i
   4.0000 + 0.0000i  -2.0000 - 0.0000i   4.0000 - 0.0000i
   5.0000 + 0.0000i  -2.0000 - 0.0000i   5.0000 + 0.0000i
   6.0000 + 0.0000i  -2.0000 + 0.0000i   6.0000 - 0.0000i
   7.0000 + 0.0000i  -2.0000 + 0.0000i   7.0000 + 0.0000i
   8.0000 + 0.0000i  -2.0000 + 0.0000i   8.0000 - 0.0000i
   9.0000 + 0.0000i  -2.0000 + 0.0000i   9.0000 + 0.0000i
   10.0000 +  0.0000i   -2.0000 -  0.0000i   10.0000 +  0.0000i
   11   -2   11
   12.0000 +  0.0000i   -2.0000 +  0.0000i   12.0000 +  0.0000i
   13.0000 +  0.0000i   -2.0000 +  0.0000i   13.0000 +  0.0000i
   14.0000 +  0.0000i   -2.0000 +  0.0000i   14.0000 -  0.0000i
   15.0000 +  0.0000i   -2.0000 +  0.0000i   15.0000 +  0.0000i
   16.0000 +  0.0000i   -2.0000 +  0.0000i   16.0000 +  0.0000i
   17.0000 +  0.0000i   -2.0000 +  0.0000i   17.0000 +  0.0000i
   18.0000 +  0.0000i   -2.0000 -  0.0000i   18.0000 -  0.0000i
   19.0000 +  0.0000i   -2.0000 +  0.0000i   19.0000 +  0.0000i
   20.0000 +  0.0000i   -2.0000 -  0.0000i   20.0000 +  0.0000i
   21.0000 +  0.0000i   -2.0000 +  0.0000i   21.0000 -  0.0000i
   22.0000 +  0.0000i   -2.0000 -  0.0000i   22.0000 +  0.0000i
   23.0000 +  0.0000i   -2.0000 -  0.0000i   23.0000 -  0.0000i
   24.0000 +  0.0000i   -2.0000 +  0.0000i   24.0000 -  0.0000i
   25.0000 +  0.0000i   -2.0000 -  0.0000i   25.0000 +  0.0000i
   26.0000 +  0.0000i   -2.0000 +  0.0000i   26.0000 +  0.0000i
   27.0000 +  0.0000i   -2.0000 -  0.0000i   27.0000 -  0.0000i
   28.0000 +  0.0000i   -2.0000 -  0.0000i   28.0000 +  0.0000i
   29.0000 +  0.0000i   -2.0000 +  0.0000i   29.0000 +  0.0000i
   30.0000 +  0.0000i   -2.0000 +  0.0000i   30.0000 -  0.0000i
   31.0000 +  0.0000i   -2.0000 -  0.0000i   31.0000 -  0.0000i
   32.0000 +  0.0000i   -2.0000 +  0.0000i   32.0000 +  0.0000i
   33.0000 +  0.0000i   -2.0000 +  0.0000i   33.0000 +  0.0000i
   34.0000 +  0.0000i   -2.0000 +  0.0000i   34.0000 +  0.0000i
   35.0000 +  0.0000i   -2.0000 +  0.0000i   35.0000 +  0.0000i
   36.0000 +  0.0000i   -2.0000 -  0.0000i   36.0000 +  0.0000i
   37.0000 +  0.0000i   -2.0000 +  0.0000i   37.0000 +  0.0000i
   38.0000 +  0.0000i   -4.0000 -  0.0000i   38.0000 +  0.0000i


Russell Carden <russell_l_carden>
Fri 11 Apr 2014 06:13:10 PM UTC, comment #6: 

Something is awfully wrong on your side. These are my results:


octave:7> A = [0 1 0; 0 0 1; 1 0 0]
A =

   0   1   0
   0   0   1
   1   0   0

octave:8> format long
octave:9> [V, D] = eig(A)
V =

   0.577350269189626 + 0.000000000000000i   0.577350269189626 - 0.000000000000000i  -0.577350269189626 + 0.000000000000000i
  -0.288675134594813 + 0.500000000000000i  -0.288675134594813 - 0.500000000000000i  -0.577350269189626 + 0.000000000000000i
  -0.288675134594813 - 0.500000000000000i  -0.288675134594813 + 0.500000000000000i  -0.577350269189626 + 0.000000000000000i

D =

Diagonal Matrix

  -0.500000000000000 + 0.866025403784439i                                        0                                        0
                                        0  -0.500000000000000 - 0.866025403784439i                                        0
                                        0                                        0   1.000000000000000 + 0.000000000000000i

octave:10> A*V -V*D
ans =

  -5.55111512312578e-17 - 6.10622663543836e-16i  -5.55111512312578e-17 + 6.10622663543836e-16i  -1.11022302462516e-16 + 0.00000000000000e+00i
  -1.11022302462516e-16 + 2.22044604925031e-16i  -1.11022302462516e-16 - 2.22044604925031e-16i  -4.44089209850063e-16 + 0.00000000000000e+00i
   2.22044604925031e-16 + 1.11022302462516e-16i   2.22044604925031e-16 - 1.11022302462516e-16i   2.22044604925031e-16 + 0.00000000000000e+00i


How did you compile Octave? Which implementation of LAPACK are you using?

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Fri 11 Apr 2014 05:06:26 PM UTC, comment #5: 

Here is a cleaner example.  Note that for this particular matrix,
a correct eigenvalue decomposition is returned for A' but not for A.  This is probably be due to the fact that A' is already upper Hessenberg, so that the problem may lie in bringing A to upper Hessenberg form.

octave:1> A=[ 0 1 0 ; 0 0 1 ; 1  0 0 ]
A =

   0   1   0
   0   0   1
   1   0   0

octave:2> format long
octave:3> [V D]=eig(A)
V =

 Columns 1 and 2:

  -0.158836098085991 - 0.540780260472586i  -0.158836098085991 + 0.540780260472586i
   0.682327803828019 + 0.000000000000000i   0.682327803828019 - 0.000000000000000i
  -0.391621714024374 + 0.251771732042861i  -0.391621714024374 - 0.251771732042861i

 Column 3:

  -0.526062999881448 + 0.000000000000000i
  -0.358947411384288 + 0.000000000000000i
  -0.770982798781043 + 0.000000000000000i

D =

Diagonal Matrix

 Columns 1 and 2:

  -0.341163901914010 + 1.161541399997252i                                        0
                                        0  -0.341163901914010 - 1.161541399997252i
                                        0                                        0

 Column 3:

                                        0
                                        0
   0.682327803828019 + 0.000000000000000i

octave:4> A*V-V*D
ans =

 Columns 1 and 2:

  -0.000000000000000 + 0.000000000000000i  -0.000000000000000 - 0.000000000000000i
  -0.158836098085990 - 0.540780260472587i  -0.158836098085990 + 0.540780260472587i
  -0.000000000000000 + 0.000000000000000i  -0.000000000000000 - 0.000000000000000i

 Column 3:

   0.000000000000000 + 0.000000000000000i
  -0.526062999881449 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i

octave:5> [V D]=eig(A')
V =

 Columns 1 and 2:

   0.577350269189626 + 0.000000000000000i   0.577350269189626 - 0.000000000000000i
  -0.288675134594813 - 0.500000000000000i  -0.288675134594813 + 0.500000000000000i
  -0.288675134594813 + 0.500000000000000i  -0.288675134594813 - 0.500000000000000i

 Column 3:

  -0.577350269189626 + 0.000000000000000i
  -0.577350269189626 + 0.000000000000000i
  -0.577350269189626 + 0.000000000000000i

D =

Diagonal Matrix

 Columns 1 and 2:

  -0.500000000000000 + 0.866025403784439i                                        0
                                        0  -0.500000000000000 - 0.866025403784439i
                                        0                                        0

 Column 3:

                                        0
                                        0
   1.000000000000000 + 0.000000000000000i

octave:6> A
A =

   0   1   0
   0   0   1
   1   0   0

octave:7> A'*V-V*D
ans =

 Column 1:

  -2.22044604925031e-16 - 5.55111512312578e-17i
   1.11022302462516e-16 + 2.77555756156289e-17i
  -1.11022302462516e-16 + 3.33066907387547e-16i

 Column 2:

  -2.22044604925031e-16 + 5.55111512312578e-17i
   1.11022302462516e-16 - 2.77555756156289e-17i
  -1.11022302462516e-16 - 3.33066907387547e-16i

 Column 3:

  -2.22044604925031e-16 + 0.00000000000000e+00i
   4.44089209850063e-16 + 0.00000000000000e+00i
   1.11022302462516e-16 + 0.00000000000000e+00i

octave:8>

Russell Carden <russell_l_carden>
Fri 11 Apr 2014 04:33:02 PM UTC, comment #4: 

Can you use format long and come up with an example again? When I attempt your example with the precision you provided,


octave:10> A = [   1.53587   0.29450   0.15366
>    0.60792  -1.93114   0.67947
>   -0.31746   1.29087  -1.14817]
A =

   1.53587   0.29450   0.15366
   0.60792  -1.93114   0.67947
  -0.31746   1.29087  -1.14817

octave:11> A
A =

   1.53587   0.29450   0.15366
   0.60792  -1.93114   0.67947
  -0.31746   1.29087  -1.14817

octave:12> [V, D] = eig(A)
V =

   0.985834   0.027681  -0.122808
   0.163509  -0.738089   0.393531
  -0.037364   0.674136   0.911072

D =

Diagonal Matrix

   1.57889         0         0
         0  -2.57454         0
         0         0  -0.54780

octave:13> A*V - V*D
ans =

  -1.1102e-15   2.2204e-16   5.5511e-17
   1.1102e-16   1.1102e-15  -3.0531e-16
   2.0817e-17   0.0000e+00   5.5511e-16


which is as expected.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Fri 11 Apr 2014 04:18:40 PM UTC, comment #3: 

Thanks for the quick reply.
Here is an example of what I am observing.
octave:1> version
ans = 3.8.0
octave:2> A=randn(3)
A =

   1.53587   0.29450   0.15366
   0.60792  -1.93114   0.67947
  -0.31746   1.29087  -1.14817

octave:3> [V D]=eig(A)
V =

   0.983566  -0.124690  -0.027442
   0.112137   0.393865   0.735934
  -0.141504   0.910671  -0.676497

D =

Diagonal Matrix

   1.54734         0         0
         0  -0.51662         0
         0         0  -2.57417

octave:4> A*V-V*D
ans =

   2.2204e-16  -6.9389e-17  -6.9389e-17
   1.1172e-01  -1.4163e-02  -3.1170e-03
   2.1394e-01  -2.7121e-02  -5.9689e-03

octave:5> eig(A')
ans =

   1.55218
  -2.53387
  -0.56175



Russell Carden <russell_l_carden>
Fri 11 Apr 2014 03:39:15 PM UTC, comment #2: 

I think this is probably specific to your system.  First, the 3.6.X series is deprecated an no longer supported.  It would be better to use the current stable release which is 3.8.X.

I used the following code as a test case which I infer is similar to what you did.


A = rand (3, 3);
v1 = eig (A);
v2 = eig (A');
v1 - v2
ans =

  -1.5543e-15
  -2.2204e-16
   5.5511e-17


This is generally within a few eps() of zero so is as good as the machine can do.


Rik <rik5>
Group administrator
Fri 11 Apr 2014 03:31:44 PM UTC, comment #1: 

O(1) doesn't make sense here. Of course everything is bounded by a constant in this context. What do you really mean? Do you mean tiny nonzero numbers?

Please show the actual example that is giving you trouble, but without further context, http://wiki.octave.org/FAQ#Why_is_this_floating_point_computation_wrong.3F

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Fri 11 Apr 2014 02:25:36 PM UTC, original submission:  

I am running Slackware 14.1 64bit.
I compiled octave using the scripts at slackbuilds.org.
I found that eig returns incorrect results for matrices with n>=3.
I found eig(A) and eig(A') as well as eig(Q'*A*Q'), where Q is unitary all return different results.  If I compute the eigenvectors as well
[V D]=eig(A)
then if I compute
A*V-V*D, I get a matrix with
zeros in the first row and O(1) entries in the remaining rows.
I noticed that eig code for make check only tests matrices of size 2, which is why this problem wasn't detected at comiple time.

Has anyone else encountered this problem?  If not I would be interested in helping to figure out what is wrong.

Thank you

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by russell_l_carden (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-04-17 jordigh StatusWorks For Me Invalid / Not an Octave Bug
        Open/ClosedOpen Closed
    2014-04-12 rik5 Release3.6.4 3.8.0
    2014-04-11 rik5 StatusNeed Info Works For Me
    2014-04-11 jordigh StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code