Hi Jan. I recently spent some time researching precisely this topic. Unfortunately I don't have the exact references at hand but two key starting points are TR433 by David Wise et al and Chris Angus "Numerical Software DEvelopment with Functional Languages". The TR433 report has code written in Gofer which I just yesterday finished converting to Haskell 98. (That is, I compiled it but I have not yet tested it.) It is attached as a starting point. Cheers Mike Thomas. ----- Original Message ----- From: "Jan Kybic" <kybic@ieee.org> To: <haskell@haskell.org> Sent: Tuesday, April 23, 2002 9:33 PM Subject: Matrix library in Haskell
Hello,
I am just discovering Haskell, so sorry if this is not the right place to ask. I want to use it for some numerical calculations. I need something higher level than C++ and faster than Python or Matlab and from the initial experiments it seems that Haskell could be the right choice. My question is: Is there any matrix/vector library available with common operations (dot products, matrix products, linear system solution etc.) ? I could not find any.
I am including her my first try on such a library, in case it might be useful for somebody. However, I am not perfectly happy with the design. In particular I would like to define MatrixClass and VectorClass so that applying a getRow operation on a matrix instance would yield automatically the correct instance of the VectorClass, but I do not know how to express this interdependency, so for the moment I have dropped the classes.
Yours,
Jan
-- ------------------------------------------------------------------------- Jan Kybic <kybic@ieee.org> Robotvis, INRIA, Sophia-Antipolis, France or <Jan.Kybic@sophia.inria.fr>,tel. work +33 492 38 7589, fax 7845 http://www-sop.inria.fr/robotvis/personnel/Jan.Kybic/
-- Module implementing vectors, matrices, and operations on them -- it uses the multiparameter class extension -- -- $Id: LinearAlgebra.hs,v 1.4 2002/04/22 11:44:41 jkybic Exp $ -- Jan Kybic, April 2002
module LinearAlgebra where
import Array import Complex import Observe
-------------------------------------------------------------------- -- now define the concrete vector and array implementations -- TODO: This could be accelerated using IArray/ UArray -- TODO: Make it a proper instance of Show
newtype ArrayVectorType a = ArrayVector (Array Int a) deriving Show newtype ArrayMatrixType a = ArrayMatrix (Array (Int,Int) a) deriving Show
listToVec l = ArrayVector (array (0,uplim) $ zip [0..uplim] l) where uplim = (length l) -1
vecToList (ArrayVector v) = [ v!i | i <- indices v ]
dot (ArrayVector a) (ArrayVector b) = sum [ a!i * b!i | i <- indices a ]
norm2 a = dot a a
norm a = sqrt (norm2 a)
(!@) (ArrayVector v) i = v!i sizeV (ArrayVector v) = rangeSize $ bounds v indicesV v= [0..(sizeV v - 1)]
scaleV (ArrayVector a) s = ArrayVector ( fmap (\x -> s * x) a )
(+@) a b = listToVec [ a!@i + b!@i | i <- indicesV a ]
(-@) a b = listToVec [ a!@i - b!@i | i <- indicesV a ]
-- matrix operations
(!@@) (ArrayMatrix m) (i,j) = m!(i,j)
listToMat l = let { bnds=((0,0),(m,n)) ; m=length l -1 ; n= length (head l) -1} in ArrayMatrix ( array bnds [ ((i,j),(l!!i)!!j) | (i,j) <- range bnds ] ) matToList a = [ [ a !@@ (i,j) | j <- range $ cbounds a ] | i <- range $ rbounds a ]
boundsM (ArrayMatrix a) = bounds a rbounds m = let z=boundsM m in (fst(fst z),fst(snd z)) cbounds m = let z=boundsM m in (snd(fst z),snd(snd z)) nrows = rangeSize . rbounds ncols = rangeSize . cbounds getRow a i = -- extract row i as vector listToVec [ a!@@(i,j) | j <- range $ cbounds a ] getCol a j = -- extract column j as vector listToVec [ a!@@(i,j) | i <- range $ cbounds a ] showMat a = "[ " ++ concat [ showVec (getRow a i) ++ " " | i <- range (rbounds a) ] ++ "]" scaleMat (ArrayMatrix a) s = ArrayMatrix ( fmap (\x -> s * x) a ) matMult a b = let bnds = transTup (rbounds a,cbounds b) in ArrayMatrix (array bnds [ ((i,j), (getRow a i) `dot` getCol b j) | (i,j) <- range bnds ]) idMat n = -- create an identity matrix let bnds=((0,0),(n-1,n-1)) in ArrayMatrix (accumArray (+) 0.0 bnds [ ((i,i),1.0) | i <- [0..n-1] ])
-- cross takes two vectors of length three and calculates their cross product -- uncomment the run-time checks if you prefer safety over speed -- the signature is a little limiting to avoid uncertainity of the resulting -- vector
cross a b -- | or [ sizeV a /=3 , sizeV b /=3 ] = -- error "Cross product needs length 3 vectors" -- | otherwise = listToVec [ a !@ 1 * b !@ 2 - a !@ 2 * b !@ 1, - a !@ 0 * b !@ 2 + a !@ 2 * b !@ 0, a !@ 0 * b !@ 1 - a !@ 1 * b !@ 0 ]
vangle a b = acos (cosvangle a b) -- angle between two vectors
cosvangle' a b = -- the cos of an angle between two vectors let mag= sqrt ( norm2 a * norm2 b ) in if mag>0.0 then (dot a b)/mag else 0.0
-- the above version had numerical problems cosvangle a b = max ( min (cosvangle' a b) 1.0 ) (-1.0)
showVec a = -- converts a vector to a string "[ " ++ concat [ show (a !@ i) ++ " " | i <- [0..sizeV a-1] ] ++ "]"
subvector a (i,j) = -- another vector with indices i..j listToVec [ a !@ k | k <- [i..j] ]
transTup ((a,b),(c,d)) = ((a,c),(b,d)) -- transpose a tuple of tuples
-- type synonymes
type DoubleVector = ArrayVectorType Double type DoubleMatrix = ArrayMatrixType Double
-- Observable instance
instance (Observable a) => Observable (ArrayVectorType a) where observer (ArrayVector a) = send "ArrayVector" (return ArrayVector << a)
instance (Observable a) => Observable (ArrayMatrixType a) where observer (ArrayMatrix a) = send "ArrayMatrix" (return ArrayMatrix << a)
-- Functor instance
instance Functor ArrayVectorType where fmap f (ArrayVector a) = ArrayVector ( fmap f a )
-- end of Linear Algebra.hs _______________________________________________ Haskell mailing list Haskell@haskell.org http://www.haskell.org/mailman/listinfo/haskell