NB. Rayleigh Quotitent Interation
NB. Cliff Reiter February 2002
NB. This script requires J from http://jsoftware.com
NB. along with the lapack addon from that site.
NB. Also the utility script raster5.ijs should be 
NB. placed into a fvj2 subdirectory of the j directory
NB. It may be obtained from 
NB. http:/www.lafayette.edu/~reiterc/J/fvj2/index.html
NB.
NB. Use J to load this script.
NB. Then the command
NB. m2 200 mk_rq (_15 9j12) 'm2.bmp'
NB. creates the figure like the first one and saves
NB. the result in the file "m2.bmp" in the J directory.

require 'addons\lapack\lapack.ijs addons\lapack\dgeev.ijs'
require 'fvj2\raster5'

NB. comutes unit vectors   
unit=:% [: %: + +/ . * ]

NB. Rayleigh Quotienbt step
rq=: 1 : 0
((+ +/ . * m.&(+/ . *)),])@:unit@:(}.%.m."_ - {.* (=i.#m.)"_)@:}: f. , >:@{:
)

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. iterate until singular or 253 steps are reached
RQ=: rq^:({: < 253"_) :: ]^:_

NB. create an array of input guesses for u0
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. create an image and save a file.
NB. M num_pixels mk_rq clur filename
NB. where clur is the centerleft and upper right
NB. pixels postions given as J complex numbers.

mk_rq=:2 : 0
:
ev=.(/: |.&.|:@:*.)>1{dgeev_jlapack_ x.
z=.m. zl_clur0 n.
'b k'=.0 1|:({.,{:)@(x. RQ)@:(,&(0,~unit 1#~#x.))"0 z
b=.{.@/:@:|@(ev&-)"0 b
(p18p;(253~:k)*1+(3|k)+3*b) writebmp8 y.
)