%***********************************************************************
%=======================================================================
%
% Code for the extension of the Matrix Package to include Sparse
% Matrices.
%
% The following code is for the functions to calculate eigenvalues and
% the rank of sparse matrices.
%
% Author: Stephen Scowcroft. Date: June 1995
%
%=======================================================================
%***********************************************************************
module spmateig;
% Based on the current eigenvalue code.
symbolic procedure spmateigen(u,eival);
% U is a matrix form, eival an indeterminate naming the eigenvalues.
% Result is a list of lists:
% {{eival-eq1,multiplicity1,eigenvector1},....},
% where eival-eq is a polynomial and eigenvector is a matrix.
% How much should we attempt to solve the eigenvalue eq.? sqfr?
% Sqfr is necessary if we want to have the full eigenspace. If there
% are multiple roots another pass through eigenvector calculation
% is needed(done).
% We should actually perform the calculations in the extension
% field generated by the eigenvalue equation(done inside).
begin scalar arbvars,exu,sgn,q,r,s,x,y,z,eivec,!*factor,!*sqfree,
!*exp,val,res,res1,rl;
integer l,cnt;
!*exp := t;
if not(getrtype u eq 'sparse) then typerr(u,"sparse matrix");
eival := !*a2k eival;
kord!* := eival . kord!*;
rl:=sprow_dim(u);
exu := spmateigen1(spmatsm u,eival);
q := car exu;
y := cadr exu;
z := caddr exu;
exu := cdddr exu;
!*sqfree := t;
for each j in cdr fctrf numr subs2(lc z ./ 1)
do if null domainp car j and mvar car j eq eival
then s := (if null red car j
then !*k2f mvar car j . (ldeg car j*cdr j)
else j) . s;
for each j in q
do (if x then rplacd(x,cdr x + cdr j)
else s := (y . cdr j) . s)
where x := assoc(y,s) where y := absf reorder car j;
l := length s;
r := 'list .
for each j in s collect
<<if null((cdr j = 1) and (l = 1)) then
<<y := 1;
for each k in exu do
if x := reduce!-mod!-eig(car j,c!:extmult(k,y))
then y := x>>;
arbvars := nil;
for each k in lpow z do
if (y=1) or null(k member lpow y) then
arbvars := (k . makearbcomplex()) . arbvars;
sgn := (y=1) or evenp length lpow y;
cnt:=0;
for each k in lpow z do
<< if x := assoc(k,arbvars) then
res:=list(cnt:=cnt+1,(1 . mvar cdr x))
else
<< val:= mkgleig(k,y,sgn := not sgn,arbvars);
if (val=simp 0) then cnt:=cnt+1
else res:=list(cnt:=cnt+1,(1 . prepsq!* val));
>>;
if res=nil then nil
else << res1:=append(res1,list(res));
res:=nil;
eivec:=mkempspmat(rl,list('spm,rl,1));
for i:=1:rl do
<<letmtr3(list(eivec,i),atsoc(i,res1),eivec,nil);>>;
>>;
>>;
res1:=nil;
list('list,prepsq!*(car j ./ 1),cdr j,eivec)>>;
kord!* := cdr kord!*;
return r
end;
symbolic procedure spmateigen1(u,eival);
begin scalar diag,q,x,y,z,w,j,res; integer l,lm,m,cc;
lm := spcol_dim(u);
z := 1;
for rr:=1:sprow_dim(u) do
<< y := 1;
diag:=nil;
cc:=findrow(u,rr);
if not (cc=nil) then
<<for each xx in cdr cc do
<<w:=simp cdr xx;
y := lcm(y,denr w);
>>;
>>;
m := lm;
l := l + 1;
x := nil;
if cc=nil then cc:=list(list(nil),list(nil));
for each xx in reverse cdr cc do
<<if xx='(nil) then <<m:=rr+1; j:=simp nil>>
else << j:=simp cdr xx; m:=car xx>>;
if not diag then
<<if m=rr then <<diag:=j; j:=simp nil>>
else if m<rr then nil
else <<diag:=simp findelem2(u,rr,rr); m:=rr>>;
x:=list m .* multf(addf(numr diag,negf multf(!*k2f eival,
denr diag)),quotf(y,denr diag)) .+ x;
if not (xx=nil) then m:=car xx;
>>;
if numr j and not (m=rr) then
x := list m .* multf(numr j,quotf(y,denr j)) .+ x;
>>;
y := z;
z := c!:extmult(if null red x then <<
q := (if p then (car p . (cdr p + 1)) . delete(p,q)
else (lc x . 1) . q) where p = assoc(lc x,q);
!*p2f lpow x>> else x,z);
res:=append(res,list x);
>>;
u:=res;
return q . y . z . u
end;
flag('(spmateigen),'opfn);
flag('(spmateigen),'noval);
% These functions calculate the rank of a sparse matrix.
% To replace current function.
% Extended for sparse matrices.
symbolic procedure rank!-eval u;
begin scalar n;
if cdr u then rerror(matrix,17,"Wrong number of arguments")
else if getrtype (u := car u) eq 'matrix
then return rank!-matrix spmatsm u
else if getrtype (u) eq 'sparse
then return sprank!-matrix spmatsm u
else if null eqcar(u := aeval u,'list) then typerr(u,"matrix")
else return rank!-matrix
for each row in cdr u collect
if not eqcar(row,'list)
then rerror(matrix,15,"list not in matrix shape")
else <<row := cdr row;
if null n then n := length row
else if n neq length row
then rerror(matrix,151,"list not in matrix shape");
for each j in row collect simp j>>
end;
put('rank,'psopfn,'rank!-eval);
% The function to calculate the rank.
% Based on the matrix rank function.
symbolic procedure sprank!-matrix u;
begin scalar x,y,z,w,j,cc; integer m,n;
z := 1;
for rr:=1:sprow_dim(u) do
<< y := 1;
cc:=findrow(u,rr);
if not (cc=nil) then
<<for each xx in cdr cc do
<<w:=simp cdr xx;
y := lcm(y,denr w);
>>;
>>;
m := 1;
x := nil;
if not (cc=nil) then
<<for each xx in cdr cc do
<<j:=simp cdr xx;
m:=car xx;
if numr j then
x := list m .* multf(numr j,quotf(y,denr j)) .+ x;
>>;
>>;
if y := c!:extmult(x,z) then <<z := y; n := n + 1>>;
>>;
return n
end;
endmodule;
end;
%***********************************************************************
%=======================================================================
%
% End of Code.
%
%=======================================================================
%***********************************************************************