9 files changed, 280 insertions, 92 deletions

diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index aaf6fbb..1235da3 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -29,7 +29,9 @@ import Internal.Conversion
 import Internal.LAPACK as LAPACK
 import Internal.Numeric
 import Data.List(foldl1')
-import Data.Array
+import qualified Data.Array as A
+import Internal.ST
+import Internal.Vectorized(range)
 
 {- | Generic linear algebra functions for double precision real and complex matrices.
 
@@ -578,11 +580,6 @@ eps =  2.22044604925031e-16
 peps :: RealFloat x => x
 peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)
 
-
--- | The imaginary unit: @i = 0.0 :+ 1.0@
-i :: Complex Double
-i = 0:+1
-
 -----------------------------------------------------------------------
 
 -- | The nullspace of a matrix from its precomputed SVD decomposition.
@@ -796,13 +793,23 @@ signlp r vals = foldl f 1 (zip [0..r-1] vals)
     where f s (a,b) | a /= b    = -s
                     | otherwise =  s
 
-swap (arr,s) (a,b) | a /= b    = (arr // [(a, arr!b),(b,arr!a)],-s)
-                   | otherwise = (arr,s)
-
-fixPerm r vals = (fromColumns $ elems res, sign)
-    where v = [0..r-1]
-          s = toColumns (ident r)
-          (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
+fixPerm r vals = (fromColumns $ A.elems res, sign)
+  where
+    v = [0..r-1]
+    t = toColumns (ident r)
+    (res,sign) = foldl swap (A.listArray (0,r-1) t, 1) (zip v vals)
+    swap (arr,s) (a,b)
+      | a /= b    = (arr A.// [(a, arr A.! b),(b,arr A.! a)],-s)
+      | otherwise = (arr,s)
+
+fixPerm' :: [Int] -> Vector I
+fixPerm' s = res $ mutable f s0
+  where
+    s0 = reshape 1 (range (length s))
+    res = flatten . fst
+    swap m i j = rowOper (SWAP i j AllCols) m
+    f :: (Num t, Element t) => (Int, Int) -> STMatrix s t -> ST s () -- needed because of TypeFamilies
+    f _ p = sequence_ $ zipWith (swap p) [0..] s
 
 triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
     where el p q = if q-p>=h then v else 1 - v
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 2843ab5..4d48594 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1398,6 +1398,65 @@ ROWOP(int64_t)
 ROWOP_MOD(int32_t,mod)
 ROWOP_MOD(int64_t,mod_l)
 
+/////////////////////////////// inplace GEMM ////////////////////////////////
+
+#define GEMM(T) int gemm_##T(VECG(T,c),VECG(int,p),MATG(T,a),MATG(T,b),MATG(T,r)) {  \
+    T a = cp[0], b = cp[1];                                    \
+    int r1a = pp[0],              c1a = pp[2],  c2a = pp[3] ;  \
+    int r1b = pp[4],              c1b = pp[6]               ;  \
+    int r1r = pp[8], r2r = pp[9], c1r = pp[10], c2r = pp[11];  \
+    int dra = r1a - r1r;                                       \
+    int dcb = c1b-c1r;                                         \
+    int nk  = c2a-c1a+1;                                       \
+    int i,j,k;                                                 \
+    T t;                                                       \
+    for (i=r1r; i<=r2r; i++) {                                 \
+        for (j=c1r; j<=c2r; j++) {                             \
+            t = 0;                                             \
+            for(k=0; k<nk; k++) {                              \
+                t += AT(a,i+dra,k+c1a) * AT(b,k+r1b,j+dcb);    \
+            }                                                  \
+            AT(r,i,j) = b*AT(r,i,j) + a*t;                     \
+        }                                                      \
+    }                                                          \
+    OK                                                         \
+}
+
+GEMM(double)
+GEMM(float)
+GEMM(TCD)
+GEMM(TCF)
+GEMM(int32_t)
+GEMM(int64_t)
+
+#define GEMM_MOD(T,M) int gemm_mod_##T(T m, VECG(T,c),VECG(int,p),MATG(T,a),MATG(T,b),MATG(T,r)) {  \
+    T a = cp[0], b = cp[1];                                    \
+    int r1a = pp[0],              c1a = pp[2],  c2a = pp[3] ;  \
+    int r1b = pp[4],              c1b = pp[6]               ;  \
+    int r1r = pp[8], r2r = pp[9], c1r = pp[10], c2r = pp[11];  \
+    int dra = r1a - r1r;                                       \
+    int dcb = c1b-c1r;                                         \
+    int nk  = c2a-c1a+1;                                       \
+    int i,j,k;                                                 \
+    T t;                                                       \
+    for (i=r1r; i<=r2r; i++) {                                 \
+        for (j=c1r; j<=c2r; j++) {                             \
+            t = 0;                                             \
+            for(k=0; k<nk; k++) {                                     \
+                t = M(t+M(AT(a,i+dra,k+c1a) * AT(b,k+r1b,j+dcb)));    \
+            }                                                         \
+            AT(r,i,j) = M(M(b*AT(r,i,j)) + M(a*t));                   \
+        }                                                      \
+    }                                                          \
+    OK                                                         \
+}
+
+#define MOD32(X) mod(X,m)
+#define MOD64(X) mod_l(X,m)
+
+GEMM_MOD(int32_t,MOD32)
+GEMM_MOD(int64_t,MOD64)
+
 ////////////////// sparse matrix-product ///////////////////////////////////////
 
 
diff --git a/packages/base/src/Internal/C/lapack-aux.h b/packages/base/src/Internal/C/lapack-aux.h
index e4d95bc..bf8c5e9 100644
--- a/packages/base/src/Internal/C/lapack-aux.h
+++ b/packages/base/src/Internal/C/lapack-aux.h
@@ -59,6 +59,7 @@ typedef short ftnlen;
 #define OQMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, complex* A##p
 #define OCMAT(A) int A##r, int A##c, int A##Xr, int A##Xc, doublecomplex* A##p
 
+#define VECG(T,A) int A##n, T* A##p
 #define MATG(T,A) int A##r, int A##c, int A##Xr, int A##Xc, T* A##p
 
 #define KIVEC(A) int A##n, const int*A##p
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index 75e92a5..8f8c219 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -274,6 +274,7 @@ class (Storable a) => Element a where
     selectV  :: Vector CInt -> Vector a -> Vector a -> Vector a -> Vector a
     remapM   :: Matrix CInt -> Matrix CInt -> Matrix a -> Matrix a
     rowOp    :: Int -> a -> Int -> Int -> Int -> Int -> Matrix a -> IO ()
+    gemm     :: Vector a -> Vector I -> Matrix a -> Matrix a -> Matrix a -> IO ()
 
 
 instance Element Float where
@@ -287,6 +288,7 @@ instance Element Float where
     selectV    = selectF
     remapM     = remapF
     rowOp      = rowOpAux c_rowOpF
+    gemm       = gemmg c_gemmF
 
 instance Element Double where
     transdata  = transdataAux ctransR
@@ -299,7 +301,7 @@ instance Element Double where
     selectV    = selectD
     remapM     = remapD
     rowOp      = rowOpAux c_rowOpD
-
+    gemm       = gemmg c_gemmD
 
 instance Element (Complex Float) where
     transdata  = transdataAux ctransQ
@@ -312,7 +314,7 @@ instance Element (Complex Float) where
     selectV    = selectQ
     remapM     = remapQ
     rowOp      = rowOpAux c_rowOpQ
-
+    gemm       = gemmg c_gemmQ
 
 instance Element (Complex Double) where
     transdata  = transdataAux ctransC
@@ -325,6 +327,7 @@ instance Element (Complex Double) where
     selectV    = selectC
     remapM     = remapC
     rowOp      = rowOpAux c_rowOpC
+    gemm       = gemmg c_gemmC
 
 instance Element (CInt) where
     transdata  = transdataAux ctransI
@@ -337,6 +340,7 @@ instance Element (CInt) where
     selectV    = selectI
     remapM     = remapI
     rowOp      = rowOpAux c_rowOpI
+    gemm       = gemmg c_gemmI
 
 instance Element Z where
     transdata  = transdataAux ctransL
@@ -349,6 +353,7 @@ instance Element Z where
     selectV    = selectL
     remapM     = remapL
     rowOp      = rowOpAux c_rowOpL
+    gemm       = gemmg c_gemmL
 
 -------------------------------------------------------------------
 
@@ -575,6 +580,21 @@ foreign import ccall unsafe "rowop_mod_int64_t" c_rowOpML :: Z -> RowOp Z
 
 --------------------------------------------------------------------------------
 
+gemmg f u v m1 m2 m3 = app5 f vec u vec v omat m1 omat m2 omat m3 "gemmg"
+
+type Tgemm x = x :> I :> x ::> x ::> x ::> Ok
+
+foreign import ccall unsafe "gemm_double"  c_gemmD :: Tgemm R
+foreign import ccall unsafe "gemm_float"   c_gemmF :: Tgemm Float
+foreign import ccall unsafe "gemm_TCD"     c_gemmC :: Tgemm C
+foreign import ccall unsafe "gemm_TCF"     c_gemmQ :: Tgemm (Complex Float)
+foreign import ccall unsafe "gemm_int32_t" c_gemmI :: Tgemm I
+foreign import ccall unsafe "gemm_int64_t" c_gemmL :: Tgemm Z
+foreign import ccall unsafe "gemm_mod_int32_t" c_gemmMI :: I -> Tgemm I
+foreign import ccall unsafe "gemm_mod_int64_t" c_gemmML :: Z -> Tgemm Z
+
+--------------------------------------------------------------------------------
+
 foreign import ccall unsafe "saveMatrix" c_saveMatrix
     :: CString -> CString -> Double ..> Ok
 
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 6c6d5c5..1cae1ac 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -34,7 +34,9 @@ import Internal.Container
 import Internal.Vectorized (prodI,sumI,prodL,sumL)
 import Internal.LAPACK (multiplyI, multiplyL)
 import Internal.Algorithms(luFact)
-import Internal.Util(Normed(..),Indexable(..),gaussElim, gaussElim_1, gaussElim_2,luST, magnit)
+import Internal.Util(Normed(..),Indexable(..),
+                     gaussElim, gaussElim_1, gaussElim_2,
+                     luST, luSolve', luPacked', magnit)
 import Internal.ST(mutable)
 import GHC.TypeLits
 import Data.Proxy(Proxy)
@@ -131,6 +133,9 @@ instance KnownNat m => Element (Mod m I)
     rowOp c a i1 i2 j1 j2 x = rowOpAux (c_rowOpMI m') c (unMod a) i1 i2 j1 j2 (f2iM x)
       where
         m' = fromIntegral . natVal $ (undefined :: Proxy m)
+    gemm u p a b c = gemmg (c_gemmMI m') (f2i u) p (f2iM a) (f2iM b) (f2iM c)
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m)
 
 instance KnownNat m => Element (Mod m Z)
   where
@@ -146,6 +151,9 @@ instance KnownNat m => Element (Mod m Z)
     rowOp c a i1 i2 j1 j2 x = rowOpAux (c_rowOpML m') c (unMod a) i1 i2 j1 j2 (f2iM x)
       where
         m' = fromIntegral . natVal $ (undefined :: Proxy m)
+    gemm u p a b c = gemmg (c_gemmML m') (f2i u) p (f2iM a) (f2iM b) (f2iM c)
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m)
 
 
 instance forall m . KnownNat m => Container Vector (Mod m I)
@@ -344,7 +352,11 @@ test = (ok, info)
     lg = (3><3) (repeat (3*10^(9::Int))) :: Matrix Z
     lgm = fromZ lg :: Matrix (Mod 10000000000 Z)
 
-    gen n = diagRect 1 (konst 5 n) n n :: Numeric t => Matrix t
+    gen  n = diagRect 1 (konst 5 n) n n :: Numeric t => Matrix t
+    
+    rgen n = gen n :: Matrix R
+    cgen n = complex (rgen n) + fliprl (complex (rgen n)) * scalar (0:+1) :: Matrix C
+    sgen n = single (cgen n)
     
     checkGen x = norm_Inf $ flatten $ invg x <> x - ident (rows x)
     
@@ -354,6 +366,11 @@ test = (ok, info)
       where
         (l,u,p,_ :: Int) = luFact $ mutable (luST okf) t
 
+    checkSolve aa = norm_Inf $ flatten (aa <> x - bb)
+       where
+         bb = flipud aa
+         x = luSolve' (luPacked' aa) bb
+
     info = do
         print v
         print m
@@ -377,9 +394,9 @@ test = (ok, info)
         print $ lgm <> lgm
         
         print (checkGen (gen 5 :: Matrix R))
-        print (checkGen (gen 5 :: Matrix C))
         print (checkGen (gen 5 :: Matrix Float))
-        print (checkGen (gen 5 :: Matrix (Complex Float)))
+        print (checkGen (cgen 5 :: Matrix C))
+        print (checkGen (sgen 5 :: Matrix (Complex Float)))
         print (invg (gen 5) :: Matrix (Mod 7 I))
         print (invg (gen 5) :: Matrix (Mod 7 Z))
         
@@ -388,11 +405,18 @@ test = (ok, info)
 
         print $ checkLU (magnit 0) (gen 5 :: Matrix R)
         print $ checkLU (magnit 0) (gen 5 :: Matrix Float)
-        print $ checkLU (magnit 0) (gen 5 :: Matrix C)
-        print $ checkLU (magnit 0) (gen 5 :: Matrix (Complex Float))
+        print $ checkLU (magnit 0) (cgen 5 :: Matrix C)
+        print $ checkLU (magnit 0) (sgen 5 :: Matrix (Complex Float))
         print $ checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 I))
         print $ checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 Z))
 
