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

Hello Martin, Sunday, March 12, 2006, 8:49:15 PM, you wrote: MP> 1. Haskell-nativeness: I have had some issues compiling and linking with gsl MP> libraries on 64-bit platforms. Also, it would be quite interesting to gauge MP> haskell's effectiveness as a scientific computing platform, in particular MP> making use unboxed arrays as representation. it's a bad idea. simple test of vector addition shows the 20x performance loss comparing GHC and GCC. so you can consider this as learning example but not as a real-world lib. btw, i've proposed making changes in ghc that will change this situation, at least for simple loops. but that is only words for this moment. if you want to allow ghc became the really good choice for writing matrix libraries, it's better to take part in ghc's own development :) MP> 2. Use of monads: in the afore-mentioned libraries, monads are either ignored, MP> or are restricted to the IO monad. IO monad is just ST monad specialized to RealWorld state. you can easily use this libs in ST monad with help of unsafeIOToST. surprised? :) you will laugh even more if i say that STUArray in Hugs implemented in just this way - by using peek/poke IO operations wrapped in unsafeIOToST. just small excerpt from this implementation: specialIOToST :: IO a -> ST s a specialIOToST = unsafeCoerce type BytePtr = ForeignPtr Word8 data MutableByteArray s = MutableByteArray !Int !BytePtr newMutableByteArray :: Int -> ST s (MutableByteArray s) newMutableByteArray size = do fp <- specialIOToST (mallocForeignPtrBytes size) return (MutableByteArray size fp) readMutableByteArray :: Storable e => MutableByteArray s -> Int -> ST s e readMutableByteArray (MutableByteArray _ fp) i = specialIOToST $ withForeignPtr fp $ \a -> peekElemOff (castPtr a) i MP> I am taking a different approach, starting MP> my code from the ST monad, and eventually this will be generalized to work with MP> the IO monad. I think the ST monad is a good monad to be able to perform MP> computations on matrices that update them, allowing efficient, in-place MP> algorithms, but it also provides the benifit of not being a one-way street like MP> the IO monad is. Being a relative newcomer to haskell, I would be interested to MP> hear any thoughts as to whether this is a good/bad idea. it's a right way (and even obvious way to one who knows how this all work). btw, i developed general i/o and serialization library that works both in IO and ST monad, and now i almost finished rewriting of arrays and references library, which internally uses the same concept of providing general code that works in both monads. just a small excerpt from my code: -- | Unboxed mutable arrays data UnboxedMutableArray s i e = UMA !i !i !(MUVec s e) instance (STorIO m s) => HasMutableBounds (UnboxedMutableArray s) m where getBounds (UMA l u _) = return (l,u) instance (STorIO m s, Unboxed e) => MArray (UnboxedMutableArray s) e m where newArray_ (l,u) = do arr <- allocUnboxed (rangeSize (l,u)) return (UMA l u arr) unsafeRead (UMA _ _ arr) index = readUnboxed arr index unsafeWrite (UMA _ _ arr) index = writeUnboxed arr index -- | Unboxed mutable arrays in ST monad type STUArray = UnboxedMutableArray -- | Unboxed mutable arrays in IO monad type IOUArray = IOSpecific3 UnboxedMutableArray my library also provides monad-independent references. i.e. you can write monadic code that works with references and this code can be runned without any problems both in IO and ST monads: -- This section demonstrates running of monad-independent algorithm -- `test_Ref` in IO and ST monads test_Ref 3 >>= print print $ runST (test_Ref 4) i can also add support for monad-independent array manipulations so that you can write code that will work in both monads. it was my old idea but i had no clients for it :) MP> Now to my question: I would like to represent a matrix as a wrapper around a MP> block, which in turn is just an unboxed array. Here are the definitions for a MP> matrix in ST and outside of a monad, respectively: MP> type MBlock s = STUArray s Int Double MP> data MMatrix s = MMatrix Int Int (MBlock s) MP> type Block = UArray Int Double MP> data Matrix = Matrix Int Int Block seems that you don't know that Haskell's indexes can be a tuples :) type Matrix a = UArray (Int,Int) a type Matrix3d a = UArray (Int,Int,Int) a MP> Now, I have started by providing some fairly low-level routines in the ST MP> monad. My problem is depressingly simple: I would like to retrieve a matrix MP> from the ST monad for use outside of it (in particular, to pretty-print it). MP> Now, this is easy for STUArrays (just use runSTUArray), but I'm not sure of how MP> to do it for a type that *encloses* an STUArray. may be you need to look inside runSTUArray's definition? ;) that is from one doc i wrote: "Operations that creates/updates immutable arrays just creates them as mutable arrays in ST monad, make all required updates on this array and then use unsafeFreeze before returning array from runST." so, you can return anything else together with array returned by unsafeFreeze. of course, because you should use tuples for indexing, it has only theoretical interest. but i strongly recommend you to read entire data.array.* sources to know about all intrinsics of arrays library and in particular freeze/thaw tale MP> runSTMatrix :: ST s (Matrix s) -> Matrix runSTMatrix :: ST s (MMatrix s) -> Matrix runSTMatrix a = runST ( do (MMatrix i j mblock) <- a block <- unsafeFreeze mblock return (Matrix i j block) ) -- Best regards, Bulat mailto:Bulat.Ziganshin@gmail.com