Raul, I have asked Igor Zhuravlov principle author of mt whether he
knows of documentation. Below is what I use. --Kip Murray
NB. This script defines four linear algebra utilities
NB. rref -- row reduced echelon form of matrix y
NB. null -- basis for the null space of matrix y
NB. eit -- table of eigenvalues and eigenvectors of square matrix y
NB. eig -- eigenvalues of square matrix y
NB. eiv -- independent eigenvectors of y belonging to eigenvalue x
NB. These routines behave well with small well-behaved problems.
NB. They are not a substitute for LAPACK!
R =: ]
L =: [: : [
R0 =: (0 {:: ]) :: ([: > ])
R1 =: 1 {:: ]
R2 =: 2 {:: ]
R3 =: 3 {:: ]
col =: ,. NB. arranges vectors and matrices into columns of a matrix
row =: ,. &. |: NB. arranges vectors and matrices into rows
the =: [:
clean =: * *@| NB. numbers with 2^_44 > | are replaced by 0
rref =: finish@(repeat^:while^:_)@start
start =: the (R0@$ ; R1@$ ; R) clean
while =: (0 < R0) *. 0 < R1
repeat =: the pivot R0 ; R1 ; R0 sort R2
sort =: (\: |)@{. , }.
pivot =: if0`if1@.(0: ~: (<0 0)"_ { R2)
if0 =: R0 ; (the <: R1) ; the shl R2
shl =: 1 |."1 R
if1 =: (the <: R0) ; (the <: R1) ; the clean the sul the pvt the prp R2
prp =: (the }. R) ,~ the (% R0) the {. R
pvt =: (the {. R) (L , R - R0"1 */ L) the }. R
sul =: the shl 1 |. R
finish =: R1 |."1 R0 |. R2
NB. columns of null y are a basis for the null space of matrix y
n =: the {: the $ ]
ones =: (] i."1 1:) -. n
free =: i.@n -. ones
nullfrom =: free (the - {"1)`ones`((the +./ [ =/ the i. n) expand the =
[) } ]
null =: the nullfrom the (-. n # 0:) rref
NB. eit "eigentable"
o =: +/ . * NB. Matrix multiply as in A o B
det =: (-/ . *)@]
trace =: +/@((<0 1)&|:)@]
eye =: =@i.@] NB. identity matrix
frame =: 3 : 0 NB. char pol coefficients
F =. I =. =i.m =. #M =. y
C =. 0
for_k. 1+i.m do.
F =. M +/ . * F + ({.C)*I
C =. (-(trace F)%k),C
end.
(}:C),1
)
NB. eig =: R1@p.@frame@] NB. eigenvalues of n by n matrix y
NB. improve repeated roots from p.
improve =: R0 ; the ; the ((the #. (3.0517578125e_5"_ > the | 2: * - %
+&|)"0/~) (the < # # +/ % #)/. ]) R1
eig =: R1@improve@p.@frame@]
eiv =: the null ] - [ * eye@#@] NB. eigenvectors of y given eigenvalue x
eit =: (the <&> the ~. eig) ([ row eiv&.>) <@] NB. "eigentable" for y
On 10/5/2013 8:13 AM, Raul Miller wrote:
This is definitely a fun exercise, and I'm having fun playing with it
(mostly looking at intermediate results and a few variations on the
theme).
That said, what I was really looking for was documentation or
explorations on J's facilities (such as mt or lapack) for working with
eigenvalues (and eigenvectors - in principle if you have the one the
other is straightforward).
Here, though, you did not use mt, nor did you use lapack. You started
with a specific set of eigenvalues and demonstrated a context where
you could plug them in.
And I certainly appreciate the effort - and it does demonstrate some
of the identities involving their use - but of course I am looking for
more!
Thanks,
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