+        print $ checkSolve (gen 5 :: Matrix R)
+        print $ checkSolve (gen 5 :: Matrix Float)
+        print $ checkSolve (cgen 5 :: Matrix C)
+        print $ checkSolve (sgen 5 :: Matrix (Complex Float))
+        print $ checkSolve (gen 5 :: Matrix (Mod 7 I))
+        print $ checkSolve (gen 5 :: Matrix (Mod 7 Z))
+
 
     ok = and
       [ toInt (m #> v) == cmod 11 (toInt m #> toInt v )
@@ -401,16 +425,22 @@ test = (ok, info)
       , am <> gaussElim   am bm == bm
       , (checkGen (gen 5 :: Matrix R)) < 1E-15
       , (checkGen (gen 5 :: Matrix Float)) < 2E-7
-      , (checkGen (gen 5 :: Matrix C)) < 1E-15
-      , (checkGen (gen 5 :: Matrix (Complex Float))) < 2E-7
+      , (checkGen (cgen 5 :: Matrix C)) < 1E-15
+      , (checkGen (sgen 5 :: Matrix (Complex Float))) < 2E-7
       , (checkGen (gen 5 :: Matrix (Mod 7 I))) == 0
       , (checkGen (gen 5 :: Matrix (Mod 7 Z))) == 0
-      , (checkLU (magnit 1E-10) (gen 5 :: Matrix R)) < 1E-15
+      , (checkLU (magnit 1E-10) (gen 5 :: Matrix R)) < 2E-15
       , (checkLU (magnit 1E-5) (gen 5 :: Matrix Float)) < 1E-6
-      , (checkLU (magnit 1E-10) (gen 5 :: Matrix C)) < 1E-15
-      , (checkLU (magnit 1E-5) (gen 5 :: Matrix (Complex Float))) < 1E-6
+      , (checkLU (magnit 1E-10) (cgen 5 :: Matrix C)) < 5E-15
+      , (checkLU (magnit 1E-5) (sgen 5 :: Matrix (Complex Float))) < 1E-6
       , (checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 I))) == 0
       , (checkLU (magnit 0) (gen 5 :: Matrix (Mod 7 Z))) == 0
