%***********************************************************************
%=======================================================================
%
% Code for the Linear Algebra for Sparse Matrices Package.
%
% Author: Stephen Scowcroft. Date: June 1995
%
%=======================================================================
%***********************************************************************
module splinalg;
% This is the beginning of a Linear Algebra Package for Sparse Matrices.
% It stems from the present Linear Algebra Package for Matrices and
% will hopefully be incorporated into it to produce a new MATRIX
% package.
switch fast_la; % If ON, then the following functions will be faster:
% spadd_columns spadd_rows spaugment_columns spcolumn_dim %
% spcopy_into spmake_identity spmatrix_augment spmatrix_stack %
% spminor spmult_column spmult_row sppivot %
% spremove_columns spremove_rows sprows_pivot spsquarep %
% spstack_rows spsub_matrix spswap_columns spswap_entries %
% spswap_rows spsymmetricp %
% This is basically done by removing some error checking and doesn't
% speed things up too much. You'll need to be making alot of calls to
% see the difference. If you get strange error messages with fast_la
% ON then thoroughly check your input.
% Creates an identity matrix of dimension size.
symbolic procedure spmake_identity(size);
begin scalar tm;
if not !*fast_la and not fixp size then
rederr "Error in spmake_identity: non integer input.";
tm:=list('sparsemat,mkvect(size),list('spm,size,size));
for i :=size step -1 until 1 do
<<letmtr3(list(tm,i,i),1,tm,nil)
>>;
return tm;
end;
flag('(spmake_identity),'opfn);
% Finds the row-wise dimension of a matrix.
symbolic procedure sprow_dim(u);
begin scalar res;
if not !*fast_la and not matrixp(u) then
rederr "Error in sprow_dim: input should be a matrix.";
res := cadr spmatlength(u);
return res;
end;
% Finds the column-wise dimension of a matrix.
symbolic procedure spcol_dim(u);
begin scalar res;
if not !*fast_la and not matrixp(u) then
rederr "Error in spcol_dim: input should be a matrix.";
res := caddr spmatlength(u);
return res;
end;
flag('(sprow_dim,spcol_dim),'opfn);
% Test if input is a matrix or sparse matrix (boolean)
% Whether id is a gerneric matrix.
symbolic procedure matrixp(u);
begin;
if not pairp u then u:=reval u;
if not eqcar(u,'mat) and not eqcar(u,'sparsemat) then return nil
else return t;
end;
flag('(matrixp),'boolean);
flag('(matrixp),'opfn);
% Test to see if input is a sparse matrix or not. (boolean)
symbolic procedure sparsematp(u);
if not eqcar(u,'sparsemat) then nil
else t;
flag('(sparsematp),'boolean);
flag('(sparsematp),'opfn);
% Tests matrix is square. (boolean).
symbolic procedure squarep(u);
begin
scalar tmp;
if not !*fast_la and not matrixp(u) then
rederr "Error in squarep: non matrix input";
if sparsematp(u) then tmp := cdr spmatlength(u)
else tmp:=size_of_matrix(u);
if car tmp neq cadr tmp
then return nil
else return t;
end;
flag('(squarep),'boolean);
flag('(squarep),'opfn);
% Checks input is symmetric. ie: transpose(A) = A. (boolean).
symbolic procedure symmetricp(u);
if eqcar(u,'sparsemat) then
if smtp (u,nil) neq u then nil else t
else
if algebraic (tp(u)) neq u then nil else t;
flag('(symmetricp),'boolean);
flag('(symmetricp),'opfn);
% Takes a constant (const) and an integer (mat_dim) and creates
% a jordan block of dimension mat_dim x mat_dim.
symbolic procedure spjordan_block(const,mat_dim);
begin scalar tm;
tm :=mkempspmat(mat_dim,list('spm,mat_dim,mat_dim));
if not fixp mat_dim then rederr
"Error in spjordan_block(second argument): should be an integer.";
for i:=mat_dim step -1 until 1 do
<<letmtr3(list(tm,i,i),const,tm,nil);
if i<mat_dim then letmtr3(list(tm,i,i+1),1,tm,nil);
>>;
return tm;
end;
flag ('(spjordan_block),'opfn);
% Removes row (row) and column (col) from in_mat.
symbolic procedure spminor(list,row,col);
begin scalar len,lena,lenb,rlist;
len:=caddr list;
rlist :=copy_vect(list,nil);
lena := cadr len;
lenb := caddr len;
if not matrixp(list) then
rederr "Error in spminor(first argument): should be a matrix.";
if not fixp row then
rederr "Error in spminor(second argument): should be an integer.";
if not fixp col then
rederr "Error in spminor(third argument): should be an integer.";
if not (row > 0 and row < lena + 1)
then rerror(matrix,20,"Row number out of range");
if not (col > 0 and col < lenb + 1)
then rerror(matrix,21,"Column number out of range");
spremrow(row,rlist);
spremcol(col,rlist);
rlist := rewrite(rlist,lena-1,row,col);
return rlist;
end;
flag('(spminor),'opfn);
% A square band matrix b is created. The elements of the diagonal
% are the middle element of elt_list. The elements to the left are
% used to fill the required number of subdiagonals and the elements
% to the right the superdiagonals.
symbolic procedure spband_matrix(elt_list,sq_size);
begin scalar tm;
integer i,j,it,no_elts,middle_pos;
tm:=mkempspmat(sq_size,list('spm,sq_size,sq_size));
if not fixp sq_size then rederr
"Error in spband_matrix(second argument): should be an integer.";
if atom elt_list then elt_list := {elt_list}
else if car elt_list = 'list then elt_list := cdr elt_list
else rederr
"Error in spband_matrix(first argument): should be single value or list.";
no_elts := length elt_list;
if evenp no_elts then rederr
"Error in spband matrix(first argument): number of elements must be odd.";
middle_pos := reval{'quotient,no_elts+1,2};
if my_reval middle_pos > sq_size then rederr
"Error in spband_matrix: too many elements. Band matrix is overflowing.";
it:=2;
for i:=1:sq_size do
<< if i<=middle_pos then j:=1
else j:=it;
while middle_pos-i+j > 0 and
j<=sq_size and middle_pos-i+j <= no_elts do
<<
letmtr3(list(tm,i,j),nth(elt_list,middle_pos-i+j),tm,nil);
j:=j+1;
>>;
if i>middle_pos then it:=it+1;
>>;
return tm;
end;
flag('(spband_matrix),'opfn);
% Stacks all rows pointed to in row_list to form a new matrix.
%
% row_list can be either an integer or a list of integers.
symbolic procedure spstack_rows(in_mat,row_list);
begin scalar tl,rlist,res,list;
integer rowdim,cnt,coldim;
list:=in_mat;
cnt:=1;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in spstack_rows(first argument): should be a matrix.";
if atom row_list then row_list := {row_list}
else if car row_list = 'list then row_list := cdr row_list
else
<<
prin2 "***** Error in spstack_rows(second argument): ";
prin2t " should be either an integer or a list of integers.";
return;
>>;
coldim := spcol_dim(in_mat);
rowdim := sprow_dim(in_mat);
tl:=list('smp,length(row_list), coldim);
res:=mkempspmat(length(row_list),tl);
for each elt in row_list do
<<
if not fixp elt then rederr
"Error in spstack_rows(second argument): contains non integer.";
if elt>rowdim or elt=0 then
<<
prin2 "***** Error in spstack_rows(second argument): ";
rederr "contains row number which is out of range for input matrix.";
>>;
rlist:= findrow(list,elt);
if rlist = nil then cnt := cnt + 1
else
<< letmtr3(list(res,cnt),rlist,res,nil);
cnt := cnt + 1
>>;
>>;
return res;
end;
% Augments all columns pointed to in col_list to form a new matrix.
%
% col_list can be either an integer or a list of integers.
symbolic procedure spaugment_columns(in_mat,col_list);
begin
integer cnt,coldim,rcnt,rowdim;
scalar tl,rlist,res,list,rrlist,colist,val,res1;
list:=in_mat;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in spaugment_columns(first argument): should be a matrix.";
if atom col_list then col_list := {col_list}
else if car col_list = 'list then col_list := cdr col_list
else
<<
prin2 "***** Error in spaugment_columns(second argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
cnt:=1;
rcnt:=1;
tl := list('spm,rowdim,length(col_list));
res:=mkempspmat(rowdim,tl);
for each elt in col_list do
<< if not fixp elt then rederr
"Error in spaugment_columns(second argument): contains non integer.";
if elt>coldim or elt=0 then
<< prin2 "***** Error in spaugment_columns(second argument): ";
rederr
"contains column number which is out of range for input matrix.";
>>;
>>;
for i:=1:rowdim do
<< rrlist:=findrow(list,i);
if rrlist then
<<for each elt in col_list do
<<colist:=atsoc(elt,rrlist);
if colist = nil then cnt := cnt + 1
else << val := cdr colist;
res1:=(cnt . val);
rlist := res1 . rlist;
cnt := cnt + 1;
>>;
>>;
if rlist then
letmtr3(list(res,i),list(nil) . reverse rlist,res,nil);
rlist := nil;
cnt := 1;
>>;
>>;
return res;
end;
flag('(spstack_rows,spaugment_columns),'opfn);
% Input is a matrix and either a single row number or a list of row
% numbers.
%
% Extracts either a single row or a number of rows and returns them
% in a list of row matrices.
symbolic procedure spget_rows(in_mat,row_list);
begin scalar tl,he,rlist,res,list,rlist1;
integer rowdim,cnt,coldim;
coldim:= spcol_dim(in_mat);
tl:=list('spm,length(row_list), coldim);
list:=in_mat;
cnt:=1;
if not matrixp(in_mat) then
rederr "Error in spget_rows(first argument): should be a matrix.";
if atom row_list then row_list := {row_list}
else if car row_list = 'list then row_list := cdr row_list
else
<<
prin2 "***** Error in spget_rows(second argument): ";
prin2t " should be either an integer or a list of integers.";
return;
>>;
rowdim := sprow_dim(in_mat);
for each elt in row_list do
<<
if not fixp elt then rederr
"Error in spget_rows(second argument): contains non integer.";
if elt>rowdim or elt=0 then
<<
prin2 "***** Error in spget_rows(second argument): ";
rederr "contains row number which is out of range for input matrix.";
>>;
rlist:= findrow(list,elt);
if rlist = nil then nil
else
<< rlist1:= mkempspmat(1,list('spm,1,coldim));
letmtr3(list(rlist1,1),rlist,rlist1,nil);
res:=append(res,{rlist1});
>>;
>>;
return 'list.res;
end;
% Input is a matrix and either a single column number or a list of
% column numbers.
%
% Extracts either a single column or a series of adjacent columns and
% returns them in a list of column matrices.
symbolic procedure spget_columns(in_mat,col_list);
begin
integer coldim,rcnt,rowdim;
scalar tl,rlist,res,list,nlist,rrlist,colist,val,res1;
rowdim:= sprow_dim(in_mat);
tl := list('spm,rowdim,length col_list);
list:=in_mat;
if not matrixp(in_mat) then
rederr "Error in spget_columns(first argument): should be a matrix.";
if atom col_list then col_list := {col_list}
else if car col_list = 'list then col_list := cdr col_list
else
<<
prin2 "***** Error in spget_columns(second argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
coldim := spcol_dim(in_mat);
rcnt:=1;
for each elt in col_list do
<< if not fixp elt then rederr
"Error in spget_columns(second argument): contains non integer.";
if elt>coldim or elt=0 then
<< prin2 "***** Error in get_columns(second argument): ";
rederr
"contains column number which is out of range for input matrix.";
>>;
>>;
for each elt in col_list do
<<rlist:=mkempspmat(coldim,list('spm,coldim,1));
for i:=1:rowdim do
<< rrlist:=findrow(list,i);
if rrlist then
<<colist:=atsoc(elt,rrlist);
if colist = nil then rcnt := rcnt + 1
else<< val := cdr colist;
res1 := {1 . val};
letmtr3(list(rlist,rcnt),
list(nil) . res1,rlist,nil);
rcnt := rcnt + 1;
>>;
>>;
>>;
res := append(res,{rlist});
rlist := nil;
rcnt := 1;
>>;
return 'list.res;
end;
flag('(spget_rows,spget_columns),'opfn);
% Removes each row in row_list from in_mat.
%
% row_list can be either an integer or a list of integers.
symbolic procedure spremove_rows(in_mat,row_list);
begin
scalar unique_row_list,list,tl;
integer rowdim,row,cnt;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in spremove_rows(first argument): non matrix input.";
if atom row_list then row_list := {row_list}
else if car row_list = 'list then row_list := cdr row_list
else
<<
prin2 "***** Error in spremove_rows(second argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
% Remove any repititions in row_list (I'm assuming here that if the
% user has inputted the same row more than once then the meaning
% is to only remove that row once).
unique_row_list := {};
for each row in row_list do
<<
if not intersection({row},unique_row_list) then
unique_row_list := append(unique_row_list,{row});
>>;
rowdim := sprow_dim(in_mat);
if not !*fast_la then
<<
for each row in unique_row_list do if not fixp row then rederr
"Error in spremove_rows(second argument): contains a non integer.";
for each row in unique_row_list do if row>rowdim or row=0 then
rederr
"Error in spremove_rows(second argument): out of range for input matrix.";
if length unique_row_list = rowdim then
<<
prin2 "***** Warning in spremove_rows:";
prin2t " all the rows have been removed. Returning nil.";
return nil;
>>;
>>;
cnt:=0;
tl:=list('spm,rowdim - length unique_row_list,spcol_dim(in_mat));
list:=copy_vect(in_mat,tl);
for each elt in unique_row_list do
<< spremrow(elt-cnt,list);
list := rewrite(list,rowdim,elt-cnt,0);
cnt := cnt + 1
>>;
return list;
end;
% Removes each column in col_list from in_mat.
%
% col_list can be either an integer or a list of integers.
symbolic procedure spremove_columns(in_mat,col_list);
begin
scalar unique_col_list,tl,list;
integer coldim,col,cnt;
if not !*fast_la and not matrixp(in_mat) then rederr
"Error in spremove_columns(first argument): non matrix input.";
if atom col_list then col_list := {col_list}
else if car col_list = 'list then col_list := cdr col_list
else
<<
prin2 "***** Error in spremove_columns(second argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
% Remove any repititions in col_list (I'm assuming here that if the
% user has inputted the same column more than once then the meaning
% is to only remove that column once).
unique_col_list := {};
for each col in col_list do
<<
if not intersection({col},unique_col_list) then
unique_col_list := append(unique_col_list,{col});
>>;
coldim := spcol_dim(in_mat);
if not !*fast_la then
<<
for each col in unique_col_list do if not fixp col then rederr
"Error in spremove_columns(second argument): contains a non integer.";
for each col in unique_col_list do if col>coldim or col=0 then
rederr
"Error in spremove_columns(second argument): out of range for matrix.";
if length unique_col_list = coldim then
<<
prin2 "***** Warning in spremove_columns: ";
prin2t " all the columns have been removed. Returning nil.";
return nil;
>>;
>>;
cnt:=0;
tl:=list('spm,sprow_dim(in_mat), coldim - length unique_col_list);
list := copy_vect(in_mat,tl);
for each elt in unique_col_list do
<< spremcol(elt-cnt,list);
list := rewrite(list,coldim,0,elt-cnt);
cnt := cnt + 1;
>>;
return list;
end;
flag('(spremove_rows,spremove_columns),'opfn);
% Creates a matrix consisting of rows*cols matrices which are taken
% sequentially from the mat_list.
symbolic procedure spblock_matrix(rows,cols,mat_list);
begin
scalar block_mat,row_list;
integer rowdim,coldim,start_row,start_col,i,j;
if not fixp rows then rederr
"Error in block_matrix(first argument): should be an integer.";
if rows=0 then
<<
prin2 "***** Error in spblock_matrix(first argument): ";
prin2t " should be an integer greater than 0.";
return;
>>;
if not fixp cols then rederr
"Error in spblock_matrix(second argument): should be an integer.";
if cols=0 then
<<
prin2 "***** Error in spblock_matrix(second argument): ";
prin2t " should be an integer greater than 0.";
return;
>>;
if matrixp mat_list then mat_list := {mat_list}
else if pairp mat_list and car mat_list = 'list then
mat_list := cdr mat_list
else
<<
prin2 "***** Error in spblock_matrix(third argument): ";
prin2t
" should be either a single matrix or a list of matrices.";
return;
>>;
if rows*cols neq length mat_list then rederr
"Error in spblock_matrix(third argument): Incorrect number of matrices.";
row_list := spcreate_row_list(rows,cols,mat_list);
rowdim := spcheck_rows(row_list);
coldim := spcheck_cols(row_list);
block_mat := mkempspmat(rowdim,list('spm,rowdim,coldim));
start_row := 1; start_col := 1;
for i:=1:length row_list do
<<
for j:=cols step -1 until 1 do
<<
block_mat := spcopy_into(nth(nth(row_list,i),j),block_mat,
start_row,start_col);
start_col := start_col + spcol_dim(nth(nth(row_list,i),j));
>>;
start_col := 1;
start_row := start_row + sprow_dim(nth(nth(row_list,i),1));
>>;
return block_mat;
end;
flag('(spblock_matrix),'opfn);
symbolic procedure spcreate_row_list(rows,cols,mat_list);
%
% Takes mat_list and creates a list of rows elements each of which is
% a list containing cols elements (ordering left to right).
% eg: create_row_list(3,2,{a,b,c,d,e,f}) will return
% {{a,b},{c,d},{e,f}}.
%
begin
scalar row_list,tmp_list,list;
integer i,j,increment;
increment := 1;
for i:=1:rows do
<<
tmp_list := {};
for j:=1:cols do
<<list:=nth(mat_list,increment);
if not sparsematp(list) then list:=sptransmat(list);
tmp_list := append(tmp_list,{list});
increment := increment + 1;
>>;
row_list := append(row_list,{tmp_list});
>>;
return row_list;
end;
symbolic procedure spcheck_rows(row_list);
%
% Checks all matrices in each element in row_list contains same
% amount of rows.
% Returns the sum of all of these row numbers (ie: number of rows
% required in the block matrix).
%
begin
integer i,listlen,rowdim,eltlen,j;
i := 1;
listlen := length(row_list);
while i<=listlen do
<<
eltlen := length nth(row_list,i);
j := 1;
while j<eltlen do
<<
if sprow_dim(nth(nth(row_list,i),j)) =
sprow_dim(nth(nth(row_list,i),j+1)) then j := j+1
else
<<
prin2 "***** Error in spblock_matrix: row dimensions of ";
rederr "matrices into spblock_matrix are not compatible";
>>;
>>;
rowdim := rowdim + sprow_dim(nth(nth(row_list,i),j));
i := i+1;
>>;
return rowdim;
end;
symbolic procedure spcheck_cols(row_list);
%
% Checks each element in row_list has same number of columns.
% Returns this number.
%
begin
integer i,listlen;
i := 1;
listlen := length(row_list);
while i<listlen
and spno_cols(nth(row_list,i)) = spno_cols(nth(row_list,i+1))
do i:=i+1;
if i=listlen then return spno_cols(nth(row_list,i))
else
<<
prin2
"***** Error in spblock_matrix: column dimensions of matrices ";
prin2t " into spblock_matrix are not compatible";
return;
>>;
end;
symbolic procedure spno_rows(mat_list);
%
% Takes list of matrices and sums the no. of rows.
%
for each mat1 in mat_list sum sprow_dim(mat1);
symbolic procedure spno_cols(mat_list);
%
% Takes list of matrices and sums the no. of columns.
%
for each mat1 in mat_list sum spcol_dim(mat1);
% Copies matrix BB into AA with BB(1,1) at AA(p,q).
symbolic procedure spcopy_into(BB,AA,p,q);
begin
scalar A,B;
integer m,n,r,c,val,j,col;
if not !*fast_la then
<<
if not matrixp(BB) then rederr
"Error in spcopy_into(first argument): should be a matrix.";
if not matrixp(AA) then rederr
"Error in spcopy_into(second argument): should be a matrix.";
if not fixp p then rederr
"Error in spcopy_into(third argument): should be an integer.";
if not fixp q then rederr
"Error in spcopy_into(fourth argument): should be an integer.";
if p = 0 or q = 0 then
<<
prin2t
"***** Error in spcopy_into: 0 is out of bounds for matrices.";
prin2t
" The top left element is labelled (1,1) and not (0,0).";
return;
>>;
>>;
if not sparsematp(BB) then BB:=sptransmat(BB);
m := sprow_dim(AA);
n := spcol_dim(AA);
r := sprow_dim(BB);
c := spcol_dim(BB);
if not !*fast_la and (r+p-1>m or c+q-1>n) then
<<
% Only print offending matrices if they're not too big.
if m*n<26 and r*c<26 then
<<
prin2t "***** Error in spcopy_into: the matrix";
myspmatpri2(BB);
prin2t " does not fit into";
myspmatpri2(AA);
prin2 " at position ";
prin2 p;
prin2 ",";
prin2 q;
prin2t ".";
return;
>>
else <<
prin2 "***** Error in spcopy_into: first matrix does not fit ";
prin2 " into second matrix at defined position.";
return;
>>;
>>;
a := copy_vect(aa,list('spm,m,n));
for i:=r step -1 until 1 do
<< col:=findrow(bb,i);
if col then
<<for each xx in cdr col do
<< val:=cdr xx;
j:=car xx;
letmtr3(list(a,p+i-1,q+j-1),val,a,nil);
>>;
>>;
>>;
return a;
end;
flag ('(spcopy_into),'opfn);
% Swaps row1 with row2.
symbolic procedure swaprow(ilist,row1,row2,len);
begin scalar r1,r2,rlist,nlist,alist,a,b,cnt,aa,bb,list;
list:=ilist;
r1:=assoc(row1,list);
r2:=assoc(row2,list);
if r1=nil and not (r2 = nil)
then << a:= row1; aa:=(row1 . cdr r2); b:=car r2>>
else if r2 = nil and not (r1=nil)
then << b:=row2; bb:=(row2 . cdr r1); a:=car r1>>
else if not(r1=nil and r2=nil) then <<b:=car r2; a:=car r1>>
else cnt :=len + 1;
cnt := 1;
while not (cnt = len +1) do
<< if list = nil then alist := list(list)
else alist:=car list;
if car alist = a then
<< if r2 = nil then
<< if cnt = b then <<nlist:=aa . nlist;
list := cdr list>>
else if cnt = a then << nlist := nlist;
list:=cdr list>>;
>>
else
<< rlist := (a . cdr r2);
nlist := rlist . nlist;
list := cdr list
>>;
>>
else if car alist = b then
<< if r1 = nil then
<< if cnt = a then <<nlist:=aa . nlist;
list := cdr list>>
else if cnt=b then << nlist := nlist;
list:=cdr list>>;
>>
else
<< rlist := (b . cdr r1);
nlist := rlist . nlist;
list := cdr list
>>;
>>
else if r1=nil and not(cnt = len+1) and cnt = a then
nlist := aa . nlist
else if r2=nil and not (cnt = len+1) and cnt = b then
nlist := bb . nlist
else
<< if alist = '(nil) then nil
else << nlist := alist . nlist;
list := cdr list>>
>>;
cnt := cnt + 1;
>>;
if nlist = nil then return ilist
else return reverse nlist;
end;
symbolic procedure spswap_rows(in_mat,row1,row2);
begin
scalar new_mat,list,pp,r1,r2;
integer rowdim;
list := copy_vect(in_mat,nil);
% if not !*fast_la then use later
<<
if not matrixp in_mat then
rederr "Error in spswap_rows(first argument): should be a matrix.";
rowdim := sprow_dim(in_mat);
if not fixp row1 then rederr
"Error in spswap_rows(second argument): should be an integer.";
if not fixp row2 then
rederr "Error in spswap_rows(third argument): should be an integer.";
if row1>rowdim or row1=0 then rederr
"Error in spswap_rows(second argument): out of range for input matrix.";
if row2>rowdim or row2=0 then rederr
"Error in spswap_rows(third argument): out of range for input matrix.";
>>;
if row1 < row2 then nil else <<pp:=row1; row1:=row2; row2:=pp;>>;
r1:=findrow(list,row1);
r2:=findrow(list,row2);
letmtr3(list(list,row1),r2,list,nil);
letmtr3(list(list,row2),r1,list,nil);
return list;
end;
% Swaps col1 with col2.
symbolic procedure swapcol(ilist,col1,col2,len);
begin scalar c1,c2,rlist,nlist,alist,a,b,aa,bb,cnt,rown,list,row;
cnt := 1;
for i:=len step -1 until 1 do
<< row:=findrow(ilist,i);
if not (row=nil) then
<<
c1:=atsoc(col1,row);
c2:=atsoc(col2,row);
if c1=nil and not (c2 = nil)
then << a:= col1; aa:=(col1 . cdr c2); b:=car c2>>
else if c2 = nil and not (c1=nil)
then << b:=col2; bb:=(col2 . cdr c1); a:=car c1>>
else if not(c1=nil and c2=nil) then <<b:=car c2; a:=car c1>>
else cnt :=len + 1;
rown:=i;
list :=cdr row;
while not (cnt = len + 1) do
<< if list = nil then alist:=list(list)
else alist:=car list;
if car alist = a then
<< if c2 = nil then
<< if cnt = b then <<nlist := bb . nlist;
list := cdr list>>
else if cnt = a then <<nlist := nlist;
list := cdr list >>;
>>
else
<< rlist := (a . cdr c2);
nlist := rlist . nlist;
list := cdr list
>>
>>
else if car alist = b then
<< if c1 = nil then
<< if cnt = a then <<nlist := aa . nlist;
list := cdr list >>
else if cnt = b then <<nlist := nlist;
list := cdr list>>;
>>
else
<< rlist := (b . cdr c1);
nlist := rlist . nlist;
list := cdr list
>>
>>
else if c1 = nil and not(cnt = len+1) and cnt = a then
nlist := aa . nlist
else if c2 = nil and not (cnt = len+1) and cnt = b then
nlist := bb . nlist
else
<< if alist = '(nil) then nil
else
<< nlist := alist . nlist;
list := cdr list>>;
>>;
cnt:=cnt + 1;
>>;
if nlist = nil then letmtr3(list(ilist,rown),list(nil) . list,ilist,nil)
else letmtr3(list(ilist,rown),list(nil) . reverse nlist,ilist,nil);
nlist := nil;
cnt :=1;
>>;
>>;
return ilist;
end;
symbolic procedure spswap_cols(in_mat,col1,col2);
begin
scalar new_mat,list,pp;
integer coldim;
list:=copy_vect(in_mat,nil);
if not !*fast_la then
<<
if not matrixp in_mat then rederr
"Error in spswap_columns(first argument): should be a matrix.";
coldim := spcol_dim(in_mat);
if not fixp col1 then rederr
"Error in spswap_columns(second argument): should be an integer.";
if not fixp col2 then rederr
"Error in spswap_columns(third argument): should be an integer.";
if col1>coldim or col1=0 then rederr
"Error in spswap_columns(second argument): out of range for matrix.";
if col2>coldim or col2=0 then rederr
"Error in spswap_columns(third argument): out of range for input matrix.";
>>;
if col1 < col2 then nil else <<pp:=col1; col1:=col2; col2:=pp;>>;
new_mat := swapcol(list,col1,col2,caddr caddr in_mat);
return new_mat;
end;
% Swaps the two entries in in_mat.
%
% entry1 and entry2 must be lists of the form
% {row position,column position}.
symbolic procedure spswap_entries(in_mat,entry1,entry2);
begin
scalar new_mat;
integer rowdim,coldim,val1,val2;
if not matrixp(in_mat) then rederr
"Error in spswap_entries(first argument): should be a matrix.";
if atom entry1 or car entry1 neq 'list or length cdr entry1 neq 2
then rederr
"Error in spswap_entries(second argument): should be list of 2 elements."
else entry1 := cdr entry1;
if atom entry2 or car entry2 neq 'list or length cdr entry2 neq 2
then rederr
"Error in spswap_entries(third argument): should be a list of 2 elements."
else entry2 := cdr entry2;
if not !*fast_la then
<<
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
if not fixp car entry1 then
<<
prin2 "***** Error in spswap_entries(second argument): ";
prin2t " first element in list must be an integer.";
return;
>>;
if not fixp cadr entry1 then
<<
prin2 "***** Error in spswap_entries(second argument): ";
prin2t " second element in list must be an integer.";
return;
>>;
if car entry1 > rowdim or car entry1 = 0 then
<<
prin2 "***** Error in spswap_entries(second argument): ";
prin2t " first element is out of range for input matrix.";
return;
>>;
if cadr entry1 > coldim or cadr entry1 = 0 then
<<
prin2 "***** Error in spswap_entries(second argument): ";
prin2t " second element is out of range for input matrix.";
return;
>>;
if not fixp car entry2 then
<<
prin2 "***** Error in spswap_entries(third argument): ";
prin2t " first element in list must be an integer.";
return;
>>;
if not fixp cadr entry2 then
<<
prin2 "***** Error in spswap_entries(third argument): ";
prin2t " second element in list must be an integer.";
return;
>>;
if car entry2 > rowdim or car entry2 = 0 then
<<
prin2 "***** Error in spswap_entries(third argument): ";
prin2t " first element is out of range for input matrix.";
return;
>>;
if cadr entry2 > coldim then
<<
prin2 "***** Error in spswap_entries(third argument): ";
prin2t " second element is out of range for input matrix.";
return;
>>;
>>;
new_mat := copy_vect(in_mat,nil);
val1:=findelem2(new_mat,car entry1,cadr entry1);
val2:=findelem2(new_mat,car entry2,cadr entry2);
% if not (val2=0) then
letmtr3(list(new_mat,car entry1,cadr entry1),val2,new_mat,nil);
% if not (val1=0) then
letmtr3(list(new_mat,car entry2,cadr entry2),val1,new_mat,nil);
return new_mat;
end;
flag('(spswap_rows,spswap_cols,spswap_entries),'opfn);
% Takes any number of matrices and joins them horizontally.
%
% Can take either a list of matrices or the matrices as seperate
% arguments.
% This function expands the columns of a matrix in ordere to augment it.
% A Further RE-WRITE function.
% Used to re-create elongated matrices.
symbolic procedure rewrite2(list,num);
begin scalar val,oldcol,newcol,nlist;
for each col in list do
<< val:=cdr col;
oldcol:=car col;
oldcol:=oldcol+num;
newcol:=(oldcol . val);
nlist:=newcol . nlist;
>>;
return reverse nlist;
end;
symbolic procedure expan2(mlist,row,list);
begin scalar rows,cols,rown,newcols,newrows,rlist,cnt,size;
for i:=1:row do
<<cnt:=0;
for each mat1 in mlist do
<< size:=spcol_dim(mat1);
rows:=findrow(mat1,i);
if rows = nil then nil
else << cols := cdr rows;
rown:= i;
if cnt=0 then <<newcols:=append(cols ,newcols);
cnt:=cnt + size>>
else <<
cols:=rewrite2(cols,cnt);
newcols:=append(newcols,cols);
cnt:=cnt + size;
>>;
>>;
>>;
if not (newcols=nil) then
<<
letmtr3(list(list,i),list(nil) . newcols,list,nil);
>>;
newcols:=nil;
>>;
return list;
end;
put('spmatrix_augment,'psopfn,'spmatrix_augment1);
symbolic procedure spmatrix_augment1(matrices);
begin
scalar mat_list,mat1,new_list,he,tl,num,row,col,list;
integer cnt;
if pairp matrices and pairp car matrices and caar matrices = 'list
then matrices := cdar matrices;
if not !*fast_la then
<<
mat_list := for each elt in matrices collect
<<if not sparsematp(list:=reval elt) then sptransmat list
else list>>;
for each elt in mat_list do
if not matrixp(elt) then
rederr "Error in spmatrix_augment: non matrix in input.";
>>;
spconst_rows_test(mat_list);
mat1:=car mat_list;
row:=sprow_dim(mat1);
for each mat1 in mat_list do
<<col:=spcol_dim(mat1);
cnt:=cnt + col;
>>;
col:=cnt;
list:=mkempspmat(row,list('spm,row,col));
new_list:=expan2(mat_list,row,list);
return new_list;
end;
% Takes any number of matrices and joins them vertically.
%
% Can take either a list of matrices or the matrices as seperate
% arguments.
put('spmatrix_stack,'psopfn,'spmatrix_stack1);
symbolic procedure spmatrix_stack1(matrices);
begin
scalar mat_list,new_list,he,tl,nam,row,col,list;
integer cnt;
if pairp matrices and pairp car matrices and caar matrices = 'list
then matrices := cdar matrices;
if not !*fast_la then
<<
mat_list := for each elt in matrices collect
<<if not sparsematp(list:=reval elt) then sptransmat list
else list>>;
for each elt in mat_list do
if not matrixp(elt) then
rederr "Error in spmatrix_stack: non matrix in input.";
>>;
spconst_columns_test(mat_list);
col:=spcol_dim(car mat_list);
for each mat1 in mat_list do
<< row:=sprow_dim(mat1);
cnt := cnt + row;
>>;
row:=cnt;
new_list:=mkempspmat(row,list('spm,row,col));
cnt:=1;
for each mat1 in mat_list do
<< row:=sprow_dim(mat1);
for i:=1:row do
<< he:=findrow(mat1,i);
if he then
letmtr3(list(new_list,cnt),he,new_list,nil);
cnt:=cnt+1;
>>;
>>;
return new_list;
end;
symbolic procedure spconst_rows_test(mat_list);
%
% Tests that each matrix in mat_list has the same number of rows
% (otherwise augmentation not possible).
%
begin
integer i,listlen,rowdim;
listlen := length(mat_list);
rowdim := sprow_dim(car mat_list);
i := 1;
while i<listlen and sprow_dim(car mat_list) = sprow_dim(cadr mat_list)
do << i := i+1; mat_list := cdr mat_list; >>;
if i=listlen then return rowdim
else
<<
prin2 "***** Error in spmatrix_augment: ";
rederr "all input matrices must have the same row dimension.";
>>;
end;
symbolic procedure spconst_columns_test(mat_list);
%
% Tests that each matrix in mat_list has the same number of columns
% (otherwise stacking not possible).
%
begin
integer i,listlen,coldim;
listlen := length(mat_list);
coldim := spcol_dim(car mat_list);
i := 1;
while i<listlen and spcol_dim(car mat_list) =
spcol_dim(cadr mat_list)
do << i := i+1; mat_list := cdr mat_list; >>;
if i=listlen then return coldim
else
<<
prin2 "***** Error in spmatrix_stack: ";
rederr
"all input matrices must have the same column dimension.";
return;
>>;
end;
% Extends in_mat by rows rows (!) and cols columns. New entries are
% initialised to entry.
symbolic procedure spextend(in_mat,rows,cols,entry);
begin
scalar ex_mat;
integer rowdim,coldim,i,j;
if not matrixp(in_mat) then
rederr "Error in spextend(first argument): should be a matrix.";
if not fixp rows then
rederr "Error in spextend(second argument): should be an integer.";
if not fixp cols then
rederr "Error in spextend(third argument): should be an integer.";
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
ex_mat := mkempspmat(rowdim+rows,list('smp,rowdim+rows,coldim+cols));
ex_mat := spcopy_into(in_mat,ex_mat,1,1);
for i:=rowdim+rows step -1 until rowdim+1 do
<<
for j:=coldim+cols step -1 until coldim+1 do
<<
letmtr3(list(ex_mat,i,j),entry,ex_mat,nil);
>>;
>>;
return ex_mat;
end;
flag('(spextend),'opfn);
% Can take either a list of arguments or the arguments seperately.
%
% Takes any number of either scalar entries or square matrices and
% creates the diagonal.
put('spdiagonal,'psopfn,'spdiagonal1); % To allow variable input.
symbolic procedure spdiagonal1(mat_list);
begin
scalar diag_mat;
if pairp mat_list and pairp car mat_list and caar mat_list = 'list
then mat_list := cdar mat_list;
mat_list := for each elt in mat_list collect
<< if not (sparsematp(aeval elt)) and not numberp elt
then sptransmat elt
else reval elt
>>;
for each elt in mat_list do
<<
if matrixp(elt) and not squarep(elt) then
<<
% Only print offending matrix if it's not too big.
if sprow_dim(elt)<5 or spcol_dim(elt)> 5 then
<<
prin2t "***** Error in spdiagonal: ";
myspmatpri2(elt);
prin2t " is not a square matrix.";
rederr "";
>>
else
rederr "Error in spdiagonal: input contains non square matrix.";
>>;
>>;
diag_mat := spdiag({mat_list});
return diag_mat;
end;
symbolic procedure spdiag(uu);
%
% Takes square or scalar matrix entries and creates a matrix with
% these matrices on the diagonal.
%
begin
scalar bigA,arg,input,u,val,a,b,col,j;
integer nargs,Aidx,stp,bigsize,smallsize;
u := car uu;
input := u;
bigsize:=0;
nargs:=length input;
for i:=1:nargs do
<<
arg:=car input;
% If scalar entry.
if algebraic length(arg) = 1 then bigsize:=bigsize+1
else
<<
bigsize:=bigsize+sprow_dim(arg);
>>;
input := cdr input;
>>;
bigA := mkempspmat(bigsize,list('spm,bigsize,bigsize));
Aidx:=1;
input := u;
for k:=1:nargs do
<<
arg:=car input;
% If scalar entry.
if algebraic length(arg) = 1 then
<<
letmtr3(list(bigA,Aidx,Aidx),arg,bigA,nil);
Aidx:=Aidx+1;
input := cdr input;
>>
else
<<
smallsize:= sprow_dim(arg);
stp:=smallsize+Aidx-1;
a:=1;
for i:=Aidx:stp do
<< col:=findrow(arg,a);
if col then
<< for each xx in cdr col do
<< val:=cdr xx;
j:=(Aidx-1)+car xx;
letmtr3(list(bigA,i,j),val,bigA,nil);
>>;
>>;
a:=a+1;
>>;
Aidx := Aidx+smallsize;
input := cdr input;
>>;
>>;
return biga;
end;
% Replaces row2 (r2) by mult1*r1 + r2.
symbolic procedure spadd_rows(in_mat,r1,r2,mult1);
begin
scalar new_mat,val,val1,val2,row1,row2;
integer i,rowdim,coldim;
coldim := spcol_dim(in_mat);
if not !*fast_la then
<<
if not matrixp in_mat then
rederr "Error in spadd_rows(first argument): should be a matrix.";
rowdim := sprow_dim(in_mat);
if not fixp r1 then
rederr "Error in spadd_rows(second argument): should be an integer.";
if not fixp r2 then
rederr "Error in spadd_rows(third argument): should be an integer.";
if r1>rowdim or r1=0 then rederr
"Error in spadd_rows(second argument): out of range for input matrix.";
if r2>rowdim or r2=0 then rederr
"Error in spadd_rows(third argument): out of range for input matrix.";
>>;
new_mat := copy_vect(in_mat,nil);
% Efficiency.
if (my_reval mult1) = 0 then return new_mat;
row1:=findrow(in_mat,r1);
row2:=findrow(in_mat,r2);
for each xx in cdr row1 do
<< i:=car xx;
val1:=cdr xx;
val2:=atsoc(i,row2);
val:=reval {'times,mult1,val1};
if val2 then
<<val:=reval {'plus,val,cdr val2};
if not (val=0) then letmtr3(list(new_mat,r2,i),val,new_mat,nil);
>>
else letmtr3(list(new_mat,r2,i),val,new_mat,nil);
>>;
return new_mat;
end;
% Replaces column2 (c2) by mult1*c1 + c2.
symbolic procedure spadd_columns(in_mat,c1,c2,mult1);
begin
scalar new_mat,val;
integer i,rowdim,coldim;
rowdim := sprow_dim(in_mat);
if not !*fast_la then
<<
if not matrixp in_mat then
rederr "Error in spadd_columns(first argument): should be a matrix.";
coldim := spcol_dim(in_mat);
if not fixp c1 then rederr
"Error in spadd_columns(second argument): should be an integer.";
if not fixp c2 then rederr
"Error in spadd_columns(third argument): should be an integer.";
if c1>coldim or c1=0 then rederr
"Error in spadd_columns(second argument): out of range for input matrix.";
if c2>rowdim or c2=0 then rederr
"Error in spadd_columns(third argument): out of range for input matrix.";
>>;
new_mat := copy_vect(in_mat,nil);
if (my_reval mult1) = 0 then return new_mat;
for i:=1:rowdim do
<<val:=reval {'plus,{'times,mult1,
findelem2(new_mat,i,c1)},findelem2(in_mat,i,c2)};
if not (val=0) then
letmtr3(list(new_mat,i,c2),val,new_mat,nil);
>>;
return new_mat;
end;
flag('(spadd_rows,spadd_columns),'opfn);
% Adds value to each element in each row in row_list.
%
% row_list can be either an integer or a list of integers.
symbolic procedure spadd_to_rows(in_mat,row_list,value);
begin
scalar new_mat,col,val;
integer i,rowdim,coldim;
if not matrixp in_mat then
rederr "Error in spadd_to_row(first argument): should be a matrix.";
if atom row_list then row_list := {row_list}
else if car row_list = 'list then row_list := cdr row_list
else
<<
prin2 "***** Error in spadd_to_rows(second argument): ";
prin2t " should be either integer or a list of integers.";
return;
>>;
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
new_mat := copy_vect(in_mat,nil);
for each row in row_list do
<<
if not fixp row then rederr
"Error in spadd_to_row(second argument): should be an integer.";
if row>rowdim or row=0 then
<<
prin2 "***** Error in spadd_to_rows(second argument): ";
rederr "contains row which is out of range for input matrix.";
>>;
>>;
for each row in row_list do
<<for i:=coldim step -1 until 1 do
letmtr3(list(new_mat,row,i),value,new_mat,nil);
col:=findrow(in_mat,row);
if col then
<<for each xx in cdr col do
<<i:=car xx;
val:=cdr xx;
letmtr3(list(new_mat,row,i),reval {'plus,val,value},new_mat,nil);
>>;
>>;
>>;
return new_mat;
end;
symbolic procedure spadd_to_columns(in_mat,col_list,value);
%
% Adds value to each element in each column in col_list.
%
% col_list can be either an integer or a list of integers.
%
begin
scalar new_mat,col,val;
integer i,rowdim,coldim;
if not matrixp in_mat then rederr
"Error in spadd_to_columns(first argument): should be a matrix.";
if atom col_list then col_list := {col_list}
else if car col_list = 'list then col_list := cdr col_list
else
<<
prin2 "***** Error in spadd_to_columns(second argument): ";
prin2t " should be either integer or list of integers.";
return;
>>;
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
new_mat := copy_vect(in_mat,nil);
for each col in col_list do
<<
if not fixp col then rederr
"Error in spadd_to_columns(second argument): should be an integer.";
if col>coldim or col=0 then
<<
prin2 "***** Error in spadd_to_columns(second argument): ";
rederr
"contains column which is out of range for input matrix.";
>>;
>>;
for i:=1:rowdim do
<<col:=findrow(in_mat,i);
for each xx in col_list do
<<val:=atsoc(xx,col);
if val then
<<letmtr3(list(new_mat,i,xx),reval
{'plus,cdr val,value},new_mat,nil);
>>
else letmtr3(list(new_mat,i,xx),value,new_mat,nil);
>>;
>>;
return new_mat;
end;
flag('(spadd_to_rows,spadd_to_columns),'opfn);
% Replaces rows specified in row_list by row * mult1.
symbolic procedure spmult_rows(in_mat,row_list,mult1);
begin
scalar new_mat,col;
integer i,rowdim,coldim,val;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in spmult_rows(first argument): should be a matrix.";
if atom row_list then row_list := {row_list}
else if car row_list = 'list then row_list := cdr row_list;
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
new_mat := copy_vect(in_mat,nil);
for each row in row_list do
<<
if not !*fast_la and not fixp row then rederr
"Error in spmult_rows(second argument): contains non integer.";
if not !*fast_la and (row>rowdim or row=0) then
<<
prin2 "***** Error in spmult_rows(second argument): ";
rederr "contains row that is out of range for input matrix.";
>>;
col:=findrow(in_mat,row);
if col then
<<for each xx in cdr col do
<<i:=car xx;
val:=cdr xx;
letmtr3(list(new_mat,row,i),reval{'times,mult1,val},
new_mat,nil);
>>;
>>;
>>;
return new_mat;
end;
% Replaces columns specified in column_list by column * mult1.
symbolic procedure spmult_columns(in_mat,column_list,mult1);
begin
scalar new_mat,col;
integer i,rowdim,coldim,val;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in spmult_columns(first argument): should be a matrix.";
if atom column_list then column_list := {column_list}
else if car column_list = 'list then column_list := cdr column_list;
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
new_mat := copy_vect(in_mat,nil);
for each column in column_list do
<<
if not !*fast_la and not fixp column then rederr
"Error in spmult_columns(second argument): contains non integer.";
if not !*fast_la and (column>coldim or column=0) then
<<
prin2 "***** Error in spmult_columns(second argument): ";
rederr "contains column that is out of range for input matrix.";
>>;
>>;
for i:=1:rowdim do
<< col:=findrow(in_mat,i);
if col then
<<for each xx in column_list do
<< val:=atsoc(xx,col);
if val then
letmtr3(list(new_mat,i,xx),reval
{'times,mult1,cdr val}, new_mat,nil);
>>;
>>;
>>;
return new_mat;
end;
flag('(spmult_rows,spmult_columns),'opfn);
% Create characteristic matrix. ie: C := lmbda*I - in_mat.
% in_ mat must be square.
symbolic procedure spchar_matrix(in_mat,lmbda);
begin
scalar charmat;
integer rowdim;
if not matrixp(in_mat) then
rederr "Error in spchar_matrix(first argument): should be a matrix.";
if not squarep(in_mat) then rederr
"Error in spchar_matrix(first argument): must be a square matrix.";
rowdim := sprow_dim(in_mat);
charmat := {'plus,{'times,lmbda,spmake_identity(rowdim)},
{'minus,in_mat}};
return charmat;
end;
% Finds characteristic polynomial of matrix in_mat.
% ie: det(lmbda*I - in_mat).
symbolic procedure spchar_poly(in_mat,lmbda);
begin
scalar chpoly,carmat;
if not matrixp(in_mat) then
rederr "Error in spchar_poly(first argument): should be a matrix.";
carmat := spchar_matrix(in_mat,lmbda);
chpoly := algebraic det(carmat);
return chpoly;
end;
flag('(spchar_matrix,spchar_poly),'opfn);
% Computes the Hermitian transpose (HT say) of in_mat.
%
% The (i,j)'th element of HT = conjugate of the (j,i)'th element of
% in__mat.
symbolic procedure sphermitian_tp(in_mat);
begin
scalar h_tp,element;
integer ii,row,col;
if not matrixp(in_mat) then
rederr "Error in sphermitian_tp: non matrix input.";
h_tp := algebraic tp(in_mat);
for row:=1:sprow_dim(h_tp) do
<<col:=findrow(h_tp,row);
if col then
<<for each xx in cdr col do
<<ii:=car xx;
element := cdr xx;
letmtr3(list(h_tp,row,ii),
algebraic (repart(element) - i*impart(element)),h_tp,nil);
>>;
>>;
>>;
return h_tp;
end;
flag('(sphermitian_tp),'opfn);
% Removes the sub_matrix from A consisting of the rows in row_list and
% the columns in col_list. (Both row_list and col_list can be single
% integer values).
symbolic procedure spsub_matrix(A,row_list,col_list);
begin
scalar new_mat;
if not !*fast_la and not matrixp(A) then rederr
"Error in spsub_matrix(first argument): should be a matrix.";
new_mat := spstack_rows(A,row_list);
new_mat := spaugment_columns(new_mat,col_list);
return new_mat;
end;
flag('(spsub_matrix),'opfn);
% Converts all elements in pivot column (apart from the one in pivot
% row) to 0.
symbolic procedure sppivot(in_mat,pivot_row,pivot_col);
begin
scalar piv_mat,ratio,val,col,val1,val2;
integer i,j,rowdim,coldim;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in sppivot(first argument): should be a matrix.";
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
if not !*fast_la then
<<
if not fixp pivot_row then
rederr "Error in sppivot(second argument): should be an integer.";
if pivot_row>rowdim or pivot_row=0 then rederr
"Error in sppivot(second argument): out of range for input matrix.";
if not fixp pivot_col then rederr
"Error in sppivot(third argument): should be an integer.";
if pivot_col>coldim or pivot_col=0 then rederr
"Error in sppivot(third argument): out of range for input matrix.";
if findelem2(in_mat,pivot_row,pivot_col) = 0 then
rederr "Error in sppivot: cannot pivot on a zero entry.";
>>;
piv_mat := copy_vect(in_mat,nil);
val2:=findelem2(in_mat,pivot_row,pivot_col);
for i:=1:rowdim do
<< col:=findrow(in_mat,i);
val1:=atsoc(pivot_col,col);
if val1 then val1:=cdr val1;
ratio := reval {'quotient,val1,val2};
if col then
<<for each xx in cdr col do
<<j:=car xx;
val1:=cdr xx;
if i = pivot_row then <<>>
else
<<val:=reval {'plus,val1,
{'minus,{'times,ratio,findelem2(in_mat,pivot_row,j)}
}};
letmtr3(list(piv_mat,i,j),val,piv_mat,nil);
>>;
>>;
>>;
>>;
return piv_mat;
end;
% Same as pivot but only rows a .. to .. b, where row_list = {a,b},
% are changed.
%
% rows_pivot will work if row_list is just an integer.
symbolic procedure sprows_pivot(in_mat,pivot_row,pivot_col,row_list);
begin
scalar piv_mat,ratio,val,col,val1,val2;
integer j,rowdim,coldim;
rowdim := sprow_dim(in_mat);
coldim := spcol_dim(in_mat);
if not !*fast_la then
<<
if not matrixp(in_mat) then
rederr "Error in sprows_pivot(first argument): should be a matrix.";
if not fixp pivot_row then
rederr "Error in sprows_pivot(second argument): should be an integer.";
if pivot_row>rowdim or pivot_row=0 then rederr
"Error in sprows_pivot(second argument): out of range for input matrix.";
if not fixp pivot_col then
rederr "Error in sprows_pivot(third argument): should be an integer.";
if pivot_col>coldim or pivot_col=0 then rederr
"Error in sprows_pivot(third argument): out of range for input matrix.";
>>;
if atom row_list then row_list := {row_list}
else if pairp row_list and car row_list = 'list
then row_list := cdr row_list
else
<<
prin2 "***** Error in sprows_pivot(fourth argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
if findelem2(in_mat,pivot_row,pivot_col) = 0 then
rederr "Error in sprows_pivot: cannot pivot on a zero entry.";
piv_mat := copy_vect(in_mat,nil);
val2:=findelem2(in_mat,pivot_row,pivot_col);
for each elt in row_list do
<<
if not !*fast_la then
<<
if not fixp elt then rederr
"Error in sprows_pivot: fourth argument contains a non integer.";
if elt>rowdim or elt=0 then
<<
prin2 "***** Error in sprows_pivot(fourth argument): ";
rederr "contains row which is out of range for input matrix.";
>>;
>>;
if elt = pivot_row then nil
else ratio := reval {'quotient,findelem2(in_mat,elt,pivot_col),val2};
col:=findrow(in_mat,elt);
if col then
<<for each xx in cdr col do
<<j:=car xx;
val1:=cdr xx;
val:=reval {'plus,val1, {'minus,
{'times,ratio,findelem2(in_mat,pivot_row,j)}}};
letmtr3(list(piv_mat,elt,j),val,piv_mat,nil);
>>;
>>;
>>;
return piv_mat;
end;
flag('(sppivot,sprows_pivot),'opfn);
% jacobian(exp,var) computes the Jacobian matrix of exp w.r.t. var.
% The (i,j)'th entry is diff(nth(exp,i),nth(var,j)).
symbolic procedure spjacobian(exp_list,var_list);
begin
scalar jac,exp1,var1,val;
integer i,j,rowdim,coldim;
if atom exp_list then exp_list := {exp_list}
else if car exp_list neq 'list then rederr
"Error in spjacobian(first argument): expressions must be in a list."
else exp_list := cdr exp_list;
if atom var_list then var_list := {var_list}
else if car var_list neq 'list then rederr
"Error in jacobian(second argument): variables must be in a list."
else var_list := cdr var_list;
rowdim := length exp_list;
coldim := length var_list;
jac := mkempspmat(rowdim,list('spm,rowdim,coldim));
for i:=1:rowdim do
<<
for j:=1:coldim do
<<
exp1 := nth(exp_list,i);
var1 := nth(var_list,j);
val:= algebraic df(exp1,var1);
if val = 0 then nil
else letmtr3(list(jac,i,j),val,jac,nil);
>>;
>>;
return jac;
end;
flag('(spjacobian),'opfn);
% variables can be either a list or a single variable.
%
% A Hessian matrix is a matrix whose (i,j)'th entry is
% df(df(poly,nth(var,i)),nth(var,j))
%
% where df is the derivative.
symbolic procedure sphessian(poly,variables);
begin
scalar hess_mat,part1,part2,elt;
integer row,col,sq_size;
if atom variables then variables := {variables}
else if car variables = 'list then variables := cdr variables
else
<<
prin2 "***** Error in sphessian(second argument): ";
prin2t
" should be either a single variable or a list of variables.";
return;
>>;
sq_size := length variables;
hess_mat := mkempspmat(sq_size,list('spm,sq_size,sq_size));
for row:=1:sq_size do
<<
for col:=1:sq_size do
<<
part1 := nth(variables,row);
part2 := nth(variables,col);
elt := algebraic df(df(poly,part1),part2);
if elt = 0 then nil
else letmtr3(list(hess_mat,row,col),elt,hess_mat,nil);
>>;
>>;
return hess_mat;
end;
flag('(sphessian),'opfn);
% Given the system of linear equations, coeff_matrix returns {A,X,b}
% s.t. AX = b.
%
% Input can be either a list of linear equations or the linear
% equations as individual arguments.
% To allow variable input.
put('spcoeff_matrix,'psopfn,'spcoeff_matrix1);
symbolic procedure spcoeff_matrix1(equation_list);
begin
scalar variable_list,A,X,b;
if pairp car equation_list and caar equation_list = 'list then
equation_list := cdar equation_list;
equation_list := remove_equals(equation_list);
variable_list := get_variable_list(equation_list);
if variable_list = nil then
rederr "Error in spcoeff_matrix: no variables in input.";
check_linearity(equation_list,variable_list);
A := spget_A(equation_list,variable_list);
X := spget_X(variable_list);
b := spget_b(equation_list,variable_list);
return {'list,A,X,b};
end;
symbolic procedure remove_equals(equation_list);
%
% If any of the equations are equalities the equalities are removed
% to leave a list of polynomials.
%
begin
equation_list := for each equation in equation_list collect
if pairp equation and car equation = 'equal then
reval{'plus,cadr equation,{'minus,caddr equation}}
else equation;
return equation_list;
end;
symbolic procedure get_variable_list(equation_list);
%
% Gets hold of all variables from the equations in equation_list.
%
begin
scalar variable_list;
for each equation in equation_list do
variable_list := union(get_coeffs(equation),variable_list);
return reverse variable_list;
end;
symbolic procedure check_linearity(equation_list,variable_list);
%
% Checks that we really are dealing with a system of linear equations.
%
for each equation in equation_list do
<<
for each variable in variable_list do
<<
if deg(equation,variable) > 1 then
rederr "Error in spcoeff_matrix: the equations are not linear.";
>>;
>>;
symbolic procedure spget_A(equation_list,variable_list);
begin
scalar A,element,var_elt,val;
integer row,col,length_equation_list,length_variable_list;
length_equation_list := length equation_list;
length_variable_list := length variable_list;
A := mkempspmat(length equation_list,
list('spm,length equation_list,length variable_list));
for row:=1:length_equation_list do
<<
for col:=1:length_variable_list do
<<
element := nth(equation_list,row);
var_elt := nth(variable_list,col);
val:=algebraic coeffn(element,var_elt,1);
if val = 0 then nil
else letmtr3(list(A,row,col),val,A,nil);
>>;
>>;
return A;
end;
symbolic procedure spget_b(equation_list,variable_list);
%
% Puts the integer parts of all the equations into a column matrix.
%
begin
scalar substitution_list,integer_list,b;
integer length_integer_list,row;
substitution_list :=
'list.for each variable in variable_list collect
{'equal,variable,0};
integer_list := for each equation in equation_list collect
algebraic sub(substitution_list,equation);
length_integer_list := length integer_list;
b := mkempspmat(length_integer_list,list('spm,length_integer_list,1));
for row:=1:length_integer_list do
letmtr3(list(b,row,1),-nth(integer_list,row),b,nil);
return b;
end;
symbolic procedure spget_X(variable_list);
begin
scalar X;
integer row,length_variable_list;
length_variable_list := length variable_list;
X := mkempspmat(length_variable_list,list('spm,length_variable_list,1));
for row := 1:length variable_list do
letmtr3(list(X,row,1),nth(variable_list,row),x,nil);
return X;
end;
symbolic procedure get_coeffs(poly);
%
% Gets all kernels in a poly.
%
begin
scalar ker_list_num,ker_list_den;
ker_list_num := kernels !*q2f simp reval num poly;
ker_list_den := kernels !*q2f simp reval den poly;
ker_list_num := union(ker_list_num,ker_list_den);
return ker_list_num;
end;
% Takes as input a monic univariate polynomial in a variable x.
% Returns a companion matrix associated with the polynomial poly(x).
%
% If C := companion(p,x) and p is a0+a1*x+...+x^n (a univariate monic
% polynomial), them C(i,n) = -coeff(p,x,i-1), C(i,i-1) = 1 (i=2..n)
% and C(i,j) = 0 for all other i and j.
symbolic procedure spcompanion(poly,x);
begin
scalar mat1;
integer n,val;
n := deg(poly,x);
if my_reval coeffn(poly,x,n) neq 1 then msgpri
("Error in spcompanion(first argument): Polynomial",
poly, "is not monic.",nil,t);
mat1 := mkempspmat(n,list('smp,n,n));
val:=coeffn(poly,x,0);
if val=0 then nil
else letmtr3(list(mat1,1,n),{'minus,val},mat1,nil);
for i:=2:n do
<<
letmtr3(list(mat1,i,i-1),1,mat1,nil);
>>;
for j:=2:n do
<< val:=coeffn(poly,x,j-1);
if val = 0 then nil
else letmtr3(list(mat1,j,n),{'minus,val},mat1,nil);
>>;
return mat1;
end;
% Given a companion matrix, find_companion will return the associated
% polynomial.
symbolic procedure spfind_companion(R,x);
begin
scalar p;
integer rowdim,k;
if not matrixp(R) then rederr
{"Error in spfind_companion(first argument): should be a matrix."};
rowdim := sprow_dim(R);
k := 2;
while k<=rowdim and findelem2(R,k,k-1)=1 do k:=k+1;
p := 0;
for j:=1:k-1 do
<<
p:={'plus,p,{'times,{'minus,findelem2(R,j,k-1)},{'expt,x,j-1}}};
>>;
p := {'plus,p,{'expt,x,k-1}};
return p;
end;
flag('(spcompanion,spfind_companion),'opfn);
endmodule;
end;
%***********************************************************************
%=======================================================================
%
% End of Code.
%
%=======================================================================
%***********************************************************************
%in "splu_decomp.red";
%in "spsvd.red";
%in "spcholesky.red";
%in "spgramshm.red";