NB. J5 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.

require 'addons\lapack\lapack.ijs addons\lapack\dgeev.ijs'
require 'addons\image3\view_m.ijs'

NB. matrix product
mp=: +/ . *

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

NB. unit vector
unit =: % elen

NB. M rqik u0
rqik=: 4 : 0"_ 0
k=.0
x=.x0_rqi
u=.y.
while. k<18 do.
  try. x=.unit x %. x. - u *I_rqi catch. break. end.
  u=. (+x) mp x. 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