module linalg; % The Linear Algebra package.
%**********************************************************************%
% %
% Author: Matt Rebbeck, March-July 1994 (at ZIB). %
% Modifications by: Walter Tietze.
% %
% Please report bugs to: Winfried Neun, %
% Konrad-Zuse-Zentrum %
% fuer Informationstechnik Berlin %
% Heilbronner Str. 10 %
% 10711 Berlin - Wilmersdorf %
% Federal Republic of Germany %
% %
% (email) neun@sc.ZIB-Berlin.de %
% %
% %
% %
% This package provides a selection of useful functions in the field %
% of linear algebra: %
% %
% add_columns add_rows add_to_columns add_to_rows %
% augment_columns band_matrix block_matrix char_matrix %
% char_poly cholesky coeff_matrix column_dim %
% companion copy_into diagonal extend %
% get_columns get_rows gram_schmidt hermitian_tp %
% hessian hilbert jacobian jordan_block %
% lu_decom make_identity matrix_augment matrixp %
% matrix_stack minor mult_column mult_row %
% pivot pseudo_inverse random_matrix remove_columns %
% remove_rows row_dim rows_pivot simplex %
% squarep stack_rows sub_matrix svd %
% swap_columns swap_entries swap_rows symmetricp %
% toeplitz vandermonde kronecker_product %
% %
% %
% %
% The package implements the following switches: %
% %
% imaginary \ %
% lower_matrix \ %
% not_negative ) for details see the random_matrix comments. %
% only_integer / %
% upper_matrix / %
% %
% fast_la (see below). %
% %
% %
% %
% For further details about the linear algebra package see the %
% linear_algebra.tex file. %
% %
% %
% %
% NB: The functions in this package are written to be used from the %
% user level. Some of them may well need a bit of fiddling with to get %
% them to work from symbolic mode. %
% %
%**********************************************************************%
load_package matrix;
create!-package('(linalg lamatrix gramschm ludecom cholesky svd
simplex tadjoint), '(contrib linalg));
switch fast_la; % If ON, then the following functions will be faster:
% add_columns add_rows augment_columns column_dim %
% copy_into make_identity matrix_augment matrix_stack %
% minor mult_column mult_row pivot %
% remove_columns remove_rows rows_pivot squarep %
% stack_rows sub_matrix swap_columns swap_entries %
% swap_rows symmetricp %
% 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.
symbolic smacro procedure my_reval(n);
%
% Only revals if it needs to.
%
if fixp(n) then n else reval(n);
symbolic procedure swap_elt(in_list,elt1,elt2);
%
% Swaps elt elt1 with elt elt2 in in_list.
%
% NB: destructive.
%
begin
scalar bucket;
bucket := nth(in_list,elt1);
nth(in_list,elt1) := nth(in_list,elt2);
nth(in_list,elt2) := bucket;
end;
symbolic procedure row_dim(in_mat);
%
% Finds row dimension of a matrix.
%
begin
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in row_dim: input should be a matrix.";
return first size_of_matrix(in_mat);
end;
symbolic procedure column_dim(in_mat);
%
% Finds column dimension of a matrix.
%
begin
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in column_dim: input should be a matrix.";
return second size_of_matrix(in_mat);
end;
flag('(row_dim,column_dim),'opfn);
symbolic procedure matrixp(A);
%
% Tests if input is a matrix (boolean).
%
if not eqcar(A,'mat) then nil else t;
flag('(matrixp),'boolean);
flag('(matrixp),'opfn);
symbolic procedure size_of_matrix(A);
%
% Takes matrix and returns list {no. of rows, no. of columns}.
%
begin
integer row_dim,column_dim;
row_dim := -1 + length A;
column_dim := length cadr A;
return {row_dim,column_dim};
end;
symbolic procedure companion(poly,x);
%
% 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.
%
begin
scalar mat1;
integer n;
n := deg(poly,x);
if my_reval coeffn(poly,x,n) neq 1 then msgpri
("Error in companion(first argument): Polynomial",
poly, "is not monic.",nil,t);
mat1 := mkmatrix(n,n);
setmat(mat1,1,n,{'minus,coeffn(poly,x,0)});
for i:=2:n do
<<
setmat(mat1,i,i-1,1);
>>;
for j:=2:n do
<<
setmat(mat1,j,n,{'minus,coeffn(poly,x,j-1)});
>>;
return mat1;
end;
symbolic procedure find_companion(R,x);
%
% Given a companion matrix, find_companion will return the associated
% polynomial.
%
begin
scalar p;
integer rowdim,k;
if not matrixp(R) then rederr
{"Error in find_companion(first argument): should be a matrix."};
rowdim := row_dim(R);
k := 2;
while k<=rowdim and getmat(R,k,k-1)=1 do k:=k+1;
p := 0;
for j:=1:k-1 do
<<
p:={'plus,p,{'times,{'minus,getmat(R,j,k-1)},{'expt,x,j-1}}};
>>;
p := {'plus,p,{'expt,x,k-1}};
return p;
end;
flag('(companion,find_companion),'opfn);
symbolic procedure jordan_block(const,mat_dim);
%
% Takes a constant (const) and an integer (mat_dim) and creates
% a jordan block of dimension mat_dim x mat_dim.
%
begin
scalar JB;
if not fixp mat_dim then rederr
"Error in jordan_block(second argument): should be an integer.";
JB := mkmatrix(mat_dim,mat_dim);
for i:=1:mat_dim do
<<
for j:=1:mat_dim do
<<
if i=j then
<<
setmat(JB,i,j,const);
if i<mat_dim then setmat(JB,i,j+1,1);
>>;
>>;
>>;
return JB;
end;
flag ('(jordan_block),'opfn);
symbolic procedure sub_matrix(A,row_list,col_list);
%
% 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).
%
begin
scalar new_mat;
if not !*fast_la and not matrixp(A) then rederr
"Error in sub_matrix(first argument): should be a matrix.";
new_mat := stack_rows(A,row_list);
new_mat := augment_columns(new_mat,col_list);
return new_mat;
end;
% flag('(sub_matrix),'opfn);
rtypecar sub_matrix;
symbolic procedure copy_into(BB,AA,p,q);
%
% Copies matrix BB into AA with BB(1,1) at AA(p,q).
%
begin
scalar A,B;
integer m,n,r,c;
if not !*fast_la then
<<
if not matrixp(BB) then rederr
"Error in copy_into(first argument): should be a matrix.";
if not matrixp(AA) then rederr
"Error in copy_into(second argument): should be a matrix.";
if not fixp p then rederr
"Error in copy_into(third argument): should be an integer.";
if not fixp q then rederr
"Error in copy_into(fourth argument): should be an integer.";
if p = 0 or q = 0 then
<<
prin2t
"***** Error in copy_into: 0 is out of bounds for matrices.";
prin2t
" The top left element is labelled (1,1) and not (0,0).";
return;
>>;
>>;
m := row_dim(AA);
n := column_dim(AA);
r := row_dim(BB);
c := column_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 copy_into: the matrix";
matpri(BB);
prin2t " does not fit into";
matpri(AA);
prin2 " at position ";
prin2 p;
prin2 ",";
prin2 q;
prin2t ".";
return;
>>
else
<<
prin2 "***** Error in copy_into: first matrix does not fit ";
prin2 " into second matrix at defined position.";
return;
>>;
>>;
A := mkmatrix(m,n);
B := mkmatrix(r,c);
for i:=1:m do
<<
for j:=1:n do
<<
setmat(A,i,j,getmat(AA,i,j));
>>;
>>;
for i:=1:r do
<<
for j:=1:c do
<<
setmat(B,i,j,getmat(BB,i,j));
>>;
>>;
for i:=1:r do
<<
for j:=1:c do
<<
setmat(A,p+i-1,q+j-1,getmat(B,i,j));
>>;
>>;
return A;
end;
flag ('(copy_into),'opfn);
symbolic procedure copy_mat(u);
if pairp u then cons (copy_mat car u, copy_mat cdr u) else u;
put('diagonal,'psopfn,'diagonal1); % To allow variable input.
symbolic procedure diagonal1(mat_list);
%
% 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.
%
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 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 row_dim(elt)<5 or column_dim(elt)> 5 then
<<
prin2t "***** Error in diagonal: ";
matpri(elt);
prin2t " is not a square matrix.";
rederr "";
>>
else
rederr "Error in diagonal: input contains non square matrix.";
>>;
>>;
diag_mat := diag({mat_list});
return diag_mat;
end;
symbolic procedure diag(uu);
%
% Takes square or scalar matrix entries and creates a matrix with
% these matrices on the diagonal.
%
% What a horrible way to do it!
%
begin
scalar bigA,arg,input,u;
integer nargs,n,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 or eqcar(arg,'quotient)
then bigsize:=bigsize+1
else
<<
bigsize:=bigsize+row_dim(arg);
>>;
input := cdr input;
>>;
bigA := mkmatrix(bigsize,bigsize);
Aidx:=1;
input := u;
for k:=1:nargs do
<<
arg:=car input;
% If scalar entry.
if algebraic length(arg) = 1 or eqcar(arg,'quotient) then
<<
setmat(bigA,Aidx,Aidx,arg);
Aidx:=Aidx+1;
input := cdr input;
>>
else
<<
smallsize:= row_dim(arg);
stp:=smallsize+Aidx-1;
for i:=Aidx:stp do
<<
for j:=Aidx:stp do
<<
arg:=car input;
% Find (i-Aidx+1)'th row.
arg := cdr arg;
<<
n:=1;
while n < (i-Aidx+1) do
<<
arg := cdr arg;
n:=n+1;
>>;
>>;
arg := car arg;
%
% Find (j-Aidx+1)'th column elt of i'th row.
%
<<
n:=1;
while n < (j-Aidx+1) do
<<
arg := cdr arg;
n:=n+1;
>>;
>>;
arg := car arg;
setmat(bigA,i,j,arg);
>>;
>>;
Aidx := Aidx+smallsize;
input := cdr input;
>>;
>>;
return biga;
end;
symbolic procedure band_matrix(elt_list,sq_size);
%
% 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.
%
begin
scalar band_matrix;
integer i,j,no_elts,middle_pos;
if not fixp sq_size then rederr
"Error in band_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 band_matrix(first argument): should be single value or list.";
no_elts := length elt_list;
if evenp no_elts then rederr
"Error in band 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 band_matrix: too many elements. Band matrix is overflowing."
else band_matrix := mkmatrix(sq_size,sq_size);
for i:=1:sq_size do
<<
for j:=1:sq_size do
<<
if middle_pos-i+j > 0 and middle_pos-i+j <= no_elts then
setmat(band_matrix,i,j,nth(elt_list,middle_pos-i+j));
>>;
>>;
return band_matrix;
end;
flag('(band_matrix),'opfn);
symbolic procedure make_identity(sq_size);
%
% Creates identity matrix.
%
if not !*fast_la and not fixp sq_size then
rederr "Error in make_identity: non integer input."
else 'mat. (for i:=1:sq_size collect
for j:=1:sq_size collect if i=j then 1 else 0);
flag('(make_identity),'opfn);
symbolic procedure squarep(in_mat);
%
% Tests matrix is square. (boolean).
%
begin
scalar tmp;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in squarep: non matrix input";
tmp := size_of_matrix(in_mat);
if first tmp neq second tmp
then return nil
else return t;
end;
flag('(squarep),'boolean);
flag('(squarep),'opfn);
symbolic procedure swap_rows(in_mat,row1,row2);
%
% Swaps row1 with rows.
%
begin
scalar new_mat;
integer rowdim;
if not !*fast_la then
<<
if not matrixp in_mat then
rederr "Error in swap_rows(first argument): should be a matrix.";
rowdim := row_dim(in_mat);
if not fixp row1 then rederr
"Error in swap_rows(second argument): should be an integer.";
if not fixp row2 then
rederr "Error in swap_rows(third argument): should be an integer.";
if row1>rowdim or row1=0 then rederr
"Error in swap_rows(second argument): out of range for input matrix.";
if row2>rowdim or row2=0 then rederr
"Error in swap_rows(third argument): out of range for input matrix.";
>>;
new_mat := copy_mat(in_mat);
swap_elt(cdr new_mat,row1,row2);
return new_mat;
end;
symbolic procedure swap_columns(in_mat,col1,col2);
%
% Swaps col1 with col2.
%
begin
scalar new_mat;
integer coldim;
if not !*fast_la then
<<
if not matrixp in_mat then rederr
"Error in swap_columns(first argument): should be a matrix.";
coldim := column_dim(in_mat);
if not fixp col1 then rederr
"Error in swap_columns(second argument): should be an integer.";
if not fixp col2 then rederr
"Error in swap_columns(third argument): should be an integer.";
if col1>coldim or col1=0 then rederr
"Error in swap_columns(second argument): out of range for matrix.";
if col2>coldim or col2=0 then rederr
"Error in swap_columns(third argument): out of range for input matrix.";
>>;
new_mat := copy_mat(in_mat);
for each row in cdr new_mat do swap_elt(row,col1,col2);
return new_mat;
end;
symbolic procedure swap_entries(in_mat,entry1,entry2);
%
% Swaps the two entries in in_mat.
%
% entry1 and entry2 must be lists of the form
% {row position,column position}.
%
begin
scalar new_mat;
integer rowdim,coldim;
if not matrixp(in_mat) then rederr
"Error in swap_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 swap_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 swap_entries(third argument): should be a list of 2 elements."
else entry2 := cdr entry2;
if not !*fast_la then
<<
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
if not fixp car entry1 then
<<
prin2 "***** Error in swap_entries(second argument): ";
prin2t " first element in list must be an integer.";
return;
>>;
if not fixp cadr entry1 then
<<
prin2 "***** Error in swap_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 swap_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 swap_entries(second argument): ";
prin2t " second element is out of range for input matrix.";
return;
>>;
if not fixp car entry2 then
<<
prin2 "***** Error in swap_entries(third argument): ";
prin2t " first element in list must be an integer.";
return;
>>;
if not fixp cadr entry2 then
<<
prin2 "***** Error in swap_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 swap_entries(third argument): ";
prin2t " first element is out of range for input matrix.";
return;
>>;
if cadr entry2 > coldim then
<<
prin2 "***** Error in swap_entries(third argument): ";
prin2t " second element is out of range for input matrix.";
return;
>>;
>>;
new_mat := copy_mat(in_mat);
setmat(new_mat,car entry1,cadr entry1,
getmat(in_mat,car entry2,cadr entry2));
setmat(new_mat,car entry2,cadr entry2,
getmat(in_mat,car entry1,cadr entry1));
return new_mat;
end;
% flag('(swap_rows,swap_columns,swap_entries),'opfn);
rtypecar swap_rows,swap_columns,swap_entries;
symbolic procedure get_rows(in_mat,row_list);
%
% 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.
%
begin
integer rowdim,coldim;
scalar ans,tmp;
if not matrixp(in_mat) then
rederr "Error in get_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 get_rows(second argument): ";
prin2t " should be either an integer or a list of integers.";
return;
>>;
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
for each elt in row_list do
<<
if not fixp elt then rederr
"Error in get_rows(second argument): contains non integer.";
if elt>rowdim or elt=0 then
<<
prin2 "***** Error in get_rows(second argument): ";
rederr
"contains row number which is out of range for input matrix.";
>>;
tmp := 'mat.{nth(cdr in_mat,elt)};
ans := append(ans,{tmp});
>>;
return 'list.ans;
end;
symbolic procedure get_columns(in_mat,col_list);
%
% 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.
%
begin
integer rowdim,coldim;
scalar ans,tmp;
if not matrixp(in_mat) then
rederr "Error in get_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 get_columns(second argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
for each elt in col_list do
<<
if not fixp elt then rederr
"Error in get_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.";
>>;
tmp := 'mat.for each row in cdr in_mat collect {nth(row,elt)};
ans := append(ans,{tmp});
>>;
return 'list.ans;
end;
flag('(get_rows,get_columns),'opfn);
symbolic procedure stack_rows(in_mat,row_list);
%
% 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.
%
begin
if not !*fast_la and not matrixp in_mat then
rederr "Error in stack_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;
return 'mat.for each elt in row_list collect nth(cdr in_mat,elt);
end;
symbolic procedure augment_columns(in_mat,col_list);
%
% 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.
%
begin
if not !*fast_la and not matrixp in_mat then rederr
"Error in augment_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;
return 'mat.for each row in cdr in_mat collect
for each elt in col_list collect nth(row,elt);
end;
% flag('(stack_rows,augment_columns),'opfn);
rtypecar stack_rows,augment_columns;
symbolic procedure add_rows(in_mat,r1,r2,mult1);
%
% Replaces row2 (r2) by mult1*r1 + r2.
%
begin
scalar new_mat;
integer i,rowdim,coldim;
coldim := column_dim(in_mat);
if not !*fast_la then
<<
if not matrixp in_mat then
rederr "Error in add_rows(first argument): should be a matrix.";
rowdim := row_dim(in_mat);
if not fixp r1 then
rederr "Error in add_rows(second argument): should be an integer.";
if not fixp r2 then
rederr "Error in add_rows(third argument): should be an integer.";
if r1>rowdim or r1=0 then rederr
"Error in add_rows(second argument): out of range for input matrix.";
if r2>rowdim or r2=0 then rederr
"Error in add_rows(third argument): out of range for input matrix.";
>>;
new_mat := copy_mat(in_mat);
% Efficiency.
if (my_reval mult1) = 0 then return new_mat;
for i:=1:coldim do
setmat(new_mat,r2,i,reval {'plus,{'times,mult1,
getmat(new_mat,r1,i)},getmat(in_mat,r2,i)});
return new_mat;
end;
symbolic procedure add_columns(in_mat,c1,c2,mult1);
%
% Replaces column2 (c2) by mult1*c1 + c2.
%
begin
scalar new_mat;
integer i,rowdim,coldim;
rowdim := row_dim(in_mat);
if not !*fast_la then
<<
if not matrixp in_mat then
rederr "Error in add_columns(first argument): should be a matrix.";
coldim := column_dim(in_mat);
if not fixp c1 then rederr
"Error in add_columns(second argument): should be an integer.";
if not fixp c2 then rederr
"Error in add_columns(third argument): should be an integer.";
if c1>coldim or c1=0 then rederr
"Error in add_columns(second argument): out of range for input matrix.";
if c2>rowdim or c2=0 then rederr
"Error in add_columns(third argument): out of range for input matrix.";
>>;
new_mat := copy_mat(in_mat);
% Why not be efficient.
if (my_reval mult1) = 0 then return new_mat;
for i:=1:rowdim do
setmat(new_mat,i,c2,{'plus,{'times,mult1,getmat(new_mat,i,c1)},
getmat(in_mat,i,c2)});
return new_mat;
end;
% flag('(add_rows,add_columns),'opfn);
rtypecar add_rows,add_columns;
symbolic procedure add_to_rows(in_mat,row_list,value);
%
% Adds value to each element in each row in row_list.
%
% row_list can be either an integer or a list of integers.
%
begin
scalar new_mat;
integer i,rowdim,coldim;
if not matrixp in_mat then
rederr "Error in add_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 add_to_rows(second argument): ";
prin2t " should be either integer or a list of integers.";
return;
>>;
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
new_mat := copy_mat(in_mat);
for each row in row_list do
<<
if not fixp row then rederr
"Error in add_to_row(second argument): should be an integer.";
if row>rowdim or row=0 then
<<
prin2 "***** Error in add_to_rows(second argument): ";
rederr "contains row which is out of range for input matrix.";
>>;
for i:=1:coldim do
setmat(new_mat,row,i,{'plus,getmat(new_mat,row,i),value});
>>;
return new_mat;
end;
symbolic procedure add_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;
integer i,rowdim,coldim;
if not matrixp in_mat then rederr
"Error in add_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 add_to_columns(second argument): ";
prin2t " should be either integer or list of integers.";
return;
>>;
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
new_mat := copy_mat(in_mat);
for each col in col_list do
<<
if not fixp col then rederr
"Error in add_to_columns(second argument): should be an integer.";
if col>coldim or col=0 then
<<
prin2 "***** Error in add_to_columns(second argument): ";
rederr
"contains column which is out of range for input matrix.";
>>;
for i:=1:rowdim do
setmat(new_mat,i,col,{'plus,getmat(new_mat,i,col),value});
>>;
return new_mat;
end;
% flag('(add_to_rows,add_to_columns),'opfn);
rtypecar add_to_rows,add_to_columns;
symbolic procedure mult_rows(in_mat,row_list,mult1);
%
% Replaces rows specified in row_list by row * mult1.
%
begin
scalar new_mat;
integer i,rowdim,coldim;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in mult_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 := row_dim(in_mat);
coldim := column_dim(in_mat);
new_mat := copy_mat(in_mat);
for each row in row_list do
<<
if not !*fast_la and not fixp row then rederr
"Error in mult_rows(second argument): contains non integer.";
if not !*fast_la and (row>rowdim or row=0) then
<<
prin2 "***** Error in mult_rows(second argument): ";
rederr "contains row that is out of range for input matrix.";
>>;
for i:=1:coldim do
<<
setmat(new_mat,row,i,reval {'times,mult1,getmat(in_mat,row,i)});
>>;
>>;
return new_mat;
end;
symbolic procedure mult_columns(in_mat,column_list,mult1);
%
% Replaces columns specified in column_list by column * mult1.
%
begin
scalar new_mat;
integer i,rowdim,coldim;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in mult_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 := row_dim(in_mat);
coldim := column_dim(in_mat);
new_mat := copy_mat(in_mat);
for each column in column_list do
<<
if not !*fast_la and not fixp column then rederr
"Error in mult_columns(second argument): contains non integer.";
if not !*fast_la and (column>coldim or column=0) then
<<
prin2 "***** Error in mult_columns(second argument): ";
rederr "contains column that is out of range for input matrix.";
>>;
for i:=1:rowdim do
<<
setmat(new_mat,i,column,
reval {'times,mult1,getmat(in_mat,i,column)});
>>;
>>;
return new_mat;
end;
% flag('(mult_rows,mult_columns),'opfn);
rtypecar mult_rows,mult_columns;
%%%%%%%%%%%%%%%%%%%%% matrix_augment/matrix_stack %%%%%%%%%%%%%%%%%%%%%%
put('matrix_augment,'psopfn,'matrix_augment1);
symbolic procedure matrix_augment1(matrices);
%
% Takes any number of matrices and joins them horizontally.
%
% Can take either a list of matrices or the matrices as seperate
% arguments.
%
begin
scalar mat_list,new_list,new_row;
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 reval elt;
for each elt in mat_list do
if not matrixp(elt) then
rederr "Error in matrix_augment: non matrix in input.";
>>;
const_rows_test(mat_list);
for i:=1:row_dim(first mat_list) do
<<
new_row := {};
for each mat1 in mat_list do
new_row := append(new_row,nth(cdr mat1,i));
new_list := append(new_list,{new_row});
>>;
return 'mat.new_list;
end;
put('matrix_stack,'psopfn,'matrix_stack1);
symbolic procedure matrix_stack1(matrices);
%
% Takes any number of matrices and joins them vertically.
%
% Can take either a list of matrices or the matrices as seperate
% arguments.
%
begin
scalar mat_list,new_list;
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 reval elt;
for each elt in mat_list do
if not matrixp(elt) then
rederr "Error in matrix_stack: non matrix in input.";
>>;
const_columns_test(mat_list);
for each mat1 in mat_list do new_list := append(new_list,cdr mat1);
return 'mat.new_list;
end;
symbolic procedure no_rows(mat_list);
%
% Takes list of matrices and sums the no. of rows.
%
for each mat1 in mat_list sum row_dim(mat1);
symbolic procedure no_cols(mat_list);
%
% Takes list of matrices and sums the no. of columns.
%
for each mat1 in mat_list sum column_dim(mat1);
symbolic procedure const_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 := row_dim(car mat_list);
i := 1;
while i<listlen and row_dim(car mat_list) = row_dim(cadr mat_list)
do << i := i+1; mat_list := cdr mat_list; >>;
if i=listlen then return rowdim
else
<<
prin2 "***** Error in matrix_augment: ";
rederr "all input matrices must have the same row dimension.";
>>;
end;
symbolic procedure const_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 := column_dim(car mat_list);
i := 1;
while i<listlen and column_dim(car mat_list) =
column_dim(cadr mat_list)
do << i := i+1; mat_list := cdr mat_list; >>;
if i=listlen then return coldim
else
<<
prin2 "***** Error in matrix_stack: ";
rederr
"all input matrices must have the same column dimension.";
return;
>>;
end;
%%%%%%%%%%%%%%%%%%%% end matrix_augment/matrix_stack %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% block_matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symbolic procedure block_matrix(rows,cols,mat_list);
%
% Creates a matrix consisting of rows*cols matrices which are taken
% sequentially from the 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 block_matrix(first argument): ";
prin2t " should be an integer greater than 0.";
return;
>>;
if not fixp cols then rederr
"Error in block_matrix(second argument): should be an integer.";
if cols=0 then
<<
prin2 "***** Error in block_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 block_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 block_matrix(third argument): Incorrect number of matrices.";
row_list := create_row_list(rows,cols,mat_list);
rowdim := check_rows(row_list);
coldim := check_cols(row_list);
block_mat := mkmatrix(rowdim,coldim);
start_row := 1; start_col := 1;
for i:=1:length row_list do
<<
for j:=1:cols do
<<
block_mat := copy_into(nth(nth(row_list,i),j),block_mat,
start_row,start_col);
start_col := start_col + column_dim(nth(nth(row_list,i),j));
>>;
start_col := 1;
start_row := start_row + row_dim(nth(nth(row_list,i),1));
>>;
return block_mat;
end;
flag('(block_matrix),'opfn);
symbolic procedure create_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;
integer i,j,increment;
increment := 1;
for i:=1:rows do
<<
tmp_list := {};
for j:=1:cols do
<<
tmp_list := append(tmp_list,{nth(mat_list,increment)});
increment := increment + 1;
>>;
row_list := append(row_list,{tmp_list});
>>;
return row_list;
end;
symbolic procedure check_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 no_cols(nth(row_list,i)) = no_cols(nth(row_list,i+1))
do i:=i+1;
if i=listlen then return no_cols(nth(row_list,i))
else
<<
prin2
"***** Error in block_matrix: column dimensions of matrices ";
prin2t " into block_matrix are not compatible";
return;
>>;
end;
symbolic procedure check_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 row_dim(nth(nth(row_list,i),j)) =
row_dim(nth(nth(row_list,i),j+1)) then j := j+1
else
<<
prin2 "***** Error in block_matrix: row dimensions of ";
rederr "matrices into block_matrix are not compatible";
>>;
>>;
rowdim := rowdim + row_dim(nth(nth(row_list,i),j));
i := i+1;
>>;
return rowdim;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% end block_matrix %%%%%%%%%%%%%%%%%%%%%%%%%%
put('vandermonde,'psopfn,'vandermonde1);
symbolic procedure vandermonde1(variables);
%
% Input can be either a list or individual arguments.
%
% Creates the Vandermonde matrix.
% ie: the square matrix in which the (i,j)'th entry is
% nth(variables,i)^(j-1).
%
begin
scalar vand,in_list;
integer i,j,sq_size;
if pairp variables and pairp car variables
and caar variables = 'list then variables := cdar variables;
in_list := for each elt in variables collect my_reval elt;
sq_size := length in_list;
vand := mkmatrix(sq_size,sq_size);
for i:=1:sq_size do
<<
for j:=1:sq_size do
<<
setmat(vand,i,j,
reval{'expt,nth(in_list,i),{'plus,j,{'minus,1}}});
>>;
>>;
return vand;
end;
put('toeplitz,'psopfn,'toeplitz1);
symbolic procedure toeplitz1(variables);
%
% Input can be either a list or individual arguments.
%
% Creates the Toeplitz matrix.
% ie: the square matrix in which the first element is placed on the
% diagonal and the nth(variables,i) element is placed on the (i-1)
% sub and super diagonals.
%
begin
scalar toep,in_list;
integer i,j,sq_size;
if pairp variables and pairp car variables
and caar variables = 'list then variables := cdar variables;
in_list := for each elt in variables collect my_reval elt;
sq_size := length in_list;
toep := mkmatrix(sq_size,sq_size);
for i:=1:sq_size do
<<
for j:=0:i-1 do
<<
setmat(toep,i,i-j,nth(in_list,j+1));
setmat(toep,i-j,i,nth(in_list,j+1));
>>;
>>;
return toep;
end;
%%%%%%%%%%%%%%%%%%%%%%%%% kronecker_product %%%%%%%%%%%%%%%%%%%%%%%%%%%%
symbolic procedure kronecker_product(AA,BB);
%
% Copies matrix BB into AA with BB(1,1) at AA(p,q).
%
begin
scalar A,B; integer m,n,r,c;
if not !*fast_la then
<<
if not matrixp(aa) then rederr
"Error in kronecker_product (first argument): should be a matrix.";
if not matrixp(bb) then rederr
"Error in kronecker_product (second argument): should be a matrix.";
>>;
m := row_dim(AA);
n := column_dim(AA);
r := row_dim(BB);
c := column_dim(BB);
A := mkmatrix(m*r,n*c);
for i:=1:m do
for j:=1:n do <<
B := getmat(AA,i,j);
for ii:=1:c do
for jj := 1 : r do
setmat(A,(i-1)*r+jj,(j-1)*c+ii,
reval list('times,b, getmat(bb,jj,ii)));
>>;
return A;
end;
% flag('(kronecker_product),'opfn);
rtypecar kronecker_product;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% minor %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symbolic procedure minor(in_mat,row,col);
%
% Removes row (row) and column (col) from in_mat.
%
begin
scalar min;
if not !*fast_la then
<<
if not matrixp(in_mat) then
rederr "Error in minor(first argument): should be a matrix.";
if not fixp row then
rederr "Error in minor(second argument): should be an integer.";
if row>row_dim(in_mat) or row=0 then rederr
"Error in minor(second argument): out of range for input matrix.";
if not fixp col then
rederr "Error in minor(third argument): should be an integer.";
if col>column_dim(in_mat) or col=0 then rederr
"Error in minor(second argument): out of range for input matrix.";
>>;
min := remove_rows(in_mat,row);
min := remove_columns(min,col);
return min;
end;
symbolic procedure remove_rows(in_mat,row_list);
%
% Removes each row in row_list from in_mat.
%
% row_list can be either an integer or a list of integers.
%
begin
scalar unique_row_list,new_list;
integer rowdim,row;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in remove_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 remove_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 := row_dim(in_mat);
if not !*fast_la then
<<
for each row in unique_row_list do if not fixp row then rederr
"Error in remove_rows(second argument): contains a non integer.";
% rowdim := row_dim(in_mat);
% coldim := column_dim(in_mat);
for each row in unique_row_list do if row>rowdim or row=0 then
rederr
"Error in remove_rows(second argument): out of range for input matrix.";
if length unique_row_list = rowdim then
<<
prin2 "***** Warning in remove_rows:";
prin2t " all the rows have been removed. Returning nil.";
return nil;
>>;
>>;
for row:=1:rowdim do if not intersection({row},unique_row_list)
then new_list := append(new_list,{nth(cdr in_mat,row)});
return 'mat.new_list;
end;
symbolic procedure remove_columns(in_mat,col_list);
%
% Removes each column in col_list from in_mat.
%
% col_list can be either an integer or a list of integers.
%
begin
scalar unique_col_list,new_list,row_list;
integer coldim,row,col;
if not !*fast_la and not matrixp(in_mat) then rederr
"Error in remove_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 remove_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 := column_dim(in_mat);
if not !*fast_la then
<<
for each col in unique_col_list do if not fixp col then rederr
"Error in remove_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 remove_columns(second argument): out of range for matrix.";
if length unique_col_list = coldim then
<<
prin2 "***** Warning in remove_columns: ";
prin2t " all the columns have been removed. Returning nil.";
return nil;
>>;
>>;
for each row in cdr in_mat do
<<
row_list := {};
for col:=1:coldim do
<<
if not intersection({col},unique_col_list)
then row_list := append(row_list,{nth(row,col)});
>> ;
new_list := append(new_list,{row_list});
>>;
return 'mat.new_list;
end;
flag('(minor,remove_rows,remove_columns),'opfn);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end minor %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% begin random_matrix/im_random_matrix %%%%%%%%%%%%%%%%%
switch imaginary; % If ON, then random_matrix creates a random
% matrix with imaginary entries.
switch not_negative; % If ON, then the random matrix functions create
% matrices with only positive entries. In the
% imaginary case, each entry x+iy will have x and
% y both > 0. Not that this really means a great
% deal mathematically apart from each guy sitting
% right up there in the top right hand corner of
% the argand plane but oh well.
switch only_integer; % If ON, then the random matrix functions will
% create matrices with only integer entries. In
% the imaginary case, each entry x+iy will have
% x and y as integers.
switch symmetric; % If ON, random matrix is symmetric.
switch upper_matrix; % If ON, then the random matrix is an upper
% triagonal matrix.
switch lower_matrix; % If ON, then the random matrix is a lower
% triagonal matrix.
symbolic procedure random_minus(limit);
%
% Creates random number in the range -limit < number < limit.
%
begin
scalar r;
r := random(limit);
if evenp random(1000) then r := {'minus,r};
return r
end;
symbolic procedure random_make_minus(u);
%
% Randomly makes u negative.
%
if evenp random(1000) then {'minus,u} else u;
symbolic procedure random_matrix(rowdim,coldim,limit);
%
% Creates an rowdim by coldim matrix with random entries in the bound
% -limit < entry < limit.
%
begin
scalar randmat,random_decimal;
integer i,j,start,current_precision;
if !*lower_matrix and !*upper_matrix then
<<
prin2 "***** Error in random matrix: ";
prin2t " both upper_matrix and lower_matrix switches are on.";
return;
>>;
if !*upper_matrix and !*symmetric then
<<
prin2 "***** Error in random_matriix: ";
prin2t " both upper_matrix and symmetric switches are on.";
return;
>>;
if !*lower_matrix and !*symmetric then
<<
prin2 "***** Error in random_matrix: ";
prin2t " both lower_matrix and symmetric switches are on.";
return;
>>;
if not fixp limit then limit := algebraic floor(abs(limit));
if not fixp rowdim then rederr
"Error in random_matrix(first argument): should be an integer.";
if rowdim=0 then rederr
"Error in random_matrix(first argument): should be integer > than 0.";
if not fixp coldim then rederr
"Error in random_matrix(second argument): should be an integer.";
if coldim=0 then
<<
prin2 "***** Error in random_matrix(second argument): ";
prin2t " should be an integer greater than 0.";
return;
>>;
current_precision := precision 0;
if !*imaginary then randmat := im_random_matrix(rowdim,coldim,limit)
else
<<
start := 1;
randmat := mkmatrix(rowdim,coldim);
for i:=1:rowdim do
<<
if !*symmetric or !*lower_matrix then coldim := i
else if !*upper_matrix then start := i;
for j:=start:coldim do
begin
scalar r1, r2;
r1 := random(limit);
r2 := random(10^current_precision);
random_decimal := {'plus,r1,{'quotient,
r2,
10^current_precision}};
if !*only_integer and !*not_negative then
setmat(randmat,i,j,random(limit))
else if !*only_integer then
setmat(randmat,i,j,random_minus(limit))
else if !*not_negative then
setmat(randmat,i,j,random_decimal)
else setmat(randmat,i,j,random_make_minus(random_decimal));
if !*symmetric then setmat(randmat,j,i,getmat(randmat,i,j));
end;
>>;
>>;
return randmat;
end;
flag('(random_matrix),'opfn);
symbolic procedure im_random_matrix(rowdim,coldim,limit);
%
% Creates an rowdim by coldim matrix with random imaginary entries.
% The entrirs are of the form x+iy where x and y are in the bound
% -limit < x,y < limit.
%
begin
scalar randmat,random_decimal,im_random_decimal;
integer i,j,start,current_precision;
start := 1;
current_precision := precision 0;
randmat := mkmatrix(rowdim,coldim);
for i:=1:rowdim do
<<
if !*symmetric or !*lower_matrix then coldim := i
else if !*upper_matrix then start := i;
for j:=start:coldim do
begin
scalar r1, r2;
r1 := random(limit);
r2 := random(10^current_precision);
random_decimal := {'plus,1,{'quotient,
r2,
10^current_precision}};
r1 := random(limit);
r2 := random(10^current_precision);
im_random_decimal := {'plus,r1,{'quotient,
r2,
10^current_precision}};
if !*only_integer and !*not_negative then <<
r1 := random(limit); r2 := random(limit);
setmat(randmat,i,j,{'plus,r1, {'times,'i,r2}}) >>
else if !*only_integer then <<
r1 := random_minus(limit); r2 := random_minus(limit);
setmat(randmat,i,j,{'plus,r1, {'times,'i,r2}}) >>
else if !*not_negative then
setmat(randmat,i,j,{'plus,random_decimal,
{'times,'i,im_random_decimal}})
else <<
r1 := random_make_minus(random_decimal);
r2 := random_make_minus(im_random_decimal);
setmat(randmat,i,j,{'plus,r1, {'times,'i,r2}}) >>;
if !*symmetric then setmat(randmat,j,i,getmat(randmat,i,j));
end;
>>;
return randmat;
end;
% flag('(im_random_matrix),'opfn);
rtypecar im_random_matrix;
%%%%%%%%%%%%%%%%%% end random_matrix/im_random_matrix %%%%%%%%%%%%%%%%%%
symbolic procedure extend(in_mat,rows,cols,entry);
%
% Extends in_mat by rows rows (!) and cols columns. New entries are
% initialised to entry.
%
begin
scalar ex_mat;
integer rowdim,coldim,i,j;
if not matrixp(in_mat) then
rederr "Error in extend(first argument): should be a matrix.";
if not fixp rows then
rederr "Error in extend(second argument): should be an integer.";
if not fixp cols then
rederr "Error in extend(third argument): should be an integer.";
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
ex_mat := mkmatrix(rowdim+rows,coldim+cols);
ex_mat := copy_into(in_mat,ex_mat,1,1);
for i:=1:rowdim+rows do
<<
for j:=1:coldim+cols do
<<
if i<=rowdim and j<=coldim then <<>>
else setmat(ex_mat,i,j,entry);
>>;
>>;
return ex_mat;
end;
flag('(extend),'opfn);
rtypecar extend;
%%%%%%%%%%%%%%%%%%%%% begin char_matrix/char_poly %%%%%%%%%%%%%%%%%%%%%%
symbolic procedure char_matrix(in_mat,lmbda);
%
% Create characteristic matrix. ie: C := lmbda*I - in_mat.
% in_ mat must be square.
%
begin
scalar carmat;
integer rowdim;
if not matrixp(in_mat) then
rederr "Error in char_matrix(first argument): should be a matrix.";
if not squarep(in_mat) then rederr
"Error in char_matrix(first argument): must be a square matrix.";
rowdim := row_dim(in_mat);
carmat := {'plus,{'times,lmbda,make_identity(rowdim)},
{'minus,in_mat}};
return carmat;
end;
symbolic procedure char_poly(in_mat,lmbda);
%
% Finds characteristic polynomial of matrix in_mat.
% ie: det(lmbda*I - in_mat).
%
begin
scalar chpoly,carmat;
if not matrixp(in_mat) then
rederr "Error in char_poly(first argument): should be a matrix.";
carmat := char_matrix(in_mat,lmbda);
chpoly := algebraic det(carmat);
return chpoly;
end;
flag('(char_matrix char_poly),'opfn);
%%%%%%%%%%%%%%%%%%%%%%%% end char_matrix/char_poly %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% begin pivot/rows_pivot %%%%%%%%%%%%%%%%%%%%%%
symbolic procedure pivot(in_mat,pivot_row,pivot_col);
%
% Converts all elements in pivot column (apart from the one in pivot
% row) to 0.
%
begin
scalar piv_mat,ratio;
integer i,j,rowdim,coldim;
if not !*fast_la and not matrixp(in_mat) then
rederr "Error in pivot(first argument): should be a matrix.";
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
if not !*fast_la then
<<
if not fixp pivot_row then
rederr "Error in pivot(second argument): should be an integer.";
if pivot_row>rowdim or pivot_row=0 then rederr
"Error in pivot(second argument): out of range for input matrix.";
if not fixp pivot_col then rederr
"Error in pivot(third argument): should be an integer.";
if pivot_col>coldim or pivot_col=0 then rederr
"Error in pivot(third argument): out of range for input matrix.";
if getmat(in_mat,pivot_row,pivot_col) = 0 then
rederr "Error in pivot: cannot pivot on a zero entry.";
>>;
piv_mat := copy_mat(in_mat);
piv_mat := copy_mat(in_mat);
for i:=1:rowdim do
<<
for j:=1:coldim do
<<
if i = pivot_row then <<>>
else
<<
ratio := {'quotient,getmat(in_mat,i,pivot_col),
getmat(in_mat,pivot_row,pivot_col)};
setmat(piv_mat,i,j,{'plus,getmat(in_mat,i,j),{'minus,
{'times,ratio,getmat(in_mat,pivot_row,j)}}});
>>;
>>;
>>;
return piv_mat;
end;
symbolic procedure rows_pivot(in_mat,pivot_row,pivot_col,row_list);
%
% 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.
%
begin
scalar piv_mat,ratio;
integer j,rowdim,coldim;
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
if not !*fast_la then
<<
if not matrixp(in_mat) then
rederr "Error in rows_pivot(first argument): should be a matrix.";
rowdim := row_dim(in_mat);
coldim := column_dim(in_mat);
if not fixp pivot_row then
rederr "Error in pivot(second argument): should be an integer.";
if pivot_row>rowdim or pivot_row=0 then rederr
"Error in rows_pivot(second argument): out of range for input matrix.";
if not fixp pivot_col then
rederr "Error in pivot(third argument): should be an integer.";
if pivot_col>coldim or pivot_col=0 then rederr
"Error in rows_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 rows_pivot(fourth argument): ";
prin2t
" should be either an integer or a list of integers.";
return;
>>;
if getmat(in_mat,pivot_row,pivot_col) = 0 then
rederr "Error in rows_pivot: cannot pivot on a zero entry.";
piv_mat := copy_mat(in_mat);
for each elt in row_list do
<<
if not !*fast_la then
<<
if not fixp elt then rederr
"Error in rows_pivot: fourth argument contains a non integer.";
if elt>rowdim or elt=0 then
<<
prin2 "***** Error in rows_pivot(fourth argument): ";
rederr "contains row which is out of range for input matrix.";
>>;
>>;
for j:=1:coldim do
<<
if elt = pivot_row then <<>>
else
<<
ratio := {'quotient,getmat(in_mat,elt,pivot_col),
getmat(in_mat,pivot_row,pivot_col)};
setmat(piv_mat,elt,j,{'plus,getmat(in_mat,elt,j),{'minus,
{'times,ratio,getmat(in_mat,pivot_row,j)}}});
>>;
>>;
>>;
return piv_mat;
end;
flag('(pivot,rows_pivot),'opfn);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end pivot/rows_pivot %%%%%%%%%%%%%%%%%%%%%
symbolic procedure jacobian(exp_list,var_list);
%
% 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)).
%
begin
scalar jac,exp1,var1;
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 jacobian(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 := mkmatrix(rowdim,coldim);
for i:=1:rowdim do
<<
for j:=1:coldim do
<<
exp1 := nth(exp_list,i);
var1 := nth(var_list,j);
setmat(jac,i,j,algebraic df(exp1,var1));
>>;
>>;
return jac;
end;
flag('(jacobian),'opfn);
symbolic procedure hessian(poly,variables);
%
% 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.
%
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 hessian(second argument): ";
prin2t
" should be either a single variable or a list of variables.";
return;
>>;
sq_size := length variables;
hess_mat := mkmatrix(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);
setmat(hess_mat,row,col,elt);
>>;
>>;
return hess_mat;
end;
flag('(hessian),'opfn);
symbolic procedure hermitian_tp(in_mat);
%
% 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.
%
begin
scalar h_tp,element;
integer row,col;
if not matrixp(in_mat) then
rederr "Error in hermitian_tp: non matrix input.";
h_tp := algebraic tp(in_mat);
for row:=1:row_dim(h_tp) do
<<
for col:=1:column_dim(h_tp) do
<<
element := getmat(h_tp,row,col);
setmat(h_tp,row,col,
algebraic (repart(element) - i*impart(element)));
>>;
>>;
return h_tp;
end;
flag('(hermitian_tp),'opfn);
symbolic procedure hilbert(sq_size,value);
%
% The Hilbert matrix is symmetric and the (i,j)'th entry in
% 1/(i+j-x).
%
begin
scalar hil_mat,denom;
integer row,col;
if not fixp sq_size or sq_size<1 then rederr
"Error in hilbert(first argument): must be a positive integer.";
hil_mat := mkmatrix(sq_size,sq_size);
for row:=1:sq_size do
<<
for col:=1:sq_size do
<<
if (denom := reval{'plus,row,col,{'minus,value}}) = 0 then
rederr "Error in hilbert: division by zero."
else setmat(hil_mat,row,col,{'quotient,1,denom});
>>;
>>;
return hil_mat;
end;
flag('(hilbert),'opfn);
%%%%%%%%%%%%%%%%%%%%%%%% begin coeff_matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%
put('coeff_matrix,'psopfn,'coeff_matrix1); % To allow variable input.
symbolic procedure coeff_matrix1(equation_list);
%
% 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.
%
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 coeff_matrix: no variables in input.";
check_linearity(equation_list,variable_list);
A := get_A(equation_list,variable_list);
X := get_X(variable_list);
b := get_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 coeff_matrix: the equations are not linear.";
>>;
>>;
symbolic procedure get_A(equation_list,variable_list);
begin
scalar A,element,var_elt;
integer row,col,length_equation_list,length_variable_list;
length_equation_list := length equation_list;
length_variable_list := length variable_list;
A := mkmatrix(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);
setmat(A,row,col,algebraic coeffn(element,var_elt,1));
>>;
>>;
return A;
end;
symbolic procedure get_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 := mkmatrix(length_integer_list,1);
for row:=1:length_integer_list do
setmat(b,row,1,-nth(integer_list,row));
return b;
end;
symbolic procedure get_X(variable_list);
begin
scalar X;
integer row,length_variable_list;
length_variable_list := length variable_list;
X := mkmatrix(length_variable_list,1);
for row := 1:length variable_list do
setmat(X,row,1,nth(variable_list,row));
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;
%%%%%%%%%%%%%%%%%%%%%%%%%% end coeff_matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smacro used in other modules.
symbolic smacro procedure my_revlis(u);
%
% As my_reval but for lists.
%
for each j in u collect my_reval(j);
endmodule; %linear algebra.
end;