From: ljdickey@watmath.waterloo.edu (L. J. Dickey)
Subject: Eigen in J now works correctly
Keywords: eigenstructure, J
Organization: University of Waterloo
Date: Sun, 5 Apr 1992 06:12:16 GMT


About a month ago Brian wrote

	Subject: Eigen in J?
	Date: Mon, 2 Mar 1992 13:33:56 GMT
	
	I am currently looking for a J function to compute the
	eigenstructure of a symmetric matrix.
	...

and Raul came back with
	Sure, why not.  It's been a while since I've studied linear algebra,
	so I don't have a real good feel for this kind of math, but bear with
	me  ...

Then Brian corrected something by saying:
	The leftmost "[" should be "]".

Thank you both.

Yesterday, I did some experiments with the code, and found that another
change was needed.  This version is now posted completely below, in two
chunks.  The first chunk is the complete code that does the job,
including one small example together with what I believe to be a
correct answer.  If you are inclined to do so, you may run this example
and compare with the results I have given.  Secondly I have a short
script that can be run repeatedly to generate random symmetric
matrices, and with these, put the code through its paces.  This section
has a line or two that checks the accuracy of the results.  I have called
this second chunk "t" and keep it in a file by that name.  Then with the
initialized (via profile) command scr, given by

	scr	=. 0!:2 & <

I can repeatedly type scr 't' to execute the code in file 't', getting
different files "aa" and not surprisingly, different execution times.  
Maybe there is a faster algorighm?

						Lee Dickey

 === begin file eigen =========================================================
 sym =. '0.5*y.+|:y.' : ''		NB. symmetrize a matrix
 
 eigen =. ] (] , (eigval ])) eigvect	NB. corrected
 dot	=. +/ .*
 eigval	=. (<0 1)&|: @ (|: @] dot [ dot ])	NB. corrected

 t=.0 0$''					NB. literal translation
 t=.t,  'y=. ident  n=. # {. a=. y.'		NB. try and use same names
 t=.t,  ' eps=. n * n * 1e_12 * +/ | (<0 1)|: a'
 t=.t,  'qr)  q=. ident  i=. n'
 t=.t,  ' r=. a'
 t=.t,  'again)  z=. ident n'
 t=.t,  ' z=. z (n #.&> {;~ (-i){. i. n) }~ house i $ r {~ {((n-i)+i. i); n-i'
 t=.t,  ' q=. q dot z'
 t=.t,  ' r=. z dot r'
 t=.t,  '$.=. $. ]^:(1< i=.i-1) again'
 t=.t,  ' a=. r dot q'
 t=.t,  ' y=. y dot q'
 t=.t,  '$.=. $. ]^:(eps < (+/,a^2)-+/((<0 1)|: a)^2) qr'
 t=.t,  ' y=. y {"1~ \: (<0 1)|: a'
 eigvect=. t : ''
 
 ident =. = @i.
 house=. # (ident@[ - 2&*@] */ ] % +/@(^&2)@]) (] - # {. +/&.(^&2))
 
 eigen 2 2 $ 1 2 2 3

NB. 0.525731 _0.850651
NB. 0.850651  0.525731
NB.  4.23607 _0.236068
 === end of file e ============================================================

 === file t ===================================================================
   aa =. sym ? 3 3 $ 100
   bb =. eigvect aa
   bb %~ aa dot bb 		NB. This sometimes gives e-vals
   (<0 1) |: (|:bb) dot aa dot bb		NB.  This gives e-vals
   aa eigval bb
   >. / ,|(aa dot bb) - (aa eigval bb)*bb	NB. accuracy 
 === end of file t ============================================================
-- 
Leroy J. Dickey, Faculty of Mathematics, U of Waterloo, Canada  N2L 3G1
Internet:      ljdickey@watmath.UWaterloo.ca
               ljdickey@watmath.waterloo.edu
UUCP:          ljdickey@watmath.UUCP