3 files changed, 171 insertions, 9 deletions

diff --git a/examples/inplace.hs b/examples/inplace.hs
new file mode 100644
index 0000000..d32933a
--- /dev/null
+++ b/examples/inplace.hs
@@ -0,0 +1,118 @@
+-- some tests of the interface for pure
+-- computations with inplace updates
+
+import Numeric.LinearAlgebra
+import Data.Packed.ST
+import Data.Packed.Convert
+
+import Data.Array.Unboxed
+import Data.Array.ST
+import Control.Monad.ST
+import Control.Monad
+
+main = sequence_[
+    print test1,
+    print test2,
+    print test3,
+    print test4,
+    test5,
+    test6,
+    print test7,
+    test0]
+
+-- helper functions
+vector l = fromList l :: Vector Double
+norm v = pnorm PNorm2 v
+zeros n m = reshape m (constant 0 (n*m))
+newMatrix n m = thawMatrix (zeros n m)
+
+
+-- hmatrix vector and matrix
+v = vector [1..10]
+m = (5><10) [1..50::Double]
+
+----------------------------------------------------------------------
+
+-- vector creation by in-place updates on a copy of the argument
+test1 = fun v
+
+fun :: Element t => Vector t -> Vector t
+fun x = runSTVector $ do
+    a <- thawVector x
+    mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]
+    return a
+
+-- another example: creation of an antidiagonal matrix from a list
+test2 = antiDiag 5 8 [1..] :: Matrix Double
+
+antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
+antiDiag r c l = runSTMatrix $ do
+    m <- newMatrix r c
+    let d = min r c - 1
+    sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l
+    return m
+
+-- using vector or matrix functions on mutable objects requires freezing:
+test3 = g1 v
+
+g1 x = runST $ do
+    a <- thawVector x
+    writeVector a (dim x -1) 0
+    b <- freezeVector a
+    return (norm b)
+
+--  another possibility:
+test4 = g2 v
+
+g2 x = runST $ do
+    a <- thawVector x
+    writeVector a (dim x -1) 0
+    t <- liftSTVector norm a
+    return t
+
+--------------------------------------------------------------
+
+-- haskell arrays
+hv = listArray (0,9) [1..10::Double]
+hm = listArray ((0,0),(4,9)) [1..50::Double]
+
+
+
+-- conversion from standard Haskell arrays
+test5 = do
+    print $ norm (vectorFromArray hv)
+    print $ norm v
+    print $ rcond (matrixFromArray hm)
+    print $ rcond m
+
+
+-- conversion to mutable ST arrays
+test6 = do
+    let y = clearColumn m 1
+    print y
+    print (matrixFromArray y)
+
+clearColumn x c = runSTUArray $ do
+    a <- mArrayFromMatrix x
+    forM_ [0..rows x-1] $ \i->
+        writeArray a (i,c) (0::Double)
+    return a
+
+-- hmatrix functions applied to mutable ST arrays
+test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)
+
+unitary v = runSTUArray $ do
+    a <- thaw v
+    n <- norm `fmap` vectorFromMArray a
+    b <- mapArray (/n) a
+    return b
+
+-------------------------------------------------
+
+-- (just to check that they are not affected)
+test0 = do
+    print v
+    print m
+    --print hv
+    --print hm
+
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 1dcc21b..5efe7a6 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -86,6 +86,7 @@ library
                         Graphics.Plot,
                         Numeric.LinearAlgebra.Tests
                         Data.Packed.Convert
