7 files changed, 143 insertions, 148 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 4d48594..cdbaab9 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1485,26 +1485,6 @@ int smTXv(KDVEC(vals),KIVEC(cols),KIVEC(rows),KDVEC(x),DVEC(r)) {
 }
 
 
-//////////////////// transpose /////////////////////////
-
-#define TRANS_IMP {                       \
-    REQUIRES(xr==tc && xc==tr,BAD_SIZE);  \
-    DEBUGMSG("trans");                    \
-    int i,j;                              \
-    for (i=0; i<tr; i++) {                \
-        for (j=0; j<tc; j++) {            \
-            tp[i*tc+j] = xp[j*xc+i];      \
-        }                                 \
-    }                                     \
-    OK }
-
-int transF(KFMAT(x),FMAT(t)) TRANS_IMP
-int transR(KDMAT(x),DMAT(t)) TRANS_IMP
-int transQ(KQMAT(x),QMAT(t)) TRANS_IMP
-int transC(KCMAT(x),CMAT(t)) TRANS_IMP
-int transI(KIMAT(x),IMAT(t)) TRANS_IMP
-int transL(KLMAT(x),LMAT(t)) TRANS_IMP
-
 //////////////////////// extract /////////////////////////////////
 
 #define EXTRACT_IMP {                      \
diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs
index 51d5686..6d86f3d 100644
--- a/packages/base/src/Internal/Element.hs
+++ b/packages/base/src/Internal/Element.hs
@@ -173,7 +173,7 @@ m ?? (e, TakeLast n) = m ?? (e, Drop (cols m - n))
 m ?? (DropLast n, e) = m ?? (Take (rows m - n), e)
 m ?? (e, DropLast n) = m ?? (e, Take (cols m - n))
 
-m ?? (er,ec) = unsafePerformIO $ extractR m moder rs modec cs
+m ?? (er,ec) = unsafePerformIO $ extractR (orderOf m) m moder rs modec cs
   where
     (moder,rs) = mkExt (rows m) er
     (modec,cs) = mkExt (cols m) ec
@@ -491,9 +491,13 @@ liftMatrix2Auto f m1 m2
 -- FIXME do not flatten if equal order
 lM f m1 m2 = matrixFromVector
                 RowMajor
