[Haskell-cafe] More STUArray questions

[Chris Kuklewicz]

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.html#v%3AunsafeFreeze

"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