NB. J utilities to creat Rayleigh Quotient iteration images
NB. Cliff Reiter, December 21, 2003
NB. Requires the use of the Lapack and Image3 addons
NB. run script; then sample usage is:
NB.
NB. m2 200 mk_rqi_pic _15 9j12                      NB. create an image
NB. m7 200 mk_rqi_pic _8 8j8
NB. (m7 200 mk_rqi _8 8j8) write_image 'm7.png'     NB. save an image
NB.
NB. updated for J601 Nov 8, 2007

require '~addons\math\lapack\lapack.ijs'
require '~addons\math\lapack\dgeev.ijs'
require '~addons\media\image3\view_m.ijs'

NB. matrix product
mp=: +/ . *

NB. Euclidean length
elen=: (+/&.:*:)@:|

NB. unit vector
unit=: % elen

NB. A rqik u0
rqik=: 4 : 0"_ 0
A=. x
k=. 0
x=. x0_rqi
u=. y
while. k<18 do.
  try. x=. unit x %. A - u *I_rqi catch. break. end.
  u=. (+x) mp A mp x
  k=. k+1
end.
({./:|EV_rqi-u),k
)

NB. Makes an image
mk_rqi=: 1 : 0
:
EV_rqi=: (/: |.&.|:@:*.)>1{dgeev_jlapack_ x
I_rqi=: =i.#x
x0_rqi=: (#x)#1
z=. m zl_clur0 y
'b k'=. 0 1|:x rqik z
p18p;(18~:k)*1+(3|k)+3*b
)

NB. show the image
mk_rqi_pic=: 1 : 0
:
view_image x m mk_rqi y
)

NB. get a mesh in the complex plane
zl_clur0=: 4 : 0
w=. -~/9 o.y
h=. -~/11 o.y
(j.h-(+:h)*(i.%<:)x) +/ ({.y)+w*(i.%<:)x
)

NB. a palette
$p18p=: 0,,/0 2|:<.1 0.72 0.45 */ Hue 0 2 4 1 3 5{(i.%])6

NB. the matrices from [1]
m1=: 1|.=i.5
m2=: 1,9 8 4,:5 6 2
m3=: 1,9j1 8 4,:5 6 2
m4=: ".;.(_2) 0 : 0
_6 7 1 _5
7 10 3 _4
1 3 8 _7
_5 _4 _7 8
)
m5=: ".;.(_2) 0 : 0
2 _9j_6 _3j_1 _3j1
_9j6 _4 _9j2 _6j6
_3j1 _9j_2 2 0j_3
_3j_1 _6j_6 0j3 2
)
m6=: ".;.(_2) 0 : 0
1 2 3 4
2 5 6 7
3 6 8 9
4 7 9 10
)
m7=: ".;.(_2) 0 : 0
_14 _10 0 0
21 12 _4 0
_15 _9 3  1
10 7 _1 _1
)
m8=: ".;.(_2) 0 : 0
12 _14 0 4
7 _6 _2  4
_10 12 _2 _2
_18 22 _3 _4
)

NB. [1] C. Reiter, Visualizing the Dynamics of the Rayleigh Quotient Iteration, Computers & Graphics, 16 3 (1992) 341-344