-                (max (rows m1) (rows m2))
-                (max (cols m1) (cols m2))
+                (max' (rows m1) (rows m2))
+                (max' (cols m1) (cols m2))
                 (f (flatten m1) (flatten m2))
+  where
+    max' 1 b = b
+    max' a 1 = a
+    max' a b = max a b
 
 compat' :: Matrix a -> Matrix b -> Bool
 compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
@@ -595,6 +599,6 @@ mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
     where
       c = cols m
 
-mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b
+mapMatrix :: (Element a, Element b) => (a -> b) -> Matrix a -> Matrix b
 mapMatrix f = liftMatrix (mapVector f)
 
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index fc9e3ad..5319e95 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -418,6 +418,8 @@ linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
 
 -----------------------------------------------------------------------------------
 
+type TMM t = t ..> t ..> Ok
+
 foreign import ccall unsafe "chol_l_H" zpotrf :: TMM C
 foreign import ccall unsafe "chol_l_S" dpotrf :: TMM R
 
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index f76b9dc..bdf2785 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -32,49 +32,13 @@ import Foreign.C.Types ( CInt(..) )
 import Foreign.C.String ( CString, newCString )
 import System.IO.Unsafe ( unsafePerformIO )
 import Control.DeepSeq ( NFData(..) )
-import Data.List.Split(chunksOf)
+import Text.Printf
 
 -----------------------------------------------------------------
 
-{- Design considerations for the Matrix Type
-   -----------------------------------------
-
-- we must easily handle both row major and column major order,
-  for bindings to LAPACK and GSL/C
-
-- we'd like to simplify redundant matrix transposes:
-   - Some of them arise from the order requirements of some functions
-   - some functions (matrix product) admit transposed arguments
-
-- maybe we don't really need this kind of simplification:
-   - more complex code
-   - some computational overhead
-   - only appreciable gain in code with a lot of redundant transpositions
-     and cheap matrix computations
-
-- we could carry both the matrix and its (lazily computed) transpose.
-  This may save some transpositions, but it is necessary to keep track of the
-  data which is actually computed to be used by functions like the matrix product
-  which admit both orders.
-
-- but if we need the transposed data and it is not in the structure, we must make
-  sure that we touch the same foreignptr that is used in the computation.
-
-- a reasonable solution is using two constructors for a matrix. Transposition just
-  "flips" the constructor. Actual data transposition is not done if followed by a
-  matrix product or another transpose.
-
--}
-
 data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
 
-transOrder RowMajor = ColumnMajor
-transOrder ColumnMajor = RowMajor
-{- | Matrix representation suitable for BLAS\/LAPACK computations.
-
-The elements are stored in a continuous memory array.
-
--}
+-- | Matrix representation suitable for BLAS\/LAPACK computations.
 
 data Matrix t = Matrix
     { irows :: {-# UNPACK #-} !Int
@@ -83,8 +47,6 @@ data Matrix t = Matrix
     , xCol  :: {-# UNPACK #-} !Int
     , xdat  :: {-# UNPACK #-} !(Vector t)
     }
--- RowMajor: preferred by C, fdat may require a transposition
--- ColumnMajor: preferred by LAPACK, cdat may require a transposition
 
 
 rows :: Matrix t -> Int
@@ -95,32 +57,55 @@ cols :: Matrix t -> Int
 cols = icols
 {-# INLINE cols #-}
 
-rowOrder m = xRow m > 1
+size m = (irows m, icols m)
+{-# INLINE size #-}
+
+rowOrder m = xCol m == 1 || cols m == 1
 {-# INLINE rowOrder #-}
 
-isSlice m = cols m < xRow m || rows m < xCol m
+colOrder m = xRow m == 1 || rows m == 1
+{-# INLINE colOrder #-}
+
+is1d (size->(r,c)) = r==1 || c==1
+{-# INLINE is1d #-}
+
+-- data is not contiguous
+isSlice m@(size->(r,c)) = (c < xRow m || r < xCol m) && min r c > 1
 {-# INLINE isSlice #-}
 
 orderOf :: Matrix t -> MatrixOrder
 orderOf m = if rowOrder m then RowMajor else ColumnMajor
 
 
+showInternal :: Storable t => Matrix t -> IO ()
+showInternal m = printf "%dx%d %s %s %d:%d (%d)\n" r c slc ord xr xc dv
+  where
+    r  = rows m
+    c  = cols m
+    xr = xRow m
+    xc = xCol m
+    slc = if isSlice m then "slice" else "full"
+    ord = if is1d m then "1d" else if rowOrder m then "rows" else "cols"
+    dv = dim (xdat m)
+
+--------------------------------------------------------------------------------
+
 -- | Matrix transpose.
 trans :: Matrix t -> Matrix t
-trans m@Matrix { irows = r, icols = c } | rowOrder m =
-             m { irows = c, icols = r, xRow = 1, xCol = c }
-trans m@Matrix { irows = r, icols = c } =
-             m { irows = c, icols = r, xRow = r, xCol = 1 }
+trans m@Matrix { irows = r, icols = c, xRow = xr, xCol = xc } =
+             m { irows = c, icols = r, xRow = xc, xCol = xr }
+
 
 cmat :: (Element t) => Matrix t -> Matrix t
-cmat m | rowOrder m = m
-cmat m@Matrix { irows = r, icols = c, xdat = d } =
-            m { xdat = transdata r d c, xRow = c, xCol = 1 }
+cmat m
+    | rowOrder m = m
+    | otherwise  = extractAll RowMajor m
+
 
 fmat :: (Element t) => Matrix t -> Matrix t
-fmat m | not (rowOrder m) = m
-fmat m@Matrix { irows = r, icols = c, xdat = d} =
-            m { xdat = transdata c d r, xRow = 1, xCol = r }
+fmat m
+    | colOrder m = m
+    | otherwise  = extractAll ColumnMajor m
 
 
 -- C-Haskell matrix adapters
@@ -157,6 +142,11 @@ a # b = apply a b
 
 --------------------------------------------------------------------------------
 
+extractAll ord m = unsafePerformIO $
+    extractR ord m
+             0 (idxs[0,rows m-1])
+             0 (idxs[0,cols m-1])
+
 {- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
 
 >>> flatten (ident 3)
@@ -164,12 +154,14 @@ fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
 
 -}
 flatten :: Element t => Matrix t -> Vector t
-flatten = xdat . cmat
+flatten m
+    | isSlice m || not (rowOrder m) = xdat (extractAll RowMajor m)
+    | otherwise                     = xdat m
 
 
 -- | the inverse of 'Data.Packed.Matrix.fromLists'
 toLists :: (Element t) => Matrix t -> [[t]]
-toLists m = chunksOf (cols m) . toList . flatten $ m
+toLists = map toList . toRows
 
 
 
@@ -205,6 +197,14 @@ fromRows vs = case compatdim (map dim vs) of
 -- | extracts the rows of a matrix as a list of vectors
 toRows :: Element t => Matrix t -> [Vector t]
 toRows m
+    | rowOrder m = map sub rowRange
+    | otherwise  = map ext rowRange
+  where
+    rowRange = [0..rows m-1]
+    sub k = subVector (k*xRow m) (cols m) (xdat m)
+    ext k = xdat $ unsafePerformIO $ extractR RowMajor m 1 (idxs[k]) 0 (idxs[0,cols m-1])
+
+{-
     | c == 0    = replicate r (fromList[])
     | otherwise = toRows' 0
   where
@@ -213,6 +213,7 @@ toRows m
     c = cols m
     toRows' k | k == r*c  = []
               | otherwise = subVector k c v : toRows' (k+c)
+-}
 
 -- | Creates a matrix from a list of vectors, as columns
 fromColumns :: Element t => [Vector t] -> Matrix t
@@ -240,7 +241,7 @@ matrixFromVector o r c v
     | r * c == dim v = m
     | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
   where
-    m | o == RowMajor = Matrix { irows = r, icols = c, xdat = v, xRow = c, xCol = 1    }
+    m | o == RowMajor = Matrix { irows = r, icols = c, xdat = v, xRow = c, xCol = 1 }
       | otherwise     = Matrix { irows = r, icols = c, xdat = v, xRow = 1, xCol = r }
 
 -- allocates memory for a new matrix
@@ -263,31 +264,26 @@ reshape :: Storable t => Int -> Vector t -> Matrix t
 reshape 0 v = matrixFromVector RowMajor 0 0 v
 reshape c v = matrixFromVector RowMajor (dim v `div` c) c v
 
---singleton x = reshape 1 (fromList [x])
 
 -- | application of a vector function on the flattened matrix elements
-liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
-liftMatrix f m@Matrix { irows = r, icols = c, xdat = d} = matrixFromVector (orderOf m) r c (f d)
+liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
+liftMatrix f m@Matrix { irows = r, icols = c, xdat = d}
+    | isSlice m = matrixFromVector RowMajor r c (f (flatten m))
+    | otherwise = matrixFromVector (orderOf m) r c (f d)
 
 -- | application of a vector function on the flattened matrices elements
 liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
-liftMatrix2 f m1 m2
-    | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
-    | otherwise = case orderOf m1 of
-        RowMajor    -> matrixFromVector RowMajor    (rows m1) (cols m1) (f (xdat m1) (flatten m2))
-        ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2))
-
-
-compat :: Matrix a -> Matrix b -> Bool
-compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
+liftMatrix2 f m1@(size->(r,c)) m2
+    | (r,c)/=size m2 = error "nonconformant matrices in liftMatrix2"
+    | rowOrder m1 = matrixFromVector RowMajor    r c (f (flatten m1) (flatten m2))
+    | otherwise   = matrixFromVector ColumnMajor r c (f (flatten (trans m1)) (flatten (trans m2)))
 
 ------------------------------------------------------------------
 
 -- | Supported matrix elements.
 class (Storable a) => Element a where
-    transdata :: Int -> Vector a -> Int -> Vector a
     constantD  :: a -> Int -> Vector a
-    extractR :: Matrix a -> CInt -> Vector CInt -> CInt -> Vector CInt -> IO (Matrix a)
+    extractR :: MatrixOrder -> Matrix a -> CInt -> Vector CInt -> CInt -> Vector CInt -> IO (Matrix a)
     setRect  :: Int -> Int -> Matrix a -> Matrix a -> IO ()
     sortI    :: Ord a => Vector a -> Vector CInt
     sortV    :: Ord a => Vector a -> Vector a
@@ -299,7 +295,6 @@ class (Storable a) => Element a where
 
 
 instance Element Float where
-    transdata  = transdataAux ctransF
     constantD  = constantAux cconstantF
     extractR   = extractAux c_extractF
     setRect    = setRectAux c_setRectF
@@ -312,7 +307,6 @@ instance Element Float where
     gemm       = gemmg c_gemmF
 
 instance Element Double where
-    transdata  = transdataAux ctransR
     constantD  = constantAux cconstantR
     extractR   = extractAux c_extractD
     setRect    = setRectAux c_setRectD
@@ -325,7 +319,6 @@ instance Element Double where
     gemm       = gemmg c_gemmD
 
 instance Element (Complex Float) where
-    transdata  = transdataAux ctransQ
     constantD  = constantAux cconstantQ
     extractR   = extractAux c_extractQ
     setRect    = setRectAux c_setRectQ
@@ -338,7 +331,6 @@ instance Element (Complex Float) where
     gemm       = gemmg c_gemmQ
 
 instance Element (Complex Double) where
-    transdata  = transdataAux ctransC
     constantD  = constantAux cconstantC
     extractR   = extractAux c_extractC
     setRect    = setRectAux c_setRectC
@@ -351,7 +343,6 @@ instance Element (Complex Double) where
     gemm       = gemmg c_gemmC
 
 instance Element (CInt) where
-    transdata  = transdataAux ctransI
     constantD  = constantAux cconstantI
     extractR   = extractAux c_extractI
     setRect    = setRectAux c_setRectI
@@ -364,7 +355,6 @@ instance Element (CInt) where
     gemm       = gemmg c_gemmI
 
 instance Element Z where
-    transdata  = transdataAux ctransL
     constantD  = constantAux cconstantL
     extractR   = extractAux c_extractL
     setRect    = setRectAux c_setRectL
@@ -378,32 +368,6 @@ instance Element Z where
 
 -------------------------------------------------------------------
 
-transdataAux fun c1 d c2 =
-    if noneed
-        then d
-        else unsafePerformIO $ do
-            -- putStrLn "T"
-            v <- createVector (dim d)
-            unsafeWith d $ \pd ->
-                unsafeWith v $ \pv ->
-                    fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
-            return v
-  where r1 = dim d `div` c1
-        r2 = dim d `div` c2
-        noneed = dim d == 0 || r1 == 1 || c1 == 1
-
-
-type TMM t = t ..> t ..> Ok
-
-foreign import ccall unsafe "transF" ctransF :: TMM Float
-foreign import ccall unsafe "transR" ctransR :: TMM Double
-foreign import ccall unsafe "transQ" ctransQ :: TMM (Complex Float)
-foreign import ccall unsafe "transC" ctransC :: TMM (Complex Double)
-foreign import ccall unsafe "transI" ctransI :: TMM CInt
-foreign import ccall unsafe "transL" ctransL :: TMM Z
-
-----------------------------------------------------------------------
-
 subMatrix :: Element a
           => (Int,Int) -- ^ (r0,c0) starting position
           -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
@@ -411,9 +375,8 @@ subMatrix :: Element a
           -> Matrix a -- ^ result
 subMatrix (r0,c0) (rt,ct) m
     | 0 <= r0 && 0 <= rt && r0+rt <= rows m &&
-      0 <= c0 && 0 <= ct && c0+ct <= cols m = unsafePerformIO $ extractR m 0 (idxs[r0,r0+rt-1]) 0 (idxs[c0,c0+ct-1])
-    | otherwise = error $ "wrong subMatrix "++
-                          show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
+      0 <= c0 && 0 <= ct && c0+ct <= cols m = unsafePerformIO $ extractR RowMajor m 0 (idxs[r0,r0+rt-1]) 0 (idxs[c0,c0+ct-1])
+    | otherwise = error $ "wrong subMatrix "++show ((r0,c0),(rt,ct))++" of "++shSize m
 
 
 sliceMatrix :: Element a
@@ -424,11 +387,12 @@ sliceMatrix :: Element a
 sliceMatrix (r0,c0) (rt,ct) m
     | 0 <= r0 && 0 <= rt && r0+rt <= rows m &&
       0 <= c0 && 0 <= ct && c0+ct <= cols m = res
-    | otherwise = error $ "wrong sliceMatrix "++
-                          show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
+    | otherwise = error $ "wrong sliceMatrix "++show ((r0,c0),(rt,ct))++" of "++shSize m
   where
-    t = r0 * xRow m + c0 * xCol m
-    res = m { irows = rt, icols = ct, xdat = subVector t (rt*ct) (xdat m) }
+    p = r0 * xRow m + c0 * xCol m
+    tot | rowOrder m = ct + (rt-1) * xRow m
+        | otherwise  = rt + (ct-1) * xCol m
+    res = m { irows = rt, icols = ct, xdat = subVector p tot (xdat m) }
 
 --------------------------------------------------------------------------
 
@@ -449,7 +413,7 @@ conformMTo (r,c) m
     | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
     | size m == (r,1) = repCols c m
     | size m == (1,c) = repRows r m
-    | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
+    | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c)
 
 conformVTo n v
     | dim v == n = v
@@ -459,9 +423,9 @@ conformVTo n v
 repRows n x = fromRows (replicate n (flatten x))
 repCols n x = fromColumns (replicate n (flatten x))
 
-size m = (rows m, cols m)
+shSize = shDim . size
 
-shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
+shDim (r,c) = "(" ++ show r ++"x"++ show c ++")"
 
 emptyM r c = matrixFromVector RowMajor r c (fromList[])
 
@@ -477,10 +441,10 @@ instance (Storable t, NFData t) => NFData (Matrix t)
 
 ---------------------------------------------------------------
 
-extractAux f m moder vr modec vc = do
+extractAux f ord m moder vr modec vc = do
     let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
         nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
-    r <- createMatrix RowMajor nr nc
+    r <- createMatrix ord nr nc
     f moder modec # vr # vc # m # r  #|"extract"
     return r
 
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 37f6e9b..c4f95d8 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -36,7 +36,7 @@ import Internal.LAPACK (multiplyI, multiplyL)
 import Internal.Algorithms(luFact)
 import Internal.Util(Normed(..),Indexable(..),
                      gaussElim, gaussElim_1, gaussElim_2,
-                     luST, luSolve', luPacked', magnit)
+                     luST, luSolve', luPacked', magnit, invershur)
 import Internal.ST(mutable)
 import GHC.TypeLits
 import Data.Proxy(Proxy)
@@ -126,9 +126,8 @@ instance forall n t . (Integral t, KnownNat n) => Num (Mod n t)
 
 instance KnownNat m => Element (Mod m I)
   where
-    transdata n v m = i2f (transdata n (f2i v) m)
     constantD x n = i2f (constantD (unMod x) n)
-    extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js
+    extractR ord m mi is mj js = i2fM <$> extractR ord (f2iM m) mi is mj js
     setRect i j m x = setRect i j (f2iM m) (f2iM x)
     sortI = sortI . f2i
     sortV = i2f . sortV . f2i
@@ -144,9 +143,8 @@ instance KnownNat m => Element (Mod m I)
 
 instance KnownNat m => Element (Mod m Z)
   where
-    transdata n v m = i2f (transdata n (f2i v) m)
     constantD x n = i2f (constantD (unMod x) n)
-    extractR m mi is mj js = i2fM <$> extractR (f2iM m) mi is mj js
+    extractR ord m mi is mj js = i2fM <$> extractR ord (f2iM m) mi is mj js
     setRect i j m x = setRect i j (f2iM m) (f2iM x)
     sortI = sortI . f2i
     sortV = i2f . sortV . f2i
@@ -293,11 +291,11 @@ f2i :: Storable t => Vector (Mod n t) -> Vector t
 f2i v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
     where (fp,i,n) = unsafeToForeignPtr v
 
-f2iM :: Storable t => Matrix (Mod n t) -> Matrix t
-f2iM = liftMatrix f2i
+f2iM :: (Element t, Element (Mod n t)) => Matrix (Mod n t) -> Matrix t
+f2iM m = m { xdat = f2i (xdat m) }
 
-i2fM :: Storable t => Matrix t -> Matrix (Mod n t)
-i2fM = liftMatrix i2f
+i2fM :: (Element t, Element (Mod n t)) => Matrix t -> Matrix (Mod n t)
+i2fM m = m { xdat = i2f (xdat m) }
 
 vmod :: forall m t. (KnownNat m, Storable t, Integral t, Numeric t) => Vector t -> Vector (Mod m t)
 vmod = i2f . cmod' m'
@@ -376,6 +374,8 @@ test = (ok, info)
          bb = flipud aa
          x = luSolve' (luPacked' aa) bb
 
+    tmm = diagRect 1 (fromList [2..6]) 5 5 :: Matrix (Mod 19 I)
+
     info = do
         print v
         print m
@@ -421,6 +421,9 @@ test = (ok, info)
         print $ checkSolve (sgen 5 :: Matrix (Complex Float))
         print $ checkSolve (gen 5 :: Matrix (Mod 7 I))
         print $ checkSolve (gen 5 :: Matrix (Mod 7 Z))
+        
+        print $ luSolve' (luPacked' tmm) (ident (rows tmm))
+        print $ invershur tmm
 
 
     ok = and
@@ -449,6 +452,7 @@ test = (ok, info)
       , 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)
+      , invershur tmm == luSolve' (luPacked' tmm) (ident (rows tmm))
       ]
 
 
diff --git a/packages/base/src/Internal/ST.hs b/packages/base/src/Internal/ST.hs
index 92654e4..73cdf0c 100644
--- a/packages/base/src/Internal/ST.hs
+++ b/packages/base/src/Internal/ST.hs
@@ -223,7 +223,7 @@ rowOper (SWAP i1 i2 r) (STMatrix m) = unsafeIOToST $ rowOp 2 0 i1' i2' j1 j2 m
     i2' = i2 `mod` (rows m)
 
 
-extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR m 0 (idxs[i1,i2]) 0 (idxs[j1,j2]))
+extractMatrix (STMatrix m) rr rc = unsafeIOToST (extractR (orderOf m) m 0 (idxs[i1,i2]) 0 (idxs[j1,j2]))
   where
     (i1,i2) = getRowRange (rows m) rr
     (j1,j2) = getColRange (cols m) rc
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index bf6c8b6..98eb4ef 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -54,7 +54,9 @@ module Internal.Util(
     -- ** 2D
     corr2, conv2, separable,
     block2x2,block3x3,view1,unView1,foldMatrix,
-    gaussElim_1, gaussElim_2, gaussElim, luST, luSolve', luSolve'', luPacked', luPacked''
+    gaussElim_1, gaussElim_2, gaussElim,
+    luST, luSolve', luSolve'', luPacked', luPacked'',
+    invershur
 ) where
 
 import Internal.Vector
@@ -829,6 +831,45 @@ luSolve' (lup,p) b = backSust lup (forwSust lup pb)
   where
     pb = b ?? (Pos (fixPerm' p), All)
 
+
+--------------------------------------------------------------------------------
+
+data MatrixView t b
+    = Elem t
+    | Block b b b b
+  deriving Show
+
+
+viewBlock' r c m
+    | (rt,ct) == (1,1) = Elem (atM' m 0 0)
+    | otherwise        = Block m11 m12 m21 m22
+  where
+    (rt,ct) = size m
+    m11 = sliceMatrix (0,0) (r,c)       m
+    m12 = sliceMatrix (0,c) (r,ct-c)    m
+    m21 = sliceMatrix (r,0) (rt-r,c)    m
+    m22 = sliceMatrix (r,c) (rt-r,ct-c) m
+
+viewBlock m = viewBlock' n n m
+  where
+    n = rows m `div` 2
+
+invershur (viewBlock -> Block a b c d) = fromBlocks [[a',b'],[c',d']]
+  where
+    r1 = invershur a
+    r2 = c <> r1
+    r3 = r1 <> b
+    r4 = c <> r3
+    r5 = r4-d
+    r6 = invershur r5
+    b' = r3 <> r6
+    c' = r6 <> r2
+    r7 = r3 <> c'
+    a' = r1-r7
+    d' = -r6
+
+invershur x = recip x
+
 --------------------------------------------------------------------------------
 
 instance Testable (Matrix I) where