[Haskell-cafe] More STUArray questions

[Martin Percossi]

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