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" <[EMAIL PROTECTED]> To: <[EMAIL PROTECTED]> 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 <[EMAIL PROTECTED]> Robotvis, INRIA, Sophia-Antipolis, France > or <[EMAIL PROTECTED]>,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 > [EMAIL PROTECTED] > http://www.haskell.org/mailman/listinfo/haskell
tr433.zip
Description: Zip compressed data
