Artifact bd29fe1894c9b0b91f910abc74341d6152271ea87eca374a23a14594a2d2a292:
- Executable file
r37/packages/sparse/splinalg.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 65754) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/sparse/splinalg.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 65754) [annotate] [blame] [check-ins using]
%*********************************************************************** %======================================================================= % % 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";