
Martin Percossi wrote:
Hello, I am trying to write a haskell-native matrix library (providing a similar set of functionality as gsl's matrix/eigenvector routines). I have had a look at other attempts at matrix libraries, and have found Hal Daume's and Alberto Ruiz's libraries to offer a good amount of functionality. There are two major reasons for wanting to write my own library:
1. Haskell-nativeness: I have had some issues compiling and linking with gsl libraries on 64-bit platforms. Also, it would be quite interesting to gauge haskell's effectiveness as a scientific computing platform, in particular making use unboxed arrays as representation.
2. Use of monads: in the afore-mentioned libraries, monads are either ignored, or are restricted to the IO monad. I am taking a different approach, starting my code from the ST monad, and eventually this will be generalized to work with the IO monad. I think the ST monad is a good monad to be able to perform computations on matrices that update them, allowing efficient, in-place algorithms, but it also provides the benifit of not being a one-way street like the IO monad is. Being a relative newcomer to haskell, I would be interested to hear any thoughts as to whether this is a good/bad idea.
Now to my question: I would like to represent a matrix as a wrapper around a block, which in turn is just an unboxed array. Here are the definitions for a matrix in ST and outside of a monad, respectively:
type MBlock s = STUArray s Int Double data MMatrix s = MMatrix Int Int (MBlock s) type Block = UArray Int Double data Matrix = Matrix Int Int Block
Now, I have started by providing some fairly low-level routines in the ST monad. My problem is depressingly simple: I would like to retrieve a matrix from the ST monad for use outside of it (in particular, to pretty-print it). Now, this is easy for STUArrays (just use runSTUArray), but I'm not sure of how to do it for a type that *encloses* an STUArray.
For example, here's a simple test:
runMatrix = do _A <- newListMatrix [[2, 0, 0], [0, 2, 0], [0, 0, 2]] _B <- newListMatrix [[1, 2, 1], [0, 1, 1], [0, 0, 1]] _C <- matMul _A _B return $ getBlock _C
main = show $ runSTUArray runMatrix
Obviously, what I get as a result of runSTUArray is a UArray, which is a pain because I'd then have to box it "by hand" (i.e. again specifying the number of rows and columns), instead of automatically, in order to have a Matrix again.
So what I'd like is a runSTMatrix routine, which would possibly have signature:
runSTMatrix :: ST s (Matrix s) -> Matrix
which would have semantics analogous to runSTUArray.
Does anyone have any idea of how I can write runSTMatrix?
Many thanks in advance, Martin
I have not used these, but I know where they are: You want Data.Array.MArray.unsafeFreeze as documented at http://www.haskell.org/ghc/docs/latest/html/libraries/base/Data-Array-MArray... "Converts an mutable array into an immutable array. The implementation may either simply cast the array from one type to the other without copying the array, or it may take a full copy of the array." This can also be used with "thaw" to make a safe copy. Untested: unsafeFreezeMatrix :: MMatrix s -> (ST s) Matrix unsafeFreezeMatrix (MMatrix x1 x2 marray) = Matrix x1 x2 (unsafeFreeze marray) thawMatrix :: Matrix -> (ST s) (MMatrix s) thawMatrix (Matrix x1 x2 iarray) = MMatrix x1 x2 (thaw iarray) addMatrix :: Matrix -> Matrix -> Matrix addMatrix a b = runST (do a' <- thawMatrix a -- This copies a ::Matrix into a' ::MMatrix -- Add to mutable a' the elements of immutable b return (unsafeFreezeMatrix a') ) Using unsafeFreeze at the end of runST is, by itself, safe. Using unsafeThaw looks much more dangerous. -- Chris