From help-octave-request at bevo dot che dot wisc dot edu  Mon Dec 21 09:30:45 1998
Subject: Re: Help...
From: lash at tellabs dot com
To: livneo at wisdom dot weizmann dot ac dot il (Livne Oren)
Cc: help-octave at bevo dot che dot wisc dot edu
Date: Mon, 21 Dec 1998 09:29:59 -0600 (CST)

> 
> Dear Guys,
> 
> I am in a great need for eigenvalue package in my Ph.D. studies. I've
> used matlab for this package ([V,D]=eig(a)) and it works, but I am
> interested in 32 digits precision whereas matlab has 16 only. Hence I
> decided to write my own routine. I got some programs (zgeev, qzhes, etc.)
> in fortran but they are all unreliable. Octave uses them too and here's an
> example that it doesn't work: (I am interested in general complex 4x4 and
> 8x8 matrices, full eigenvalues & eigenvectors):
> 
> a =
> 
>  Column 1:
> 
>   0.660674095153809 + 0.304533690214157i
>   0.218023106455803 + 0.675618231296539i
>   0.258169531822205 + 0.914756238460541i
>   0.301558434963226 + 0.222921922802925i
> 
>  Column 2:
> 
>   0.283803135156631 + 0.052084900438786i
>   0.474917590618134 + 0.666901826858521i
>   0.590481042861938 + 0.390818536281586i
>   0.554883301258087 + 0.487036257982254i
> 
>  Column 3:
> 
>   0.548641562461853 + 0.745089948177338i
>   0.714440703392029 + 0.638778746128082i
>   0.148456916213036 + 0.686463475227356i
>   0.106128662824631 + 0.144878074526787i
> 
>  Column 4:
> 
>   0.144677996635437 + 0.371638923883438i
>   0.906892597675323 + 0.958721995353699i
>   0.942158818244934 + 0.547215640544891i
>   0.406394332647324 + 0.574961066246033i
> 
> octave:66> eig(a)
> error: zgeev failed to converge
> error: evaluating index expression near line 66, column 1
> octave:66> 
> 
> 
> 
> Since matlab works for any matrix I suspect the routines I have and which
> octave uses are fine but need some driver of balancing the matrix. Can
> any of you guys kindly help me? I would be of course interested in source
> in C, but even fortran routines are fine. Maybe it just needs a call to
> some balancing routine before zgeev, otherwise no convergence for the
> eigenvectors. I have this problem already for many months but didn't find
> yet ansewr :( I could make a lot of progress in the Ph.D. instead of
> calling the slow matlab interface I wrote and/or go to the high precision.
> 
> If you find an answer please return to me at livneo at wisdom dot weizmann dot ac dot il
> - Thanks! Oren.
> 

Which version of Octave, and on what platform?  I just tried your example
on Octave 2.0.13 on a Sun Ultra 1 under solaris, and I get the results:

octave:11> version
ans = 2.0.13
octave:12> a
a =

 Column 1:

  0.660674095153809 + 0.304533690214157i
  0.218023106455803 + 0.675618231296539i
  0.258169531822205 + 0.914756238460541i
  0.301558434963226 + 0.222921922802925i

 Column 2:

  0.283803135156631 + 0.052084900438786i
  0.474917590618134 + 0.666901826858521i
  0.590481042861938 + 0.390818536281586i
  0.554883301258087 + 0.487036257982254i

 Column 3:

  0.548641562461853 + 0.745089948177338i
  0.714440703392029 + 0.638778746128082i
  0.148456916213036 + 0.686463475227356i
  0.106128662824631 + 0.144878074526787i

 Column 4:

  0.144677996635437 + 0.371638923883438i
  0.906892597675323 + 0.958721995353699i
  0.942158818244934 + 0.547215640544891i
  0.406394332647324 + 0.574961066246033i

octave:13> eig(a)
ans =

   1.806166293907551 + 2.030664098006237i
  -0.136299627666124 - 0.401122045040072i
   0.169073650200598 + 0.273930811396951i
  -0.148497381809724 + 0.329387194182949i

octave:14> [v,d]=eig(a)
v =

 Column 1:

   0.375114592535907 - 0.068746998143506i
   0.655757003642081 + 0.000000000000000i
   0.546031202268586 + 0.003354453649958i
   0.353494938156497 - 0.037757745891452i

 Column 2:

  -0.209700339740800 + 0.503409326295536i
   0.572606616085979 + 0.000000000000000i
   0.143043661588861 - 0.460961450832140i
  -0.375967021752729 + 0.020696164045430i

 Column 3:

   0.648813510859040 + 0.000000000000000i
  -0.271997939078370 + 0.215281719781302i
  -0.240438864070340 + 0.374995359785745i
  -0.053748923044789 - 0.507336787989142i

 Column 4:

  -0.071558944673890 + 0.353911410227957i
  -0.389480689912547 + 0.086669056222871i
  -0.407382573405088 - 0.399879298863578i
   0.620125215785105 + 0.000000000000000i

d =

 Column 1:

   1.806166293907551 + 2.030664098006237i
   0.000000000000000 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i

 Column 2:

   0.000000000000000 + 0.000000000000000i
  -0.136299627666124 - 0.401122045040072i
   0.000000000000000 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i

 Column 3:

   0.000000000000000 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i
   0.169073650200598 + 0.273930811396951i
   0.000000000000000 + 0.000000000000000i

 Column 4:

   0.000000000000000 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i
   0.000000000000000 + 0.000000000000000i
  -0.148497381809724 + 0.329387194182949i


I also don't think that having the fortran or C code to these routines
will (directly) help you get better precision.  Their precision is
mainly due to the precision of the floating point representation of the
machine.  You might take a look at some of the arbitrary precision
arithmetic packages that exist.  I don't know if any of them have
eigenvalues or especially complex eigenvalues, but its worth a shot.

The only one that I can remember off the top of my head is called
"pari".

Good Luck,


Bill Lash
lash at tellabs dot com