Didn't Ric Sherlock reply to this already?
He pointed out that the matrix is not positive definite.
However, this appears not to stymie LAPACK:
load '~addons/lapack/lapack.ijs'
load '~addons/lapack/dgeev.ijs'
coinsert 'jlapack'
A=: 33 16 72 , _24 _10 _57 ,: _8 _4 _17
'leigv eigv reigv'=. dgeev A
leigv
0.62469505 _0.24253563 0
0.62469505 4.3325343e_15 _0.31622777
0.46852129 _0.9701425 0.9486833
eigv
3 1 2
reigv
0.78446454 0.76447079 0.76190476
_0.58834841 _0.61157663 _0.61904762
_0.19611614 _0.20385888 _0.19047619
On 8/11/07, Alfonso Salazar < [EMAIL PROTECTED]> wrote:
>
>
>
> Friends:
>
> Consider in Eugene McDonnell's algorithm:
>
> IR=:@([EMAIL PROTECTED]@])
>
> Does somebody know why for some matriz fail the verb
> below?
> Does somebody know another verb in J in order to
> calculate eigenvalues?
>
> The verb in question is:
>
> NB. ¡Oh No!, ¡Not Eigenvalues Again!
> NB. script of Eugene McDonnell
>
>
> NB. Nota:
> NB. Se tarda mucho con matrices mayores de 9
> renglones
> NB. y muchos de los resultados no concuerdan con
> los
> NB. del programa de las utiler¡as de APL.
>
> NB. En este caso el resultado si coincide con el de
> APL
> A=: 33 16 72 , _24 _10 _57 ,: _8 _4 _17
>
> NB. Sintaxis:
> NB. eigen A
> NB. 3 2 1
>
> eigen=: 3 : 0
> IR=:@([EMAIL PROTECTED]@]) NB. Here is yet
> B1=:<"1
> B0=:<"0
> d=:(<0 1)&|:
> dif=:-/@,:
> difb=:dif&.>
> ppr=:+//.@(*/)
> pprb=:ppr&.>
> detp=:>@(difb/ .pprb)
> db=:B1@(d ,. _1:) d@([EMAIL PROTECTED]@])} B0
> cm=:>@{:@[EMAIL PROTECTED]@db
> (<'eigenvalores con el m ximo'),:<(<<>./r),<"0 r=.cm
> y.
> )
>
> ¡S A L U D!
>
> Alfonso Salazar
>
>
>
>
> ____________________________________________________________________________________
>
> ¡Sé un mejor ambientalista!
> Encuentra consejos para cuidar el lugar donde vivimos.
> http://mx.yahoo.com/promos/mejorambientalista.html
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
--
Devon McCormick, CFA
^me^ at acm.
org is my
preferred e-mail
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
