; 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))