+                        Data.Packed.ST
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Vector,
diff --git a/lib/Data/Packed/Convert.hs b/lib/Data/Packed/Convert.hs
index 55fdf8f..206fa06 100644
--- a/lib/Data/Packed/Convert.hs
+++ b/lib/Data/Packed/Convert.hs
@@ -1,3 +1,5 @@
+{-# OPTIONS -XTypeOperators -XRank2Types  -XFlexibleContexts #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Data.Packed.Convert
@@ -9,18 +11,26 @@
 -- Portability :  portable
 --
 -- Conversion of Vectors and Matrices to and from the standard Haskell arrays.
+-- (provisional)
 --
 -----------------------------------------------------------------------------
 
 module Data.Packed.Convert (
-    vectorToStorableArray,
-    storableArrayToVector,
---    unsafeUpdateVector,
+    arrayFromVector, vectorFromArray,
+    mArrayFromVector, vectorFromMArray,
+    vectorToStorableArray, storableArrayToVector,
+    arrayFromMatrix, matrixFromArray,
+    mArrayFromMatrix, matrixFromMArray,
+--    matrixToStorableArray, storableArrayToMatrix
 ) where
 
 import Data.Packed.Internal
 import Data.Array.Storable
 import Foreign
+import Control.Monad.ST
+import Data.Array.ST
+import Data.Array.IArray
+import Data.Array.Unboxed
 
 -- | Creates a StorableArray indexed from 0 to dim -1.
 --   (Memory is efficiently copied, so you can then freely modify the obtained array)
@@ -49,11 +59,44 @@ unsafeVectorToStorableArray v = unsafeForeignPtrToStorableArray (fptr v) (0,dim
 --unsafeStorableArrayToVector s = undefined -- the foreign ptr of Storable Array is not available?
 
 -----------------------------------------------------------------
+-- provisional, we need Unboxed arrays for Complex Double
+
+
+unsafeFreeze' :: (MArray a e m, Ix i) => a i e -> m (Array i e)
+unsafeFreeze' = unsafeFreeze
+
+-- | creates an immutable Array from an hmatrix Vector (to do: unboxed)
+arrayFromVector :: (Storable t) => Vector t -> Array Int t
+arrayFromVector x = runSTArray (mArrayFromVector x)
+
+-- | creates a mutable array from an hmatrix Vector (to do: unboxed)
+mArrayFromVector :: (MArray b t (ST s), Storable t) => Vector t -> ST s (b Int t)
+mArrayFromVector v = unsafeThaw =<< unsafeIOToST ( unsafeFreeze' =<< (vectorToStorableArray $ v))
+
+
+-- (creates an hmatrix Vector from an immutable array (to do: unboxed))
+vectorFromArray :: (Storable t) => Array Int t -> Vector t
+vectorFromArray a = unsafePerformIO $ storableArrayToVector =<< unsafeThaw a
+
+-- | creates a mutable Array from an hmatrix Vector for manipulation with runSTUArray (to do: unboxed)
+vectorFromMArray :: (Storable t, MArray a t (ST s)) => a Int t -> ST s (Vector t)
+vectorFromMArray x = fmap vectorFromArray (unsafeFreeze' x)
+
+--------------------------------------------------------------------
+-- provisional
+
+matrixFromArray :: UArray (Int, Int) Double -> Matrix Double
+matrixFromArray m = reshape c . fromList . elems $ m
+    where ((r1,c1),(r2,c2)) = bounds m
+          r = r2-r1+1
+          c = c2-c1+1
+
+arrayFromMatrix :: Matrix Double -> UArray (Int, Int) Double
+arrayFromMatrix m = listArray ((0,0),(rows m -1, cols m -1)) (toList $ flatten m)
+
 
+mArrayFromMatrix :: (MArray b Double m) => Matrix Double -> m (b (Int, Int) Double)
+mArrayFromMatrix = unsafeThaw . arrayFromMatrix
 
-unsafeUpdateVector :: Storable t => Int -- ^ position
-                           -> (t->t)    -- ^ modification function
-                           -> Vector t  -- ^ source
-                           -> IO ()     -- ^ result
-unsafeUpdateVector k h v = do
-    withForeignPtr (fptr v) $ \s -> pokeElemOff s k (h (v`at'`k))
+matrixFromMArray :: (MArray a Double (ST s)) => a (Int,Int) Double -> ST s (Matrix Double)
+matrixFromMArray x = fmap matrixFromArray (unsafeFreeze x)