simplified (but wrong)

7 files changed, 188 insertions, 69 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 48652f3..6ba2d06 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -27,6 +27,10 @@ import Data.Maybe(fromJust)
 
 ----------------------------------------------------------------
 
+-- the condition Storable a => Field a means that we can only put
+-- in Field types that are in Storable, and therefore Storable a
+-- is not required in signatures if we have a Field a.
+
 class Storable a => Field a where
     constant :: a -> Int -> Vector a
     transdata :: Int -> Vector a -> Int -> Vector a
@@ -36,7 +40,6 @@ class Storable a => Field a where
               -> Matrix a -> Matrix a
     diag :: Vector a -> Matrix a
 
-
 instance Field Double where
     constant  = constantR
     transdata = transdataR
@@ -78,12 +81,40 @@ foreign import ccall safe "aux.h transC"
 
 transdataG c1 d c2 = fromList . concat . transpose . partit c1 . toList $ d
 
+{- 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. Therefore, maybe it is better to have something like
+  viewC and viewF, which may actually perform a transpose if required.
 
+- 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. Access
+  to such pointer cannot be made by creating a new vector.
+
+-}
 
 data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
 
+{-
+
+
+
 data Matrix t = M { rows    :: Int
                   , cols    :: Int
                   , dat     :: Vector t
@@ -91,28 +122,26 @@ data Matrix t = M { rows    :: Int
                   , isTrans :: Bool
                   , order   :: MatrixOrder
                   } -- deriving Typeable
+-}
 
+data Matrix t = MC { rows :: Int, cols :: Int, dat :: Vector t } -- row major order
+              | MF { rows :: Int, cols :: Int, dat :: Vector t } -- column major order
 
-data NMat t = MC { rws, cls :: Int, dtc :: Vector t}
-            | MF { rws, cls :: Int, dtf :: Vector t}
-            | Tr (NMat t)
-
-ntrans (Tr m) = m
-ntrans m = Tr m
+-- transposition just changes the data order
+trans :: Matrix t -> Matrix t
+trans MC {rows = r, cols = c, dat = d} = MF {rows = c, cols = r, dat = d}
+trans MF {rows = r, cols = c, dat = d} = MC {rows = c, cols = r, dat = d}
 
 viewC m@MC{} = m
-viewF m@MF{} = m
+viewC MF {rows = r, cols = c, dat = d} = MC {rows = r, cols = c, dat = transdata r d c}
 
-fortran m = order m == ColumnMajor
+viewF m@MF{} = m
+viewF MC {rows = r, cols = c, dat = d} = MF {rows = r, cols = c, dat = transdata c d r}
 
-cdat m = if fortran m `xor` isTrans m then tdat m else dat m
-fdat m = if fortran m `xor` isTrans m then dat m else tdat m
+--fortran m = order m == ColumnMajor
 
-trans :: Matrix t -> Matrix t
-trans m = m { rows = cols m
-            , cols = rows m
-            , isTrans = not (isTrans m)
-            }
+cdat m = dat (viewC m)
+fdat m = dat (viewF m)
 
 type Mt t s = Int -> Int -> Ptr t -> s
 -- not yet admitted by my haddock version
@@ -120,11 +149,14 @@ type Mt t s = Int -> Int -> Ptr t -> s
 -- type t ::> s = Mt t s
 
 mat d m f = f (rows m) (cols m) (ptr (d m))
+--mat m f = f (rows m) (cols m) (ptr (dat m))
+--matC m f = f (rows m) (cols m) (ptr (cdat m))
+
 
-toLists :: (Storable t) => Matrix t -> [[t]]
+--toLists :: (Storable t) => Matrix t -> [[t]]
 toLists m = partit (cols m) . toList . cdat $ m
 
-instance (Show a, Storable a) => (Show (Matrix a)) where
+instance (Show a, Field a) => (Show (Matrix a)) where
     show m = (sizes++) . dsp . map (map show) . toLists $ m
         where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"
 
@@ -136,6 +168,7 @@ dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unw
         pad n str = replicate (n - length str) ' ' ++ str
         unwords' = concat . intersperse ", "
 
+{-
 matrixFromVector RowMajor c v =
     M { rows = r
       , cols = c
@@ -147,8 +180,6 @@ matrixFromVector RowMajor c v =
               r | m==0 = d
                 | otherwise = error "matrixFromVector"
 
--- r = dim v `div` c -- TODO check mod=0
-
 matrixFromVector ColumnMajor c v =
     M { rows = r
       , cols = c
@@ -160,6 +191,23 @@ matrixFromVector ColumnMajor c v =
               r | m==0 = d
                 | otherwise = error "matrixFromVector"
 
+-}
+
+matrixFromVector RowMajor c v = MC { rows = r, cols = c, dat = v}
+    where (d,m) = dim v `divMod` c
+          r | m==0 = d
+            | otherwise = error "matrixFromVector"
+
+matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, dat = v}
+    where (d,m) = dim v `divMod` c
+          r | m==0 = d
+            | otherwise = error "matrixFromVector"
+
+
+
+
+
+
 createMatrix order r c = do
     p <- createVector (r*c)
     return (matrixFromVector order c p)
@@ -178,10 +226,10 @@ reshape c v = matrixFromVector RowMajor c v
 
 singleton x = reshape 1 (fromList [x])
 
-liftMatrix :: (Field a, Field b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
+--liftMatrix :: (Field a, Field b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
 liftMatrix f m = reshape (cols m) (f (cdat m))
 
-liftMatrix2 :: (Field t) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+--liftMatrix2 :: (Field t) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
 liftMatrix2 f m1 m2 | compat m1 m2 = reshape (cols m1) (f (cdat m1) (cdat m2))
                     | otherwise    = error "nonconformant matrices in liftMatrix2"
 ------------------------------------------------------------------
@@ -203,6 +251,7 @@ multiplyG a b = matrixFromVector RowMajor (cols b) $ fromList $ concat $ multipl
 
 ------------------------------------------------------------------
 
+{-
 gmatC m f | fortran m =
                 if (isTrans m)
                     then f 0 (rows m) (cols m) (ptr (dat m))
@@ -211,7 +260,11 @@ gmatC m f | fortran m =
                 if isTrans m
                     then f 1 (cols m) (rows m) (ptr (dat m))
                     else f 0 (rows m) (cols m) (ptr (dat m))
+-}
 
+gmatC MF {rows = r, cols = c, dat = d} f = f 1 c r (ptr d)
+gmatC MC {rows = r, cols = c, dat = d} f = f 0 r c (ptr d)
+{-# INLINE gmatC #-}
 
 multiplyAux fun order a b = unsafePerformIO $ do
     when (cols a /= rows b) $ error $ "inconsistent dimensions in contraction "++
@@ -219,6 +272,7 @@ multiplyAux fun order a b = unsafePerformIO $ do
     r <- createMatrix order (rows a) (cols b)
     fun // gmatC a // gmatC b // mat dat r // check "multiplyAux" [dat a, dat b]
     return r
+{-# INLINE multiplyAux #-}
 
 foreign import ccall safe "aux.h multiplyR"
     cmultiplyR :: Int -> Int -> Int -> Ptr Double
@@ -234,13 +288,15 @@ foreign import ccall safe "aux.h multiplyC"
 
 multiply :: (Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a
 multiply RowMajor a b    = multiplyD RowMajor a b
-multiply ColumnMajor a b = m {rows = cols m, cols = rows m, order = ColumnMajor}
-    where m = multiplyD RowMajor (trans b) (trans a)
+multiply ColumnMajor a b = MF {rows = c, cols = r, dat = d}
+    where MC {rows = r, cols = c, dat = d } = multiplyD RowMajor (trans b) (trans a)
 
 
-multiplyR = multiplyAux cmultiplyR
+multiplyR = multiplyAux cmultiplyR'
 multiplyC = multiplyAux cmultiplyC
 
+cmultiplyR' p1 p2 p3 p4 q1 q2 q3 q4 r1 r2 r3 = {-# SCC "mulR" #-} cmultiplyR p1 p2 p3 p4 q1 q2 q3 q4 r1 r2 r3
+
 ----------------------------------------------------------------------
 
 -- | extraction of a submatrix of a real matrix
@@ -249,7 +305,7 @@ subMatrixR :: (Int,Int) -- ^ (r0,c0) starting position
            -> Matrix Double -> Matrix Double
 subMatrixR (r0,c0) (rt,ct) x = unsafePerformIO $ do
     r <- createMatrix RowMajor rt ct
-    c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat cdat r // check "subMatrixR" [dat r]
+    c_submatrixR r0 (r0+rt-1) c0 (c0+ct-1) // mat cdat x // mat dat r // check "subMatrixR" [dat r]
     return r
 foreign import ccall "aux.h submatrixR" c_submatrixR :: Int -> Int -> Int -> Int -> TMM
 
@@ -278,8 +334,8 @@ subMatrixG (r0,c0) (rt,ct) x = reshape ct $ fromList $ concat $ map (subList c0
 
 diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
     m <- createMatrix RowMajor n n
-    fun // vec v // mat dat m // check msg [dat m]
-    return m {tdat = dat m}
+    fun // vec v // mat cdat m // check msg [dat m]
+    return m -- {tdat = dat m}
 
 -- | diagonal matrix from a real vector
 diagR :: Vector Double -> Matrix Double
@@ -305,13 +361,13 @@ diagG v = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..
                     | otherwise = 0
 
 -- | creates a Matrix from a list of vectors
-fromRows :: Field t => [Vector t] -> Matrix t
+--fromRows :: Field t => [Vector t] -> Matrix t
 fromRows vs = case common dim vs of
     Nothing -> error "fromRows applied to [] or to vectors with different sizes"
     Just c  -> reshape c (join vs)
 
 -- | extracts the rows of a matrix as a list of vectors
-toRows :: Storable t => Matrix t -> [Vector t]
+--toRows :: Storable t => Matrix t -> [Vector t]
 toRows m = toRows' 0 where
     v = cdat m
     r = rows m
@@ -324,16 +380,25 @@ fromColumns :: Field t => [Vector t] -> Matrix t
 fromColumns m = trans . fromRows $ m
 
 -- | Creates a list of vectors from the columns of a matrix
-toColumns :: Storable t => Matrix t -> [Vector t]
+toColumns :: Field t => Matrix t -> [Vector t]
 toColumns m = toRows . trans $ m
 
 
 -- | Reads a matrix position.
 (@@>) :: Storable t => Matrix t -> (Int,Int) -> t
 infixl 9 @@>
-m@M {rows = r, cols = c} @@> (i,j)
+--m@M {rows = r, cols = c} @@> (i,j)
+--    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
+--    | otherwise   = cdat m `at` (i*c+j)
+
+MC {rows = r, cols = c, dat = v} @@> (i,j)
+    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
+    | otherwise                  = v `at` (i*c+j)
+
+MF {rows = r, cols = c, dat = v} @@> (i,j)
     | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
-    | otherwise   = cdat m `at` (i*c+j)
+    | otherwise                  = v `at` (j*r+i)
+
 
 ------------------------------------------------------------------
 
diff --git a/lib/Data/Packed/Internal/Tensor.hs b/lib/Data/Packed/Internal/Tensor.hs
index 6876685..dea1636 100644
--- a/lib/Data/Packed/Internal/Tensor.hs
+++ b/lib/Data/Packed/Internal/Tensor.hs
@@ -92,7 +92,7 @@ tensor dssig vec = T d v `withIdx` seqind where
 tensorFromVector :: IdxType -> Vector t -> Tensor t
 tensorFromVector tp v = T {dims = [IdxDesc (dim v) tp "1"], ten = v}
 
-tensorFromMatrix :: IdxType -> IdxType -> Matrix t -> Tensor t
+tensorFromMatrix :: Field t => IdxType -> IdxType -> Matrix t -> Tensor t
 tensorFromMatrix tpr tpc m = T {dims = [IdxDesc (rows m) tpr "1",IdxDesc (cols m) tpc "2"]
                                , ten = cdat m}
 
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 2e8cb3d..45aaaba 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -77,7 +77,7 @@ diagRect s r c
     | r > c     = joinVert  [diag s , zeros (r-c,c)]
     where zeros (r,c) = reshape c $ constant 0 (r*c)
 
-takeDiag :: (Storable t) => Matrix t -> Vector t
+takeDiag :: (Field t) => Matrix t -> Vector t
 takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
 
 ident :: (Num t, Field t) => Int -> Matrix t
@@ -119,7 +119,7 @@ dropColumns n mat = subMatrix (0,n) (rows mat, cols mat - n) mat
 @\> flatten ('ident' 3)
 9 # [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
 -}
-flatten :: Matrix t -> Vector t
+flatten :: Field t => Matrix t -> Vector t
 flatten = cdat
 
 -- | Creates a 'Matrix' from a list of lists (considered as rows).
diff --git a/lib/GSL.hs b/lib/GSL.hs
index 23865b0..e6ff6df 100644
--- a/lib/GSL.hs
+++ b/lib/GSL.hs
@@ -21,7 +21,7 @@ module LinearAlgebra.Algorithms,
 module LAPACK,
 module GSL.Integration,
 module GSL.Differentiation,
-module GSL.Special,
+--module GSL.Special,
 module GSL.Fourier,
 module GSL.Polynomials,
 module GSL.Minimization,
diff --git a/lib/GSL/Compat.hs b/lib/GSL/Compat.hs
index 809a1f5..1d6f7b9 100644
--- a/lib/GSL/Compat.hs
+++ b/lib/GSL/Compat.hs
@@ -38,7 +38,7 @@ adaptScalar f1 f2 f3 x y
     | dim y == 1 = f3 x (y@>0)
     | otherwise = f2 x y
 
-liftMatrix2' :: (Field t) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
 liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (cdat m1) (cdat m2))
                      | otherwise    = error "nonconformant matrices in liftMatrix2'"
 
@@ -63,6 +63,7 @@ instance (Eq a, Field a) => Eq (Matrix a) where
 
 instance (Field a, Linear Vector a) => Num (Matrix a) where
     (+) = liftMatrix2' (+)
+    (-) = liftMatrix2' (-)
     negate = liftMatrix negate
     (*) = liftMatrix2' (*)
     signum = liftMatrix signum
diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs
index 26c5e2a..15710df 100644
--- a/lib/GSL/Matrix.hs
+++ b/lib/GSL/Matrix.hs
@@ -46,13 +46,14 @@ import Foreign.C.String
 
 -}
 eigSg :: Matrix Double -> (Vector Double, Matrix Double)
-eigSg (m@M {rows = r})
+eigSg m
     | r == 1 = (fromList [cdat m `at` 0], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix RowMajor r r
         c_eigS // mat cdat m // vec l // mat dat v // check "eigSg" [cdat m]
         return (l,v)
+  where r = rows m
 foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
 
 ------------------------------------------------------------------
@@ -76,13 +77,14 @@ foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
 
 -}
 eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double))
-eigHg (m@M {rows = r})
+eigHg m
     | r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix RowMajor r r
         c_eigH // mat cdat m // vec l // mat dat v // check "eigHg" [cdat m]
         return (l,v)
+  where r = rows m
 foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
 
 
@@ -108,16 +110,18 @@ foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
 
 -}
 svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdg x@M {rows = r, cols = c} = if r>=c
+svdg x = if rows x >= cols x
     then svd' x
     else (v, s, u) where (u,s,v) = svd' (trans x)
 
-svd' x@M {rows = r, cols = c} = unsafePerformIO $ do
+svd' x = unsafePerformIO $ do
     u <- createMatrix RowMajor r c
     s <- createVector c
     v <- createMatrix RowMajor c c
     c_svd // mat cdat x // mat dat u // vec s // mat dat v // check "svdg" [cdat x]
     return (u,s,v)
+  where r = rows x
+        c = cols x
 foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
 
 {- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/.
@@ -138,11 +142,13 @@ foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
 
 -}
 qr :: Matrix Double -> (Matrix Double, Matrix Double)
-qr x@M {rows = r, cols = c} = unsafePerformIO $ do
+qr x = unsafePerformIO $ do
     q <- createMatrix RowMajor r r
     rot <- createMatrix RowMajor r c
     c_qr // mat cdat x // mat dat q // mat dat rot // check "qr" [cdat x]
     return (q,rot)
+  where r = rows x
+        c = cols x
 foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
 
 {- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
@@ -159,11 +165,11 @@ foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
 
 -}
 chol :: Matrix Double -> Matrix Double
---chol x@(M r _ p) = createM [p] "chol" r r $ m c_chol x
-chol x@M {rows = r} = unsafePerformIO $ do
+chol x = unsafePerformIO $ do
     res <- createMatrix RowMajor r r
     c_chol // mat cdat x // mat dat res // check "chol" [cdat x]
     return res
+  where r = rows x
 foreign import ccall "gsl-aux.h chol" c_chol :: TMM
 
 --------------------------------------------------------
@@ -171,43 +177,53 @@ foreign import ccall "gsl-aux.h chol" c_chol :: TMM
 {- -| efficient multiplication by the inverse of a matrix (for real matrices)
 -}
 luSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-luSolveR  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
+luSolveR a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix RowMajor r c
         c_luSolveR // mat cdat a // mat cdat b // mat dat s // check "luSolveR" [cdat a, cdat b]
         return s
     | otherwise = error "luSolveR of nonsquare matrix"
-
+  where n1 = rows a
+        n2 = cols a
+        r  = rows b
+        c  = cols b
 foreign import ccall "gsl-aux.h luSolveR" c_luSolveR ::  TMMM
 
 {- -| efficient multiplication by the inverse of a matrix (for complex matrices). 
 -}
 luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-luSolveC  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
+luSolveC a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix RowMajor r c
         c_luSolveC // mat cdat a // mat cdat b // mat dat s // check "luSolveC" [cdat a, cdat b]
         return s
     | otherwise = error "luSolveC of nonsquare matrix"
-
+  where n1 = rows a
+        n2 = cols a
+        r  = rows b
+        c  = cols b
 foreign import ccall "gsl-aux.h luSolveC" c_luSolveC ::  TCMCMCM
 
 {- | lu decomposition of real matrix (packed as a vector including l, u, the permutation and sign)
 -}
 luRaux  :: Matrix Double -> Vector Double
-luRaux x@M {rows = r, cols = c} = unsafePerformIO $ do
+luRaux x = unsafePerformIO $ do
     res <- createVector (r*r+r+1)
     c_luRaux // mat cdat x // vec res // check "luRaux" [cdat x]
     return res
+  where r = rows x
+        c = cols x
 foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV
 
 {- | lu decomposition of complex matrix (packed as a vector including l, u, the permutation and sign)
 -}
 luCaux  :: Matrix (Complex Double) -> Vector (Complex Double)
-luCaux x@M {rows = r, cols = c} = unsafePerformIO $ do
+luCaux x = unsafePerformIO $ do
     res <- createVector (r*r+r+1)
     c_luCaux // mat cdat x // vec res // check "luCaux" [cdat x]
     return res
+  where r = rows x
+        c = cols x
 foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV
 
 {- | The LU decomposition of a square matrix. Is based on /gsl_linalg_LU_decomp/ and  /gsl_linalg_complex_LU_decomp/ as described in <http://www.gnu.org/software/gsl/manual/gsl-ref_13.html#SEC223>.
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index b0008b1..ba72681 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -37,16 +37,19 @@ foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM
 --
 -- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
 svdR :: Matrix Double -> (Matrix Double, Matrix Double, Matrix Double)
-svdR x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdR' x
+svdR x = (u, diagRect s r c, v) where (u,s,v) = svdR' x
+                                      r = rows x
+                                      c = cols x
 
 svdR' :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
+svdR' x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
     dgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdR" [fdat x]
     return (u,s,trans v)
-
+  where r = rows x
+        c = cols x
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM
 
@@ -54,15 +57,19 @@ foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM
 --
 -- @(u,s,v)=svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@.
 svdRdd :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
-svdRdd x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdRdd' x
+svdRdd x = (u, diagRect s r c, v) where (u,s,v) = svdRdd' x
+                                        r = rows x
+                                        c = cols x
 
 svdRdd' :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdRdd' x@M {rows = r, cols = c} = unsafePerformIO $ do
+svdRdd' x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
     dgesdd // mat fdat x // mat dat u // vec s // mat dat v // check "svdRdd" [fdat x]
     return (u,s,trans v)
+  where r = rows x
+        c = cols x
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
@@ -71,15 +78,20 @@ foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
 --
 -- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@.
 svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double))
-svdC x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdC' x
+svdC x = (u, diagRect s r c, v) where (u,s,v) = svdC' x
+                                      r = rows x
+                                      c = cols x
 
 svdC' :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
+svdC' x = unsafePerformIO $ do
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
     zgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdC" [fdat x]
     return (u,s,trans v)
+  where r = rows x
+        c = cols x
+
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
@@ -91,7 +103,7 @@ foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
 -- The eigenvectors are the columns of v.
 -- The eigenvalues are not sorted.
 eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
-eigC (m@M {rows = r}) 
+eigC m
     | r == 1 = (fromList [cdat m `at` 0], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
@@ -99,6 +111,7 @@ eigC (m@M {rows = r})
         dummy <- createMatrix ColumnMajor 1 1
         zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m]
         return (l,v)
+  where r = rows m
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
@@ -110,14 +123,15 @@ foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
 -- The eigenvectors are the columns of v.
 -- The eigenvalues are not sorted.
 eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
-eigR (m@M {rows = r}) = (s', v'')
+eigR m = (s', v'')
     where (s,v) = eigRaux m
           s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
           v' = toRows $ trans v
           v'' = fromColumns $ fixeig (toList s') v'
+          r = rows m
 
 eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
-eigRaux (m@M {rows = r})
+eigRaux m
     | r == 1 = (fromList [(cdat m `at` 0):+0], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
@@ -125,6 +139,7 @@ eigRaux (m@M {rows = r})
         dummy <- createMatrix ColumnMajor 1 1
         dgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigR" [fdat m]
         return (l,v)
+  where r = rows m
 
 fixeig  []  _ =  []
 fixeig [r] [v] = [comp v]
@@ -148,13 +163,14 @@ eigS m = (s', fliprl v)
     where (s,v) = eigS' m
           s' = fromList . reverse . toList $  s
 
-eigS' (m@M {rows = r})
+eigS' m
     | r == 1 = (fromList [cdat m `at` 0], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
         dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m]
         return (l,v)
+  where r = rows m
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
@@ -170,37 +186,46 @@ eigH m = (s', fliprl v)
     where (s,v) = eigH' m
           s' = fromList . reverse . toList $  s
 
-eigH' (m@M {rows = r})
+eigH' m
     | r == 1 = (fromList [realPart (cdat m `at` 0)], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
         zheev // mat fdat m // vec l // mat dat v // check "eigH" [fdat m]
         return (l,v)
+  where r = rows m
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
 
 -- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition.
 linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-linearSolveR  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
+linearSolveR  a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix ColumnMajor r c
         dgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveR" [fdat a, fdat b]
         return s
     | otherwise = error "linearSolveR of nonsquare matrix"
+  where n1 = rows a
+        n2 = cols a
+        r  = rows b
+        c  = cols b
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
 
 -- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition.
 linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-linearSolveC  a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
+linearSolveC  a b
     | n1==n2 && n1==r = unsafePerformIO $ do
         s <- createMatrix ColumnMajor r c
         zgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveC" [fdat a, fdat b]
         return s
     | otherwise = error "linearSolveC of nonsquare matrix"
+  where n1 = rows a
+        n2 = cols a
+        r  = rows b
+        c  = cols b
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM
@@ -209,10 +234,13 @@ foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM
 linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
 linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSR_l a b
 
-linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
+linearSolveLSR_l a b = unsafePerformIO $ do
     r <- createMatrix ColumnMajor (max m n) nrhs
     dgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSR" [fdat a, fdat b]
     return r
+  where m = rows a
+        n = cols a
+        nrhs = cols b
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM
@@ -221,10 +249,13 @@ foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM
 linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
 linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSC_l a b
 
-linearSolveLSC_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
+linearSolveLSC_l a b = unsafePerformIO $ do
     r <- createMatrix ColumnMajor (max m n) nrhs
     zgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSC" [fdat a, fdat b]
     return r
+  where m = rows a
+        n = cols a
+        nrhs = cols b
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM
@@ -237,10 +268,13 @@ linearSolveSVDR :: Maybe Double   -- ^ rcond
 linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDR_l rcond a b
 linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) a b
 
-linearSolveSVDR_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
+linearSolveSVDR_l rcond a b = unsafePerformIO $ do
     r <- createMatrix ColumnMajor (max m n) nrhs
     dgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDR" [fdat a, fdat b]
     return r
+  where m = rows a
+        n = cols a
+        nrhs = cols b
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
@@ -253,8 +287,11 @@ linearSolveSVDC :: Maybe Double            -- ^ rcond
 linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDC_l rcond a b
 linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) a b
 
-linearSolveSVDC_l rcond  a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
+linearSolveSVDC_l rcond  a b = unsafePerformIO $ do
     r <- createMatrix ColumnMajor (max m n) nrhs
     zgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDC" [fdat a, fdat b]
     return r
+  where m = rows a
+        n = cols a
+        nrhs = cols b