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

 
 

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

Submitted by:  None
Submitted on:  Fri 11 Apr 2014 02:25:36 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: InvalidAssigned to: None
Originator Name: Russell CardenOriginator Email: -unavailable-
Open/Closed: ClosedRelease: 3.8.0
Operating System: GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission 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>
Project Administrator
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>
Project 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>
Project Administrator
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>
Project 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>
Project 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>
Project 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>
Project 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>
Project 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>
Project Administrator
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>
Project 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>
Project Administrator
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>
Project Administrator
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:

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

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
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> AV-VD
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-VD
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,

which is as expected.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
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> AV-VD
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.

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

Rik <rik5>
Project 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>
Project Administrator
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'AQ'), where Q is unitary all return different results. If I compute the eigenvectors as well
[V D]=eig(A)
then if I compute
AV-VD, 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 File(s):
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 5 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Thu 17 Apr 2014 02:20:35 PM UTCjordighStatusWorks For Me=>Invalid
      Open/ClosedOpen=>Closed
    Sat 12 Apr 2014 03:54:22 PM UTCrik5Release3.6.4=>3.8.0
    Fri 11 Apr 2014 03:39:15 PM UTCrik5StatusNeed Info=>Works For Me
    Fri 11 Apr 2014 03:31:44 PM UTCjordighStatusNone=>Need Info

    Back to the top


    Powered by Savane 3.1-cleanup1