On Sun, Mar 12, 2006 at 10:37:45PM +0300, Bulat Ziganshin wrote:
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 :)
What a pity :-( Probably the best thing to do is to extend Alberto's library to do in-monad versions of GSL's updating algorithms. They would have to be in the IO monad, as I'm calling foreign code, but I could use the unsafeIOToST trick you mention.
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? :)
No, because I did know about this. But precisely because IO monad is ST specialized to RealWorld, the "correct" (i.e. elegant) way of writing the library would be as I suggest; from ST to IO, not the other way around.
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.
Now this *is* funny! ;-)
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
It's true that your version is slightly cleaner.
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.
I *do* use tuples for indexing (even though I don't use them [yet] for *specification* of the dimension of the matrix). Indeed: infixr 2 !@ infixr 3 <-- (?@) :: MMatrix s -> (Int, Int) -> ST s Double m ?@ b = readArray (getMBlock m) ((mrows m)*((fst b)-1) + snd b) (!@) :: MMatrix s -> ((Int, Int), Double) -> ST s () m !@ ((i,j), m_ij) = let n = mrows m in writeArray (getMBlock m) (n*(i-1) + j) m_ij (<--) :: (Int, Int) -> Double -> ((Int, Int), Double) i <-- x = (i, x) which allows code like: runMatrix = do _A <- newListMatrix [[2, 0, 0], [0, 2, 0], [0, 0, 2]] _B <- newListMatrix [[1, 2, 1], [0, 1, 1], [0, 0, 1]] b_12 <- _B ?@ (1,2) _B !@ (1,1) <-- 2*b_12 _C <- matMul _A _B return _C which I think looks quite readable -- and without needing to "cheat" by using outside parsers ;-)
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) )
I tried this implementation, but I still get an error message, which looks quite similar to my previous implementations' errors: matrix.hs:138:27: Couldn't match the rigid variable `s' against the rigid variable `s1' `s' is bound by the polymorphic type `forall s. ST s a' at matrix.hs:(138,16)-(141,22) `s1' is bound by the type signature for `runSTMatrix' Expected type: ST s Inferred type: ST s1 In a 'do' expression: (MMatrix i j mblock) <- a In the first argument of `runST', namely `(do (MMatrix i j mblock) <- a block <- unsafeFreeze mblock return (Matrix i j block))'
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
I will do, and thanks for your excellent advice! Martin