+      , checkSolve (gen 5 :: Matrix R) < 2E-15
+      , checkSolve (gen 5 :: Matrix Float) < 1E-6
+      , checkSolve (cgen 5 :: Matrix C) < 4E-15
+      , checkSolve (sgen 5 :: Matrix (Complex Float)) < 1E-6
+      , checkSolve (gen 5 :: Matrix (Mod 7 I)) == 0
+      , checkSolve (gen 5 :: Matrix (Mod 7 Z)) == 0
       , prodElements (konst (9:: Mod 10 I) (12::Int)) == product (replicate 12 (9:: Mod 10 I))
       , gm <> gm == konst 0 (3,3)
       , lgm <> lgm == konst 0 (3,3)
diff --git a/packages/base/src/Internal/ST.hs b/packages/base/src/Internal/ST.hs
index 434fe63..d1defda 100644
--- a/packages/base/src/Internal/ST.hs
+++ b/packages/base/src/Internal/ST.hs
@@ -1,5 +1,6 @@
 {-# LANGUAGE Rank2Types    #-}
 {-# LANGUAGE BangPatterns  #-}
+{-# LANGUAGE ViewPatterns #-}
 
 -----------------------------------------------------------------------------
 -- |
@@ -15,14 +16,14 @@
 -----------------------------------------------------------------------------
 
 module Internal.ST (
+    ST, runST,
     -- * Mutable Vectors
     STVector, newVector, thawVector, freezeVector, runSTVector,
     readVector, writeVector, modifyVector, liftSTVector,
     -- * Mutable Matrices
     STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
     readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
---  axpy, scal, swap, rowOp,
-    mutable, extractMatrix, setMatrix, rowOper, RowOper(..), RowRange(..), ColRange(..),
+    mutable, extractMatrix, setMatrix, rowOper, RowOper(..), RowRange(..), ColRange(..), gemmm, Slice(..),
     -- * Unsafe functions
     newUndefinedVector,
     unsafeReadVector, unsafeWriteVector,
@@ -70,13 +71,13 @@ unsafeWriteVector  (STVector x) k = unsafeIOToST . ioWriteV x k
 modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s ()
 modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k
 
-liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a
+liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s t -> ST s a
 liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x
 
-freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
+freezeVector :: (Storable t) => STVector s t -> ST s (Vector t)
 freezeVector v = liftSTVector id v
 
-unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
+unsafeFreezeVector :: (Storable t) => STVector s t -> ST s (Vector t)
 unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
 
 {-# INLINE safeIndexV #-}
@@ -139,14 +140,14 @@ unsafeWriteMatrix  (STMatrix x) r c = unsafeIOToST . ioWriteM x r c
 modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
 modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c
 
-liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a
+liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s t -> ST s a
 liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x
 
-unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
+unsafeFreezeMatrix :: (Storable t) => STMatrix s t -> ST s (Matrix t)
 unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
 
 
-freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
+freezeMatrix :: (Storable t) => STMatrix s t -> ST s (Matrix t)
 freezeMatrix m = liftSTMatrix id m
 
 cloneMatrix (Matrix r c d o) = cloneVector d >>= return . (\d' -> Matrix r c d' o)
@@ -227,6 +228,18 @@ extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR m 0 (idxs[i1,i2]) 0 (i
     (i1,i2) = getRowRange (rows m) rr
     (j1,j2) = getColRange (cols m) rc
 
+-- | r0 c0 height width
+data Slice s t = Slice (STMatrix s t) Int Int Int Int
+
+slice (Slice (STMatrix m) r0 c0 nr nc) = (m, idxs[r0,r0+nr-1,c0,c0+nc-1])
+
+gemmm :: Element t => t -> Slice s t -> t -> Slice s t -> Slice s t -> ST s ()
+gemmm beta (slice->(r,pr)) alpha (slice->(a,pa)) (slice->(b,pb)) = res
+  where
+    res = unsafeIOToST (gemm u v a b r)
+    u = fromList [alpha,beta]
+    v = vjoin[pa,pb,pr]
+    
 
 mutable :: Storable t => (forall s . (Int, Int) -> STMatrix s t -> ST s u) -> Matrix t -> (Matrix t,u)
 mutable f a = runST $ do
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index f08f710..079663d 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -1,6 +1,5 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FunctionalDependencies #-}
 {-# LANGUAGE ViewPatterns #-}
@@ -55,7 +54,7 @@ module Internal.Util(
     -- ** 2D
     corr2, conv2, separable,
     block2x2,block3x3,view1,unView1,foldMatrix,
-    gaussElim_1, gaussElim_2, gaussElim, luST
+    gaussElim_1, gaussElim_2, gaussElim, luST, luSolve', luSolve'', luPacked'
 ) where
 
 import Internal.Vector
@@ -65,7 +64,7 @@ import Internal.Element
 import Internal.Container
 import Internal.Vectorized
 import Internal.IO
-import Internal.Algorithms hiding (i,Normed,swap,linearSolve')
+import Internal.Algorithms hiding (Normed,linearSolve',luSolve', luPacked')
 import Numeric.Matrix()
 import Numeric.Vector()
 import Internal.Random
@@ -73,7 +72,7 @@ import Internal.Convolution
 import Control.Monad(when,forM_)
 import Text.Printf
 import Data.List.Split(splitOn)
-import Data.List(intercalate,sortBy)
+import Data.List(intercalate,sortBy,foldl')
 import Control.Arrow((&&&))
 import Data.Complex
 import Data.Function(on)
@@ -687,6 +686,115 @@ luST ok (r,_) x = do
     v <- unsafeFreezeVector p
     return (toList v)
 
+{- | Experimental implementation of 'luPacked'
+     for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
+
+>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
+(5><5)
+ [  1,  1,  2,  3,  4
+ ,  5,  7,  7,  8,  9
+ , 10, 11, 13, 13, 14
+ , 15, 16,  0,  2,  2
+ ,  3,  4,  5,  6,  8 ]
+
+>>> let (l,u,p,s) = luFact $ luPacked' m
+>>> l
+(5><5)
+ [  1,  0, 0,  0, 0
+ ,  6,  1, 0,  0, 0
+ , 12,  7, 1,  0, 0
+ ,  7, 10, 7,  1, 0
+ ,  8,  2, 6, 11, 1 ]
+>>> u
+(5><5)
+ [ 15, 16,  0,  2,  2
+ ,  0, 13,  7, 13, 14
+ ,  0,  0, 15,  0, 11
+ ,  0,  0,  0, 15, 15
+ ,  0,  0,  0,  0,  1 ]
+
+-}
+luPacked' x = mutable (luST (magnit 0)) x
+
+--------------------------------------------------------------------------------
+
+rowRange m = [0..rows m -1]
+
+at k = Pos (idxs[k])
+
+backSust' lup rhs = foldl' f (rhs?[]) (reverse ls)
+  where
+    ls  = [ (d k , u k , b k) | k <- rowRange lup ]
+      where
+        d k = lup ?? (at k, at k)
+        u k = lup ?? (at k, Drop (k+1))
+        b k = rhs ?? (at k, All)
+
+    f x (d,u,b) = (b - u<>x) / d
+                       ===
+                        x
+
+
+forwSust' lup rhs = foldl' f (rhs?[]) ls
+  where
+    ls  = [ (l k , b k) | k <- rowRange lup ]
+      where
+        l k = lup ?? (at k, Take k)
+        b k = rhs ?? (at k, All)
+
+    f x (l,b) =     x
+                   ===
+                (b - l<>x)
+
+
+luSolve'' (lup,p) b = backSust' lup (forwSust' lup pb)
+  where
+    pb = b ?? (Pos (fixPerm' p), All)
+
+--------------------------------------------------------------------------------
+
+forwSust lup rhs = fst $ mutable f rhs
+  where
+    f (r,c) x = do
+        l <- unsafeThawMatrix lup
+        let go k = gemmm 1 (Slice x k 0 1 c) (-1) (Slice l k 0 1 k) (Slice x 0 0 k c)
+        mapM_ go [0..r-1]
+
+
+backSust lup rhs = fst $ mutable f rhs
+  where
+    f (r,c) m = do
+        l <- unsafeThawMatrix lup
+        let d k = recip (lup `atIndex` (k,k))
+            u k = Slice l k (k+1) 1 (r-1-k)
+            b k = Slice m k 0 1 c
+            x k = Slice m (k+1) 0 (r-1-k) c
+            scal k = rowOper (SCAL (d k) (Row k) AllCols) m
+
+            go k = gemmm 1 (b k) (-1) (u k) (x k) >> scal k
+        mapM_ go [r-1,r-2..0]
+
+
+{- | Experimental implementation of 'luSolve' for any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
+
+>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
+(2><2)
+ [ 1, 2
+ , 3, 5 ]
+>>> b
+(2><3)
+ [ 5, 1, 3
+ , 8, 6, 3 ]
+
+>>> luSolve' (luPacked' a) b
+(2><3)
+ [ 4,  7, 4
+ , 7, 10, 6 ]
+
+-}
+luSolve' (lup,p) b = backSust lup (forwSust lup pb)
+  where
+    pb = b ?? (Pos (fixPerm' p), All)
 
 --------------------------------------------------------------------------------
 
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 2e6b8ca..0b8abbb 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -77,10 +77,10 @@ module Numeric.LinearAlgebra (
     linearSolveLS,
     linearSolveSVD,
     luSolve,
+    luSolve',
     cholSolve,
     cgSolve,
     cgSolve',
-    linearSolve',
 
     -- * Inverse and pseudoinverse
     inv, pinv, pinvTol,
@@ -122,7 +122,7 @@ module Numeric.LinearAlgebra (
     schur,
 
     -- * LU
-    lu, luPacked, luFact, luPacked',
+    lu, luPacked, luPacked', luFact,
 
     -- * Matrix functions
     expm,
@@ -158,14 +158,13 @@ import Numeric.Vector()
 import Internal.Matrix
 import Internal.Container hiding ((<>))
 import Internal.Numeric hiding (mul)
-import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve')
+import Internal.Algorithms hiding (linearSolve,Normed,orth,luPacked',linearSolve',luSolve')
 import qualified Internal.Algorithms as A
 import Internal.Util
 import Internal.Random
 import Internal.Sparse((!#>))
 import Internal.CG
 import Internal.Conversion
-import Internal.ST(mutable)
 
 {- | infix synonym of 'mul'
 
@@ -240,53 +239,4 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
 -- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
 orth m = orthSVD (Left (1*eps)) m (leftSV m)
 
-{- | Experimental implementation of LU factorization
-     working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
-
->>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
-(5><5)
- [  1,  1,  2,  3,  4
- ,  5,  7,  7,  8,  9
- , 10, 11, 13, 13, 14
- , 15, 16,  0,  2,  2
- ,  3,  4,  5,  6,  8 ]
-
->>> let (l,u,p,s) = luFact $ luPacked' m
->>> l
-(5><5)
- [  1,  0, 0,  0, 0
- ,  6,  1, 0,  0, 0
- , 12,  7, 1,  0, 0
- ,  7, 10, 7,  1, 0
- ,  8,  2, 6, 11, 1 ]
->>> u
-(5><5)
- [ 15, 16,  0,  2,  2
- ,  0, 13,  7, 13, 14
- ,  0,  0, 15,  0, 11
- ,  0,  0,  0, 15, 15
- ,  0,  0,  0,  0,  1 ]
-
--}
-luPacked' x = mutable (luST (magnit 0)) x
-
-{- | Experimental implementation of gaussian elimination
-     working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
-
->>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
-(2><2)
- [ 1, 2
- , 3, 5 ]
->>> b
-(2><3)
- [ 5, 1, 3
- , 8, 6, 3 ]
-
->>> let x = linearSolve' a b
-(2><3)
- [ 4,  7, 4
- , 7, 10, 6 ]
-
--}
-linearSolve' x y = gaussElim x y
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
index 36c5f03..db4236b 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
@@ -43,7 +43,7 @@ module Numeric.LinearAlgebra.Devel(
     -- ** Mutable Matrices
     STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
     readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
-    mutable, extractMatrix, setMatrix, rowOper, RowOper(..), RowRange(..), ColRange(..),
+    mutable, extractMatrix, setMatrix, rowOper, RowOper(..), RowRange(..), ColRange(..), gemmm, Slice(..),
     -- ** Unsafe functions
     newUndefinedVector,
     unsafeReadVector, unsafeWriteVector,