;  produce an animated graph showing the progress of the
;  power method for calculation of eigen values
 
(def dim 20)
(def mat (read-data-file "dat.20"))
(def y (repeat 0 dim))
; the matrix mat, read in by rows, must be converted to a matrix
; with the make-array function before calling eigen
; eigen returns a list containing the eigenvalues and vectors
; the select function pulls off the first element, the eigenvectors
(def ee (select (eigen (make-array (list dim dim) :initial-contents (list mat)))
 0))
 
; dynamically load the C routine  and call the powr function
; to create the executable use
;    cc -c -O -KPIC pwr.c
;    cc -c -O -KPIC gs.c
;    cc --f77 -G -o pwr.so pwr.o gs.o -R/usr/local/lib -lnag
(dyn-load "./pwr.so")
; call-cfun returns a list containing its (possibly updated)
; arguments - we want the second one, i.e. (select ... 1)
(def res (call-cfun "powr" (float mat) (float y) dim ))
 
(def e1 (select res 1))
; produce an initial plot of sequence number versus true eigenvalue
(def pl (plot-lines (iseq 1 dim) ee))
 
(def e1 (+ 1 ee))
; keep reevaluating the power method until the maximum difference
; between the estimate and the true value is 1e-3
(loop   (if (< (max (abs (- e1 ee))) 1.e-3) (return))
        (def res (call-cfun "powr" (float mat) (float y) dim ))
(def e1 (select res 1))
;  pause so that the points don't move too fast
        (pause 5)
;  clear the points and add the new estimates to the graph
;  note that clear-points does not affect the original line which
;  was drawn with plot-lines
;  the call to (gc) forces garbage collection, which seems to be
;  necessary in iterative programs
        (send pl :clear-points)
        (send pl :add-points (iseq 1 dim) e1)
        (gc))