[Octave-bug-tracker] [bug #42097] 64bit Slackware, eig returning incorrect results

[Russell Carden]

Follow-up Comment #22, bug #42097 (project octave):

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








    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?42097>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/