On Mon, 11 Dec 2006, sen1 at math.msu.edu wrote:
whoops, in the last post I had some other aliases (e.g. "flt" for "float",
"prt" for "print".
Here is a new attachment (I hope I got rid of them all).
-sen
> Here is the program as an attachment.
>
> The file (which is a text file) is called "Esparse_test.mac"
>
> It loads the program called "Esparse" which is described in the
> beginning comments.
>
> The offending 27x27 matrix "A" is also loaded so one can try "Esparse"
> on that matrix.
>
> Thanks for any comments.
>
> -sen
>
-------------- next part --------------
/****************** File Esparse_test.mac *****************/
/* finds the largest eigenvalue of a sparse matrix A
Typically, the matrix is entered as
A: zeromatrix(SIZE,SIZE);
Then, the entries can be modified by specifying the A[i,j] terms.
The variables MAX and ERR are the number of iterations, and error
desired.
Example: A: zeromatrix(4,4);
A[1,1]:1.7;
A[1,2]:2.8;
A[3,2]:4;
Esparse(A,4,30,1e-10);
*/
ratprint: false;
load("eigen")$
Norm(list):= flt(sqrt(inprod(list,list)))$
Esparse(A,SIZE,MAX,ERR):=
block( [ran,nran, v_old, iter:0, nv_new, v_new, expan, s:SIZE,B],
/* initialize the normalized random s-vector */
ran: makelist(random(1.0),i,1,s),
nran: ran/Norm(ran),
v_old: nran,
while iter < MAX do
( v_new: makelist( sum(A[i,j]*v_old[j],j,1,s),i,1,s),
expan: Norm(v_new),
nv_new: v_new/expan,
if Norm(nv_new - v_old) < ERR then
( v_new: makelist( sum(A[i,j]*v_old[j],j,1,s),i,1,s),
expan: Norm(v_new), return(expan))
else
(v_old: nv_new, iter: iter+1,
if iter=MAX -1 then
prt("MAX reached and ERR not achieved; try decreasing
ERR or increasing MAX ")
)
)
)$
/********** End of Esparse_test.mac ******************/
/******** 27x27 matrix which causes problems with default cmucl *****/
A: matrix([0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0],
[2,0,2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0],
[2,0,2,0,2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0])$