[Haskell-cafe] More STUArray questions

[Bulat Ziganshin]

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 at gmail.com