sph@doc.ic.ac.uk (Simon Hughes) (12/04/90)
Can anybody email me (a reference for) functional programs for matrix manipulation? E.g. transposition, multiplication etc. Preferably in a form in which the matrices are represented by lists of lists. My research involves static analysis of functional programs to improve their store use efficiency. I suspect a matrix manipulation package would be a good application of the analyses. Thanks. -- ------------------------------------------------------------------------------ Simon Hughes, Dept. of Computing, Imperial College. Email: sph@doc.ic.ac.uk
sph@iti.org (Simon Hughes) (01/11/91)
Some time ago I posted to this newsgroup, asking people to
email me matrix manipulation programs. Here's what I received.
Thanks to everyone who responded (including those not detailed
below).
Simon
P.S. Larry's programs will appear in the following book.
Lawrence C Paulson.
ML for the Working Programmer.
(Cambridge University Press, 1991, in press).
------------------------------------------------------------------
From jeroen%ruuinf.cs.ruu.nl@nsf.ac.uk Wed Dec 5 09:19:17 1990
This is my version, written in Miranda.
It stores matrices as lists of lists.
It avoids the list-indexing operator (!); instead
it uses functions like map, inits and tails as much as possible.
This satisfies my sense of "functional beauty", although
I must admit that functions like inv may be difficult to
read because of it.
--
Jeroen Fokker | jeroen@cs.ruu.nl
dept.of Computer Science, Utrecht University | tel.+31-30-534129
PObox 80089, 3508TB Utrecht, the Netherlands | fax.+31-30-513791
--
scalar == num
vector == [scalar]
matrix == [vector]
inprod :: vector -> vector -> scalar
apply :: matrix -> vector -> vector
mul :: matrix -> matrix -> matrix
diag :: matrix -> vector
transp :: matrix -> matrix
det :: matrix -> scalar
inv :: matrix -> matrix
solve :: matrix -> vector -> vector
inprod a b = sum (map2 (*) a b)
apply m v = map (inprod v) m
mul as bs = map f as
where f a = map (inprod a)(transp bs)
diag [] = []
diag (l:ls) = hd l : diag (map tl ls)
antidiag = diag . reverse
transp [xs] = map single xs
transp (xs:xss) = map2 (:) xs (transp xss)
gaps xs = map2 (++) (init(inits xs)) (tl(tails xs))
altern = inprod plusmin
det [[x]] = x
det (v:m) = altern (map2 (*) v (map det (gaps (transp m))))
inv [[x]] = [[1/x]]
inv m = transp (zipp (*) c (mapp det (map gaps (map transp (gaps m)))))
where c = mapp (/(det m)) (iterate tl plusmin)
solve m v = map ((/(det m)).det) (repls v (transp m))
where repls v m = map2 (insert v) (init(inits m)) (tl(tails m))
insert v a b = a++ [v] ++b
m = [[1,3,1,1],[2,1,5,2],[1,-1,2,3],[4,1,-3,7]]
n= [[2,1],[4,3]]
mapp f = map (map f)
zipp f = map2 (map2 f)
inits [] = [[]]
inits (x:xs) = [] : map (x:) (inits xs)
tails [] = [[]]
tails (x:xs) = (x:xs) : (tails xs)
plusmin = 1 : -1 : plusmin
single x = [x]
negvec = map neg
a $vecminus b = map2 (-) a b
From Larry.Paulson%uk.ac.cam.cl@cl.cam.ac.uk Wed Dec 5 10:37:03 1990
I enclose ML code for matrix transpose, multiplication, and Gaussian
Elimination. It comes from my forthcoming book on Standard ML. (The code is
entirely first-order because higher-order functions come later in the book!)
Please report any errors you find; other comments are also welcome. Hope you
find the code useful!
Larry Paulson
------
(** Matrix transpose **)
fun headcol [] = []
| headcol ((x::_) :: rows) = x :: headcol rows;
fun tailcols [] = []
| tailcols ((_::xs) :: rows) = xs :: tailcols rows;
val matrix =
[ ["a","b","c"],
["d","e","f"] ];
headcol matrix;
tailcols matrix;
fun transp ([]::rows) = []
| transp rows = headcol rows :: transp (tailcols rows);
transp matrix;
(** Matrix product **)
fun dotprod([], []) = 0.0
| dotprod(x::xs,y::ys) = x*y + dotprod(xs,ys);
fun rowprod(row, []) = []
| rowprod(row, col::cols) =
dotprod(row,col) :: rowprod(row,cols);
fun rowlistprod([], cols) = []
| rowlistprod(row::rows, cols) =
rowprod(row,cols) :: rowlistprod(rows,cols);
fun matprod(Arows,Brows) = rowlistprod(Arows, transp Brows);
val Arows = [ [2.0, 0.0],
[3.0, ~1.0],
[0.0, 1.0],
[1.0, 1.0] ];
val Brows = [ [1.0, 0.0, 2.0],
[4.0, ~1.0, 0.0] ];
matprod(Arows,Brows);
(** Gaussian Elimination (Sedgewick Chapter 5) **)
(*the first row with absolutely greatest head*)
fun pivotrow [row] = row : real list
| pivotrow (row1::row2::rows) =
if abs(hd row1) >= abs(hd row2)
then pivotrow(row1::rows)
else pivotrow(row2::rows);
(*the matrix excluding the first row with head p*)
fun delrow (p, []) = []
| delrow (p, row::rows) =
if p = hd row then rows
else row :: delrow(p, rows);
fun scalarprod(k, []) = [] : real list
| scalarprod(k, x::xs) = k*x :: scalarprod(k,xs);
fun vectorsum ([], []) = [] : real list
| vectorsum (x::xs,y::ys) = x+y :: vectorsum(xs,ys);
fun gausselim [row] = [row]
| gausselim rows =
let val p::prow = pivotrow rows
fun elimcol [] = []
| elimcol ((x::xs)::rows) =
vectorsum(xs, scalarprod(~x/p, prow))
:: elimcol rows
in (p::prow) :: gausselim(elimcol(delrow(p,rows)))
end;
fun solutions [] = [~1.0]
| solutions((x::xs)::rows) =
let val solns = solutions rows
in ~(dotprod(solns,xs)/x) :: solns end;
(*Sedgewick's example*)
gausselim [[ 1.0, 3.0, ~4.0, 8.0],
[ 0.0, 2.0, ~2.0, 6.0],
[~1.0, ~2.0, 5.0, ~1.0]];
solutions it;
(*should be [1,5,2] with ~1 at end*)
gausselim [[ 0.0, 1.0, 2.0, 7.0, 7.0],
[~4.0, 0.0, 0.0, ~5.0, ~17.0],
[ 6.0, ~1.0, 0.0, 0.0, 28.0],
[~2.0, 0.0, 2.0, 8.0, 12.0]];
solutions it;
(*should be [3,~10,5,1] with ~1 at end*)
------
From ahsnippe%praxis.cs.ruu.nl@nsf.ac.uk Fri Dec 7 08:20:10 1990
I've done some programming in Miranda which resulted in a library
for wire-frame-graphics. Part of it is a matrix-manipulation program.
A better version is owned by jeroen@cs.ruu.nl (Jeroen Fokker).
:
vector == [num]
matrix == [vector]
dotprodukt [] [] = 0
dotprodukt (a:l) (b:m) = a*b + dotprodukt l m
apply xs y = map (dotprodukt y) xs
ident n = map (line n) [0..(n-1)]
where line n i = (rep i 0)++[1]++(rep (n-i-1) 0)
scalevec n x = map (*n) x
scalemat n = map (scalevec n)
mul x y = transpose (map (apply x) (transpose y))
negvec = scalevec (-1)
vecminus x y = map2 (-) x y
++++++++++++++++++++++++++
This is the one used in the graphics package
Next has some additonal functions :
++++++++++++++++++++++++++++++++
vector == [num]
matrix == [vector]
inprodukt :: vector -> vector -> num
inprodukt [] [] = 0
inprodukt (a:l) (b:m) = a*b + inprodukt l m
matapply :: matrix -> vector -> vector
matapply xs y = map (inprodukt y) xs
transp :: matrix -> matrix
transp m = (kop m) : (transp (rest m)) , if (vol m)
= [] , otherwise
vol [] = False
vol (x:xs) = (x~=[])
kop = map (hd)
rest = map (tl)
ident :: num -> matrix
ident n = map (line n) [0..(n-1)]
line n i = (rep i 0)++[1]++(rep (n-i-1) 0)
scalvec :: num -> vector -> vector
scalvec n x = map (*n) x
scalmat :: num -> matrix -> matrix
scalmat n = map (scalvec n)
matmul :: matrix -> matrix -> matrix
matmul x = map (matapply x)
det :: matrix -> num
det m = 0 , if ~(vol m)
= hd (hd m) , if #m=1
= som (#m) m , otherwise
som i m = (-1)^i * (((m)!0)!(i-1)) * (det (slicedmat i m)) + (som (i-1) m) , if i >0
= 0 , otherwise
slicedmat n [] = []
slicedmat n (m:ms) = plak (stuk n ms) (rest2 n ms)
stuk n [] = []
stuk n (m:ms) = (take (n-1) m ) : stuk (n-1) ms
rest2 n [] = []
rest2 n (m:ms) = (drop n m) : rest2 n ms
plak [] [] = []
plak (m:ms) (n:ns) = (m ++ n) : plak ms ns
I hope the program is bug-free since I haven't used it for a long time.
If you have any questions, feel free to ask.
--[Arnold]--
From S.Clayman@cs.ucl.ac.uk Fri Dec 7 18:57:24 1990
I suspect you've had a bundle of replies. Anyway here's
a bit of Miranda I wrote ages ago
______________________________________________________________
|| some matrix things
%export scalar_mult_vector scalar_mult_matrix multmatrix
scalar == num
vector == [num]
matrix == [[num]]
|| this does s.v => [s.v1,s.v2 ...]
scalar_mult_vector::scalar->vector->vector
scalar_mult_vector s v = map ((*)s) v
|| this does s.M => [[s.m11,s.m12...],[s.m21,s.m22...]...]
scalar_mult_matrix::scalar->matrix->matrix
scalar_mult_matrix s m = map (scalar_mult_vector s) m
|| this is used in multmatrix to multiply
|| each row of matrix1 by matrix2
vector_mult_matrix::vector->matrix->vector
vector_mult_matrix v m
= domult v m, #v=#m
= error "vector_mult_matrix: vector and matrix have different lengths", otherwise
where
domult v ([]:any) = []
domult v m = sum (multi (map (*) v) (map hd m)) : domult v (map tl m)
multmatrix::matrix->matrix->matrix
multmatrix m1 m2 = multi (map vector_mult_matrix m1) (repeat m2)
|| multi - applies a list of functions to
|| a list of args
|| multi [(*)1.5,(*)2,(+)3] [2,4,6]
|| should return [3,8,9]
multi::[*->**]->[*]->[**]
multi = map2 apply
where
apply f x = f x
______________________________________________________________
I hope it's not too obscure
Stuart
From pierre%ra.mcrcim.mcgill.edu@nsf.ac.uk Sat Dec 8 00:48:38 1990
I have done some work using MIT's eager functional language Id, to
perform LU decomposition. My code is pretty simple-minded, and thus
inefficient, but code given by Arvind and Nikhil in a draft of a
forthcoming book on Id is far better written. Perhaps they could be
convinced to mail you a copy. One of their e-mail address is:
bug-id-world@mc.mit.edu
You can also look in J. Parallel and Distrib. Computing, vol 5 (no 5),
p. 460--493, for some of their earlier work.
Pierre
--
------------------------------------------------------------------------------
Simon Hughes, Dept. of Computing, Imperial College. Email: sph@doc.ic.ac.uk