11 files changed, 113 insertions, 193 deletions

diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index c4f1a60..7ead5b4 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -30,6 +30,7 @@ import Internal.LAPACK as LAPACK
 import Internal.Numeric
 import Data.List(foldl1')
 import qualified Data.Array as A
+import qualified Data.Vector.Storable as Vector
 import Internal.ST
 import Internal.Vectorized(range)
 import Control.DeepSeq
@@ -291,9 +292,13 @@ fromList [35.18264833189422,1.4769076999800903]
 
 -}
 compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
-compactSVD m = (u', subVector 0 d s, v') where
+compactSVD = compactSVDTol 1
+
+-- | @compactSVDTol r@ is similar to 'compactSVD', but uses tolerance @tol=r*g*eps*(max rows cols)@ to distinguish nonzero singular values, where @g@ is the greatest singular value.
+compactSVDTol :: Field t  => Double -> Matrix t -> (Matrix t, Vector Double, Matrix t)
+compactSVDTol r m = (u', subVector 0 d s, v') where
     (u,s,v) = thinSVD m
-    d = rankSVD (1*eps) m s `max` 1
+    d = rankSVD (r*eps) m s `max` 1
     u' = takeColumns d u
     v' = takeColumns d v
 
@@ -475,9 +480,14 @@ instance (NFData t, Numeric t) => NFData (QR t)
 -- | QR factorization.
 --
 -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
+-- Note: the current implementation is very slow for large matrices. 'thinQR' is much faster.
 qr :: Field t => Matrix t -> (Matrix t, Matrix t)
 qr = {-# SCC "qr" #-} unpackQR . qr'
 
+-- | A version of 'qr' which returns only the @min (rows m) (cols m)@ columns of @q@ and rows of @r@.
+thinQR :: Field t => Matrix t -> (Matrix t, Matrix t)
+thinQR = {-# SCC "thinQR" #-} thinUnpackQR . qr'
+
 -- | Compute the QR decomposition of a matrix in compact form.
 qrRaw :: Field t => Matrix t -> QR t
 qrRaw m = QR x v
@@ -494,9 +504,17 @@ qrgr n (QR a t)
 -- | RQ factorization.
 --
 -- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
+-- Note: the current implementation is very slow for large matrices. 'thinRQ' is much faster.
 rq :: Field t => Matrix t -> (Matrix t, Matrix t)
-rq m =  {-# SCC "rq" #-} (r,q) where
-    (q',r') = qr $ trans $ rev1 m
+rq = {-# SCC "rq" #-} rqFromQR qr
+
+-- | A version of 'rq' which returns only the @min (rows m) (cols m)@ columns of @r@ and rows of @q@.
+thinRQ :: Field t => Matrix t -> (Matrix t, Matrix t)
+thinRQ = {-# SCC "thinQR" #-} rqFromQR thinQR
+
+rqFromQR :: Field t => (Matrix t -> (Matrix t, Matrix t)) -> Matrix t -> (Matrix t, Matrix t)
+rqFromQR qr0 m = (r,q) where
+    (q',r') = qr0 $ trans $ rev1 m
     r = rev2 (trans r')
     q = rev2 (trans q')
     rev1 = flipud . fliprl
@@ -724,6 +742,12 @@ unpackQR (pq, tau) =  {-# SCC "unpackQR" #-} (q,r)
           hs = zipWith haussholder (toList tau) vs
           q = foldl1' mXm hs
 
+thinUnpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
+thinUnpackQR (pq, tau) = (q, r)
+  where mn = uncurry min $ size pq
+        q = qrgr mn $ QR pq tau
+        r = fromRows $ zipWith (\i v -> Vector.replicate i 0 Vector.++ Vector.drop i v) [0..mn-1] (toRows pq)
+
 unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
 unpackHess hf m
     | rows m == 1 = ((1><1)[1],m)
diff --git a/packages/base/src/Internal/Devel.hs b/packages/base/src/Internal/Devel.hs
index 92b5604..3887663 100644
--- a/packages/base/src/Internal/Devel.hs
+++ b/packages/base/src/Internal/Devel.hs
@@ -80,8 +80,8 @@ class TransArray c
   where
     type Trans c b
     type TransRaw c b
-    apply      :: (Trans c b) -> c -> b
-    applyRaw   :: (TransRaw c b) -> c -> b
+    apply      :: c -> (b -> IO r) -> (Trans c b) -> IO r
+    applyRaw   :: c -> (b -> IO r) -> (TransRaw c b) -> IO r
     infixl 1 `apply`, `applyRaw`
 
 instance Storable t => TransArray (Vector t)
@@ -92,4 +92,3 @@ instance Storable t => TransArray (Vector t)
     {-# INLINE apply #-}
     applyRaw = avec
     {-# INLINE applyRaw #-}
-
diff --git a/packages/base/src/Internal/Foreign.hs b/packages/base/src/Internal/Foreign.hs
deleted file mode 100644
index ea071a4..0000000
--- a/packages/base/src/Internal/Foreign.hs
+++ /dev/null
@@ -1,102 +0,0 @@
-{-# LANGUAGE MagicHash, UnboxedTuples #-}
--- | FFI and hmatrix helpers.
---
--- Sample usage, to upload a perspective matrix to a shader.
---
--- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) 
--- @
---
-
-module Internal.Foreign 
-    ( app
-    , appVector, appVectorLen
-    , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
-    , unsafeMatrixToVector, unsafeMatrixToForeignPtr
-    ) where
-
-import Foreign.C.Types(CInt)
-import Internal.Vector
-import Internal.Matrix
-import qualified Data.Vector.Storable as S
-import Foreign (Ptr, ForeignPtr, Storable)
-import GHC.Base (IO(..), realWorld#)
-
-{-# INLINE unsafeInlinePerformIO #-}
--- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be
--- unecessary calls to unsafePerformIO or its internals.
-unsafeInlinePerformIO :: IO a -> a
-unsafeInlinePerformIO (IO f) = case f realWorld# of
-    (# _, x #) -> x
-
-{-# INLINE app #-}
--- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
--- e.g.
---
--- @
--- someFunction
---     \`appMatrixLen\` m
---     \`appVectorLen\` v
---     \`app\` other
---     \`app\` arguments
---     \`app\` go here
--- @
---
--- One could also write:
---
--- @
--- (someFunction 
---     \`appMatrixLen\` m
---     \`appVectorLen\` v) 
---     other 
---     arguments 
---     (go here)
--- @
---
-app :: (a -> b) -> a -> b
-app f = f
-
-{-# INLINE appVector #-}
-appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
-appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
-
-{-# INLINE appVectorLen #-}
-appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
-appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
-
-{-# INLINE appMatrix #-}
-appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
-appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
-
-{-# INLINE appMatrixLen #-}
-appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
-  where
-    r = fromIntegral (rows x)
-    c = fromIntegral (cols x)
-
-{-# INLINE appMatrixRaw #-}
-appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
-appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
-
-{-# INLINE appMatrixRawLen #-}
-appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
-  where
-    r = fromIntegral (rows x)
-    c = fromIntegral (cols x)
-
-infixl 1 `app`
-infixl 1 `appVector`
-infixl 1 `appMatrix`
-infixl 1 `appMatrixRaw`
-
-{-# INLINE unsafeMatrixToVector #-}
--- | This will disregard the order of the matrix, and simply return it as-is. 
--- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
-unsafeMatrixToVector :: Matrix a -> Vector a
-unsafeMatrixToVector = xdat
-
-{-# INLINE unsafeMatrixToForeignPtr #-}
-unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
-unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
-
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index c2c140b..231109a 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -18,7 +18,7 @@ module Internal.LAPACK where
 
 import Internal.Devel
 import Internal.Vector
-import Internal.Matrix hiding ((#))
+import Internal.Matrix hiding ((#), (#!))
 import Internal.Conversion
 import Internal.Element
 import Foreign.Ptr(nullPtr)
@@ -28,10 +28,13 @@ import System.IO.Unsafe(unsafePerformIO)
 
 -----------------------------------------------------------------------------------
 
-infixl 1 #
+infixr 1 #
 a # b = apply a b
 {-# INLINE (#) #-}
 
+a #! b = a # b # id
+{-# INLINE (#!) #-}
+
 -----------------------------------------------------------------------------------
 
 type TMMM t = t ::> t ::> t ::> Ok
@@ -56,7 +59,7 @@ multiplyAux f st a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++
                                        show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    f (isT a) (isT b) # (tt a) # (tt b) # s #| st
+    ((tt a) # (tt b) #! s) (f (isT a) (isT b)) #| st
     return s
 
 -- | Matrix product based on BLAS's /dgemm/.
@@ -80,7 +83,7 @@ multiplyI m a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $
         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    c_multiplyI m # a # b # s #|"c_multiplyI"
+    (a # b #! s) (c_multiplyI m) #|"c_multiplyI"
     return s
 
 multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
@@ -88,7 +91,7 @@ multiplyL m a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $
         "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b
     s <- createMatrix ColumnMajor (rows a) (cols b)
-    c_multiplyL m # a # b # s #|"c_multiplyL"
+    (a # b #! s) (c_multiplyL m) #|"c_multiplyL"
     return s
 
 -----------------------------------------------------------------------------
@@ -121,7 +124,7 @@ svdAux f st x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
-    f # a # u # s # v #| st
+    (a # u # s #! v) f #| st
     return (u,s,v)
   where
     r = rows x
@@ -149,7 +152,7 @@ thinSVDAux f st x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r q
     s <- createVector q
     v <- createMatrix ColumnMajor q c
-    f # a # u # s # v #| st
+    (a # u # s #! v) f #| st
     return (u,s,v)
   where
     r = rows x
@@ -176,7 +179,7 @@ svCd = svAux zgesdd "svCd"
 svAux f st x = unsafePerformIO $ do
     a <- copy ColumnMajor x
     s <- createVector q
-    g # a # s #| st
+    (a #! s) g #| st
     return s
   where
     r = rows x
@@ -197,7 +200,7 @@ rightSVAux f st x = unsafePerformIO $ do
     a <- copy ColumnMajor x
     s <- createVector q
     v <- createMatrix ColumnMajor c c
-    g # a # s # v #| st
+    (a # s #! v) g #| st
     return (s,v)
   where
     r = rows x
@@ -218,7 +221,7 @@ leftSVAux f st x = unsafePerformIO $ do
     a <- copy ColumnMajor x
     u <- createMatrix ColumnMajor r r
     s <- createVector q
-    g # a # u # s #| st
+    (a # u #! s) g #| st
     return (u,s)
   where
     r = rows x
@@ -237,7 +240,7 @@ eigAux f st m = unsafePerformIO $ do
     a <- copy ColumnMajor m
     l <- createVector r
     v <- createMatrix ColumnMajor r r
-    g # a # l # v #| st
+    (a # l #! v) g #| st
     return (l,v)
   where
     r = rows m
@@ -252,7 +255,7 @@ eigC = eigAux zgeev "eigC"
 eigOnlyAux f st m = unsafePerformIO $ do
     a <- copy ColumnMajor m
     l <- createVector r
-    g # a # l #| st
+    (a #! l) g #| st
     return l
   where
     r = rows m
@@ -277,7 +280,7 @@ eigRaux m = unsafePerformIO $ do
     a <- copy ColumnMajor m
     l <- createVector r
     v <- createMatrix ColumnMajor r r
-    g # a # l # v #| "eigR"
+    (a # l #! v) g #| "eigR"
     return (l,v)
   where
     r = rows m
@@ -305,7 +308,7 @@ eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR"
 eigSHAux f st m = unsafePerformIO $ do
     l <- createVector r
     v <- copy ColumnMajor m
-    f # l # v #| st
+    (l #! v) f #| st
     return (l,v)
   where
     r = rows m
@@ -356,7 +359,7 @@ linearSolveSQAux g f st a b
     | n1==n2 && n1==r = unsafePerformIO . g $ do
         a' <- copy ColumnMajor a
         s  <- copy ColumnMajor b
-        f # a' # s #| st
+        (a' #! s) f #| st
         return s
     | otherwise = error $ st ++ " of nonsquare matrix"
   where
@@ -387,7 +390,7 @@ foreign import ccall unsafe "cholSolveC_l" zpotrs  :: C ::> C ::> Ok
 linearSolveSQAux2 g f st a b
     | n1==n2 && n1==r = unsafePerformIO . g $ do
         s <- copy ColumnMajor b
-        f # a # s #| st
+        (a #! s) f #| st
         return s
     | otherwise = error $ st ++ " of nonsquare matrix"
   where
@@ -415,7 +418,7 @@ linearSolveAux f st a b
         a' <- copy ColumnMajor a
         r  <- createMatrix ColumnMajor (max m n) nrhs
         setRect 0 0 b r
-        f # a' # r #| st
+        (a' #! r) f #| st
         return r
     | otherwise = error $ "different number of rows in linearSolve ("++st++")"
   where
@@ -458,7 +461,7 @@ foreign import ccall unsafe "chol_l_S" dpotrf :: R ::> Ok
 
 cholAux f st a = do
     r <- copy ColumnMajor a
-    f # r #| st
+    (r # id) f #| st
     return r
 
 -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
@@ -495,7 +498,7 @@ qrC = qrAux zgeqr2 "qrC"
 qrAux f st a = unsafePerformIO $ do
     r <- copy ColumnMajor a
     tau <- createVector mn
-    f # tau # r #| st
+    (tau #! r) f #| st
     return (r,tau)
   where
     m = rows a
@@ -514,7 +517,7 @@ qrgrC = qrgrAux zungqr "qrgrC"
 
 qrgrAux f st n (a, tau) = unsafePerformIO $ do
     res <- copy ColumnMajor (subMatrix (0,0) (rows a,n) a)
-    f # (subVector 0 n tau') # res #| st
+    ((subVector 0 n tau') #! res) f #| st
     return res
   where
     tau' = vjoin [tau, constantD 0 n]
@@ -534,7 +537,7 @@ hessC = hessAux zgehrd "hessC"
 hessAux f st a = unsafePerformIO $ do
     r <- copy ColumnMajor a
     tau <- createVector (mn-1)
-    f # tau # r #| st
+    (tau #! r) f #| st
     return (r,tau)
   where
     m = rows a
@@ -556,7 +559,7 @@ schurC = schurAux zgees "schurC"
 schurAux f st a = unsafePerformIO $ do
     u <- createMatrix ColumnMajor n n
     s <- copy ColumnMajor a
-    f # u # s #| st
+    (u #! s) f #| st
     return (u,s)
   where
     n = rows a
@@ -576,7 +579,7 @@ luC = luAux zgetrf "luC"
 luAux f st a = unsafePerformIO $ do
     lu <- copy ColumnMajor a
     piv <- createVector (min n m)
-    f # piv # lu #| st
+    (piv #! lu) f #| st
     return (lu, map (pred.round) (toList piv))
   where
     n = rows a
@@ -598,7 +601,7 @@ lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv b
 lusAux f st a piv b
     | n1==n2 && n2==n =unsafePerformIO $ do
          x <- copy ColumnMajor b
-         f # a # piv' # x #| st
+         (a # piv' #! x) f #| st
          return x
     | otherwise = error st
   where
@@ -622,7 +625,7 @@ ldlC = ldlAux zhetrf "ldlC"
 ldlAux f st a = unsafePerformIO $ do
     ldl <- copy ColumnMajor a
     piv <- createVector (rows a)
-    f # piv # ldl #| st
+    (piv #! ldl) f #| st
     return (ldl, map (pred.round) (toList piv))
 
 -----------------------------------------------------------------------------------
@@ -637,4 +640,3 @@ ldlsR a piv b = lusAux dsytrs "ldlsR" (fmat a) piv b
 -- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/.
 ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
 ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b
-
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index c47c625..0135288 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -22,7 +22,7 @@ module Internal.Matrix where
 
 import Internal.Vector
 import Internal.Devel
-import Internal.Vectorized hiding ((#))
+import Internal.Vectorized hiding ((#), (#!))
 import Foreign.Marshal.Alloc ( free )
 import Foreign.Marshal.Array(newArray)
 import Foreign.Ptr ( Ptr )
@@ -110,15 +110,15 @@ fmat m
 
 -- C-Haskell matrix adapters
 {-# INLINE amatr #-}
-amatr :: Storable a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-amatr f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c))
+amatr :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> Ptr a -> f) -> IO r
+amatr x f g = unsafeWith (xdat x) (f . g r c)
   where
     r  = fi (rows x)
     c  = fi (cols x)
 
 {-# INLINE amat #-}
-amat :: Storable a => (CInt -> CInt -> CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
-amat f x = inlinePerformIO (unsafeWith (xdat x) (return . f r c sr sc))
+amat :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> CInt -> CInt -> Ptr a -> f) -> IO r
+amat x f g = unsafeWith (xdat x) (f . g r c sr sc)
   where
     r  = fi (rows x)
     c  = fi (cols x)
@@ -135,10 +135,13 @@ instance Storable t => TransArray (Matrix t)
     applyRaw = amatr
     {-# INLINE applyRaw #-}
 
-infixl 1 #
+infixr 1 #
 a # b = apply a b
 {-# INLINE (#) #-}
 
+a #! b = a # b # id
+{-# INLINE (#!) #-}
+
 --------------------------------------------------------------------------------
 
 copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])
@@ -426,7 +429,8 @@ 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 ord nr nc
-    f moder modec # vr # vc # m # r  #|"extract"
+    (vr # vc # m #! r) (f moder modec)  #|"extract"
+
     return r
 
 type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))
@@ -440,7 +444,7 @@ foreign import ccall unsafe "extractL" c_extractL :: Extr Z
 
 ---------------------------------------------------------------
 
-setRectAux f i j m r = f (fi i) (fi j) # m # r #|"setRect"
+setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"
 
 type SetRect x = I -> I -> x ::> x::> Ok
 
@@ -455,7 +459,7 @@ foreign import ccall unsafe "setRectL" c_setRectL :: SetRect Z
 
 sortG f v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # v # r #|"sortG"
+    (v #! r) f #|"sortG"
     return r
 
 sortIdxD = sortG c_sort_indexD
@@ -482,7 +486,7 @@ foreign import ccall unsafe "sort_valuesL" c_sort_valL :: Z :> Z :> Ok
 
 compareG f u v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # u # v # r #|"compareG"
+    (u # v #! r) f #|"compareG"
     return r
 
 compareD = compareG c_compareD
@@ -499,7 +503,7 @@ foreign import ccall unsafe "compareL" c_compareL :: Z :> Z :> I :> Ok
 
 selectG f c u v w = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # c # u # v # w # r #|"selectG"
+    (c # u # v # w #! r) f #|"selectG"
     return r
 
 selectD = selectG c_selectD
@@ -522,7 +526,7 @@ foreign import ccall unsafe "chooseL" c_selectL :: Sel Z
 
 remapG f i j m = unsafePerformIO $ do
     r <- createMatrix RowMajor (rows i) (cols i)
-    f # i # j # m # r #|"remapG"
+    (i # j # m #! r) f #|"remapG"
     return r
 
 remapD = remapG c_remapD
@@ -545,7 +549,7 @@ foreign import ccall unsafe "remapL" c_remapL :: Rem Z
 
 rowOpAux f c x i1 i2 j1 j2 m = do
     px <- newArray [x]
-    f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2) # m #|"rowOp"
+    (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp"
     free px
 
 type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok
@@ -561,7 +565,7 @@ foreign import ccall unsafe "rowop_mod_int64_t" c_rowOpML :: Z -> RowOp Z
 
 --------------------------------------------------------------------------------
 
-gemmg f v m1 m2 m3 = f # v # m1 # m2 # m3 #|"gemmg"
+gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"
 
 type Tgemm x = x :> x ::> x ::> x ::> Ok
 
@@ -589,10 +593,9 @@ saveMatrix
 saveMatrix name format m = do
     cname   <- newCString name
     cformat <- newCString format
-    c_saveMatrix cname cformat # m #|"saveMatrix"
+    (m # id) (c_saveMatrix cname cformat) #|"saveMatrix"
     free cname
     free cformat
     return ()
 
 --------------------------------------------------------------------------------
-
diff --git a/packages/base/src/Internal/Sparse.hs b/packages/base/src/Internal/Sparse.hs
index 1604e7e..1ff3f57 100644
--- a/packages/base/src/Internal/Sparse.hs
+++ b/packages/base/src/Internal/Sparse.hs
@@ -144,13 +144,13 @@ gmXv :: GMatrix -> Vector Double -> Vector Double
 gmXv SparseR { gmCSR = CSR{..}, .. } v = unsafePerformIO $ do
     dim v /= nCols ~!~ printf "gmXv (CSR): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
     r <- createVector nRows
-    c_smXv # csrVals # csrCols # csrRows # v # r #|"CSRXv"
+    (csrVals # csrCols # csrRows # v #! r) c_smXv #|"CSRXv"
     return r
 
 gmXv SparseC { gmCSC = CSC{..}, .. } v = unsafePerformIO $ do
     dim v /= nCols ~!~ printf "gmXv (CSC): incorrect sizes: (%d,%d) x %d" nRows nCols (dim v)
     r <- createVector nRows
-    c_smTXv # cscVals # cscRows # cscCols # v # r #|"CSCXv"
+    (cscVals # cscRows # cscCols # v #! r) c_smTXv #|"CSCXv"
     return r
 
 gmXv Diag{..} v
@@ -211,4 +211,3 @@ instance Transposable GMatrix GMatrix
     tr (Diag v n m) = Diag v m n
     tr (Dense a n m) = Dense (tr a) m n
     tr' = tr
-
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index cf42961..7f52a9c 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -255,6 +255,7 @@ norm :: Vector Double -> Double
 -- ^ 2-norm of real vector
 norm = pnorm PNorm2
 
+-- | p-norm for vectors, operator norm for matrices
 class Normed a
   where
     norm_0   :: a -> R
@@ -319,10 +320,11 @@ instance Normed (Vector (Complex Float))
     norm_2 = norm_2 . double
     norm_Inf = norm_Inf . double
 
-
+-- | Frobenius norm (Schatten p-norm with p=2)
 norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
 norm_Frob = norm_2 . flatten
 
+-- | Sum of singular values (Schatten p-norm with p=1)
 norm_nuclear :: Field t => Matrix t -> R
 norm_nuclear = sumElements . singularValues
 
diff --git a/packages/base/src/Internal/Vector.hs b/packages/base/src/Internal/Vector.hs
index b4e235c..c4a310d 100644
--- a/packages/base/src/Internal/Vector.hs
+++ b/packages/base/src/Internal/Vector.hs
@@ -66,9 +66,8 @@ dim = Vector.length
 
 -- C-Haskell vector adapter
 {-# INLINE avec #-}
-avec :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
-avec f v = inlinePerformIO (unsafeWith v (return . f (fromIntegral (Vector.length v))))
-infixl 1 `avec`
+avec :: Storable a => Vector a -> (f -> IO r) -> ((CInt -> Ptr a -> f) -> IO r)
+avec v f g = unsafeWith v $ \ptr -> f (g (fromIntegral (Vector.length v)) ptr)
 
 -- allocates memory for a new vector
 createVector :: Storable a => Int -> IO (Vector a)
@@ -199,7 +198,7 @@ takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (
 
 ---------------------------------------------------------------
 
--- | transforms a complex vector into a real vector with alternating real and imaginary parts 
+-- | transforms a complex vector into a real vector with alternating real and imaginary parts
 asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
 asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
     where (fp,i,n) = unsafeToForeignPtr v
@@ -244,7 +243,7 @@ zipVectorWith f u v = unsafePerformIO $ do
 {-# INLINE zipVectorWith #-}
 
 -- | unzipWith for Vectors
-unzipVectorWith :: (Storable (a,b), Storable c, Storable d) 
+unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
                    => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
 unzipVectorWith f u = unsafePerformIO $ do
       let n = dim u
@@ -255,7 +254,7 @@ unzipVectorWith f u = unsafePerformIO $ do
               unsafeWith w $ \pw -> do
                   let go (-1) = return ()
                       go !k   = do z <- peekElemOff pu k
-                                   let (x,y) = f z 
+                                   let (x,y) = f z
                                    pokeElemOff      pv k x
                                    pokeElemOff      pw k y
                                    go (k-1)
@@ -303,11 +302,11 @@ mapVectorM f v = do
     return w
     where mapVectorM' w' !k !t
               | k == t               = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f x
                                        return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
               | otherwise            = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f x
                                        _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
                                        mapVectorM' w' (k+1) t
@@ -322,7 +321,7 @@ mapVectorM_ f v = do
                                     x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     f x
               | otherwise         = do
-                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     _ <- f x
                                     mapVectorM' (k+1) t
 {-# INLINE mapVectorM_ #-}
@@ -336,11 +335,11 @@ mapVectorWithIndexM f v = do
     return w
     where mapVectorM' w' !k !t
               | k == t               = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f k x
                                        return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
               | otherwise            = do
-                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                        y <- f k x
                                        _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
                                        mapVectorM' w' (k+1) t
@@ -355,7 +354,7 @@ mapVectorWithIndexM_ f v = do
                                     x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     f k x
               | otherwise         = do
-                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k 
+                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                     _ <- f k x
                                     mapVectorM' (k+1) t
 {-# INLINE mapVectorWithIndexM_ #-}
@@ -454,4 +453,3 @@ unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vect
 unzipVector = unzipVectorWith id
 
 -------------------------------------------------------------------
-
diff --git a/packages/base/src/Internal/Vectorized.hs b/packages/base/src/Internal/Vectorized.hs
index 03bcf90..a410bb2 100644
--- a/packages/base/src/Internal/Vectorized.hs
+++ b/packages/base/src/Internal/Vectorized.hs
@@ -27,10 +27,13 @@ import Foreign.C.String
 import System.IO.Unsafe(unsafePerformIO)
 import Control.Monad(when)
 
-infixl 1 #
+infixr 1 #
 a # b = applyRaw a b
 {-# INLINE (#) #-}
 
+a #! b = a # b # id
+{-# INLINE (#!) #-}
+
 fromei x = fromIntegral (fromEnum x) :: CInt
 
 data FunCodeV = Sin
@@ -103,7 +106,7 @@ sumL m = sumg (c_sumL m)
 
 sumg f x = unsafePerformIO $ do
     r <- createVector 1
-    f # x # r #| "sum"
+    (x #! r) f #| "sum"
     return $ r @> 0
 
 type TVV t = t :> t :> Ok
@@ -139,7 +142,7 @@ prodL = prodg . c_prodL
 
 prodg f x = unsafePerformIO $ do
     r <- createVector 1
-    f # x # r #| "prod"
+    (x #! r) f #| "prod"
     return $ r @> 0
 
 
@@ -154,24 +157,24 @@ foreign import ccall unsafe "prodL" c_prodL :: Z -> TVV Z
 
 toScalarAux fun code v = unsafePerformIO $ do
     r <- createVector 1
-    fun (fromei code) # v # r #|"toScalarAux"
+    (v #! r) (fun (fromei code)) #|"toScalarAux"
     return (r @> 0)
 
 vectorMapAux fun code v = unsafePerformIO $ do
     r <- createVector (dim v)
-    fun (fromei code) # v # r #|"vectorMapAux"
+    (v #! r) (fun (fromei code)) #|"vectorMapAux"
     return r
 
 vectorMapValAux fun code val v = unsafePerformIO $ do
     r <- createVector (dim v)
     pval <- newArray [val]
-    fun (fromei code) pval # v # r #|"vectorMapValAux"
+    (v #! r) (fun (fromei code) pval) #|"vectorMapValAux"
     free pval
     return r
 
 vectorZipAux fun code u v = unsafePerformIO $ do
     r <- createVector (dim u)
-    fun (fromei code) # u # v # r #|"vectorZipAux"
+    (u # v #! r) (fun (fromei code)) #|"vectorZipAux"
     return r
 
 ---------------------------------------------------------------------
@@ -368,7 +371,7 @@ randomVector :: Seed
              -> Vector Double
 randomVector seed dist n = unsafePerformIO $ do
     r <- createVector n
-    c_random_vector (fi seed) ((fi.fromEnum) dist) # r #|"randomVector"
+    (r # id) (c_random_vector (fi seed) ((fi.fromEnum) dist)) #|"randomVector"
     return r
 
 foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> Double :> Ok
@@ -377,7 +380,7 @@ foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> D
 
 roundVector v = unsafePerformIO $ do
     r <- createVector (dim v)
-    c_round_vector # v # r #|"roundVector"
+    (v #! r) c_round_vector #|"roundVector"
     return r
 
 foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
@@ -391,7 +394,7 @@ foreign import ccall unsafe "round_vector" c_round_vector :: TVV Double
 range :: Int -> Vector I
 range n = unsafePerformIO $ do
     r <- createVector n
-    c_range_vector # r #|"range"
+    (r # id) c_range_vector #|"range"
     return r
 
 foreign import ccall unsafe "range_vector" c_range_vector :: CInt :> Ok
@@ -431,7 +434,7 @@ long2intV = tog c_long2int
 
 tog f v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # v # r #|"tog"
+    (v #! r) f #|"tog"
     return r
 
 foreign import ccall unsafe "float2double" c_float2double :: Float :> Double :> Ok
@@ -450,7 +453,7 @@ foreign import ccall unsafe "long2int"    c_long2int    :: Z :> I :> Ok
 
 stepg f v = unsafePerformIO $ do
     r <- createVector (dim v)
-    f # v # r #|"step"
+    (v #! r) f #|"step"
     return r
 
 stepD :: Vector Double -> Vector Double
@@ -475,7 +478,7 @@ foreign import ccall unsafe "stepL" c_stepL :: TVV Z
 
 conjugateAux fun x = unsafePerformIO $ do
     v <- createVector (dim x)
-    fun # x # v #|"conjugateAux"
+    (x #! v) fun #|"conjugateAux"
     return v
 
 conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
@@ -493,7 +496,7 @@ cloneVector v = do
         let n = dim v
         r <- createVector n
         let f _ s _ d =  copyArray d s n >> return 0
-        f # v # r #|"cloneVector"
+        (v #! r) f #|"cloneVector"
         return r
 
 --------------------------------------------------------------------------------
@@ -501,7 +504,7 @@ cloneVector v = do
 constantAux fun x n = unsafePerformIO $ do
     v <- createVector n
     px <- newArray [x]
-    fun px # v #|"constantAux"
+    (v # id) (fun px) #|"constantAux"
     free px
     return v
 
@@ -515,4 +518,3 @@ foreign import ccall unsafe "constantI" cconstantI :: TConst CInt
 foreign import ccall unsafe "constantL" cconstantL :: TConst Z
 
 ----------------------------------------------------------------------
-
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index b5bacf4..626ee52 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -117,6 +117,7 @@ module Numeric.LinearAlgebra (
     svd,
     thinSVD,
     compactSVD,
+    compactSVDTol,
     singularValues,
     leftSV, rightSV,
 
@@ -126,7 +127,7 @@ module Numeric.LinearAlgebra (
     geigSH,
 
     -- * QR
-    qr, rq, qrRaw, qrgr,
+    qr, thinQR, rq, thinRQ, qrRaw, qrgr,
 
     -- * Cholesky
     chol, mbChol,
diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
index 57a68e7..941b597 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs
@@ -12,13 +12,6 @@ The library can be easily extended using the tools in this module.
 --------------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Devel(
-    -- * FFI helpers
-    -- | Sample usage, to upload a perspective matrix to a shader.
-    --
-    -- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
-    -- @
-    module Internal.Foreign,
-
     -- * FFI tools
     -- | See @examples/devel@ in the repository.
     
@@ -66,7 +59,6 @@ module Numeric.LinearAlgebra.Devel(
 
 ) where
 
-import Internal.Foreign
 import Internal.Devel
 import Internal.ST
 import Internal.Vector