5 files changed, 155 insertions, 187 deletions

diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
index 5a5c19e..e803c53 100644
--- a/lib/Numeric/GSL/Matrix.hs
+++ b/lib/Numeric/GSL/Matrix.hs
@@ -44,12 +44,14 @@ import Complex
 
 -}
 eigSg :: Matrix Double -> (Vector Double, Matrix Double)
-eigSg m
+eigSg = eigSg . cmat
+
+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]
+        c_eigS // matc m // vec l // matc v // check "eigSg" [cdat m]
         return (l,v)
   where r = rows m
 foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
@@ -75,12 +77,14 @@ foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
 
 -}
 eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double))
-eigHg m
+eigHg = eigHg' . cmat
+
+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]
+        c_eigH // matc m // vec l // matc v // check "eigHg" [cdat m]
         return (l,v)
   where r = rows m
 foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
@@ -109,14 +113,14 @@ foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
 -}
 svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
 svdg x = if rows x >= cols x
-    then svd' x
-    else (v, s, u) where (u,s,v) = svd' (trans x)
+    then svd' (cmat x)
+    else (v, s, u) where (u,s,v) = svd' (cmat (trans x))
 
 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]
+    c_svd // matc x // matc u // vec s // matc v // check "svdg" [cdat x]
     return (u,s,v)
   where r = rows x
         c = cols x
@@ -140,30 +144,36 @@ foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
 
 -}
 qr :: Matrix Double -> (Matrix Double, Matrix Double)
-qr x = unsafePerformIO $ do
+qr = qr' . cmat
+
+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]
+    c_qr // matc x // matc q // matc 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
 
 qrPacked :: Matrix Double -> (Matrix Double, Vector Double)
-qrPacked x = unsafePerformIO $ do
+qrPacked = qrPacked' . cmat
+
+qrPacked' x = unsafePerformIO $ do
     qr <- createMatrix RowMajor r c
     tau <- createVector (min r c)
-    c_qrPacked // mat cdat x // mat dat qr // vec tau // check "qrUnpacked" [cdat x]
+    c_qrPacked // matc x // matc qr // vec tau // check "qrUnpacked" [cdat x]
     return (qr,tau)
   where r = rows x
         c = cols x
 foreign import ccall "gsl-aux.h QRpacked" c_qrPacked :: TMMV
 
 unpackQR :: (Matrix Double, Vector Double) -> (Matrix Double, Matrix Double)
-unpackQR (qr,tau) = unsafePerformIO $ do
+unpackQR (qr,tau) = unpackQR' (cmat qr, tau)
+
+unpackQR' (qr,tau) = unsafePerformIO $ do
     q <- createMatrix RowMajor r r
     rot <- createMatrix RowMajor r c
-    c_qrUnpack // mat cdat qr // vec tau // mat dat q // mat dat rot // check "qrUnpack" [cdat qr,tau]
+    c_qrUnpack // matc qr // vec tau // matc q // matc rot // check "qrUnpack" [cdat qr,tau]
     return (q,rot)
   where r = rows qr
         c = cols qr
@@ -183,17 +193,21 @@ type TMVMM = Int -> Int -> PD -> Int -> PD -> TMM
 
 -}
 cholR :: Matrix Double -> Matrix Double
-cholR x = unsafePerformIO $ do
+cholR = cholR' . cmat
+
+cholR' x = unsafePerformIO $ do
     res <- createMatrix RowMajor r r
-    c_cholR // mat cdat x // mat dat res // check "cholR" [cdat x]
+    c_cholR // matc x // matc res // check "cholR" [cdat x]
     return res
   where r = rows x
 foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
 
 cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
-cholC x = unsafePerformIO $ do
+cholC = cholC' . cmat
+
+cholC' x = unsafePerformIO $ do
     res <- createMatrix RowMajor r r
-    c_cholC // mat cdat x // mat dat res // check "cholC" [cdat x]
+    c_cholC // matc x // matc res // check "cholC" [cdat x]
     return res
   where r = rows x
 foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
@@ -204,10 +218,12 @@ foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
 {- -| efficient multiplication by the inverse of a matrix (for real matrices)
 -}
 luSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-luSolveR a b
+luSolveR a b = luSolveR' (cmat a) (cmat b)
+
+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]
+        c_luSolveR // matc a // matc b // matc s // check "luSolveR" [cdat a, cdat b]
         return s
     | otherwise = error "luSolveR of nonsquare matrix"
   where n1 = rows a
@@ -219,10 +235,12 @@ 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 b
+luSolveC a b = luSolveC' (cmat a) (cmat b)
+
+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]
+        c_luSolveC // matc a // matc b // matc s // check "luSolveC" [cdat a, cdat b]
         return s
     | otherwise = error "luSolveC of nonsquare matrix"
   where n1 = rows a
@@ -234,9 +252,11 @@ 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 = unsafePerformIO $ do
+luRaux = luRaux' . cmat
+
+luRaux' x = unsafePerformIO $ do
     res <- createVector (r*r+r+1)
-    c_luRaux // mat cdat x // vec res // check "luRaux" [cdat x]
+    c_luRaux // matc x // vec res // check "luRaux" [cdat x]
     return res
   where r = rows x
         c = cols x
@@ -245,9 +265,11 @@ 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 = unsafePerformIO $ do
+luCaux = luCaux' . cmat
+
+luCaux' x = unsafePerformIO $ do
     res <- createVector (r*r+r+1)
-    c_luCaux // mat cdat x // vec res // check "luCaux" [cdat x]
+    c_luCaux // matc x // vec res // check "luCaux" [cdat x]
     return res
   where r = rows x
         c = cols x
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index e93f8cb..f523849 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -186,12 +186,12 @@ foreign import ccall "wrapper"
 
 aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO())
 aux_vTov f n p r = g where
-    V {fptr = pr, ptr = t} = f x
+    v@V {fptr = pr} = f x
     x = createV n copy "aux_vTov" []
     copy n q = do
         copyArray q p n
         return 0
-    g = withForeignPtr pr $ \_ -> copyArray r t n
+    g = withForeignPtr pr $ \_ -> copyArray r (ptr v) n
 
 --------------------------------------------------------------------
 
@@ -202,10 +202,10 @@ createV n fun msg ptrs = unsafePerformIO $ do
 
 createM r c fun msg ptrs = unsafePerformIO $ do
     r <- createMatrix RowMajor r c
-    fun // mat cdat r // check msg ptrs
+    fun // matc r // check msg ptrs
     return r
 
 createMIO r c fun msg ptrs = do
     r <- createMatrix RowMajor r c
-    fun // mat cdat r // check msg ptrs
+    fun // matc r // check msg ptrs
     return r
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 628d4f8..315be17 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -1,4 +1,4 @@
-{-# OPTIONS_GHC -fglasgow-exts #-}
+{-# OPTIONS_GHC #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.LinearAlgebra.LAPACK
@@ -36,54 +36,52 @@ import Foreign
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM
+foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
+foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM
 
 -- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
 --
 -- @(u,s,v)=full svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
 svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-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
+svdR = svdAux dgesvd "svdR" . fmat
 
 -- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
 --
 -- @(u,s,v)=full svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@.
 svdRdd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-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
+svdRdd = svdAux dgesdd "svdRdd" . fmat
 
 -- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
 --
 -- @(u,s,v)=full svdC m@ so that @m=u \<\> comp s \<\> 'trans' v@.
 svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-svdC x = unsafePerformIO $ do
+svdC = svdAux zgesvd "svdC" . fmat
+
+svdAux f st 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]
+    f // matf x // matf u // vec s // matf v // check st [fdat x]
     return (u,s,trans v)
   where r = rows x
         c = cols x
 
-
 -----------------------------------------------------------------------------
+eigAux f st m
+    | r == 1 = (fromList [flatten m `at` 0], singleton 1)
+    | otherwise = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix ColumnMajor r r
+        dummy <- createMatrix ColumnMajor 1 1
+        f // matf m // matf dummy // vec l // matf v // check st [fdat m]
+        return (l,v)
+  where r = rows m
+
+
 foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
 
 -- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
 --
@@ -92,18 +90,9 @@ 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
-    | r == 1 = (fromList [cdat m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r 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
+eigC = eigAux zgeev "eigC" . fmat
 
 -----------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
 
 -- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
 --
@@ -113,7 +102,7 @@ foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
 -- The eigenvalues are not sorted.
 eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
 eigR m = (s', v'')
-    where (s,v) = eigRaux m
+    where (s,v) = eigRaux (fmat m)
           s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
           v' = toRows $ trans v
           v'' = fromColumns $ fixeig (toList s') v'
@@ -121,12 +110,12 @@ eigR m = (s', v'')
 
 eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
 eigRaux m
-    | r == 1 = (fromList [(cdat m `at` 0):+0], singleton 1)
+    | r == 1 = (fromList [(flatten m `at` 0):+0], singleton 1)
     | otherwise = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
         dummy <- createMatrix ColumnMajor 1 1
-        dgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigR" [fdat m]
+        dgeev // matf m // matf dummy // vec l // matf v // check "eigR" [fdat m]
         return (l,v)
   where r = rows m
 
@@ -138,7 +127,6 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
   where scale = vectorMapValR Scale
 
 -----------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
 
 -- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
 --
@@ -148,20 +136,19 @@ foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
 -- The eigenvalues are sorted in descending order (use eigS' for ascending order).
 eigS :: Matrix Double -> (Vector Double, Matrix Double)
 eigS m = (s', fliprl v)
-    where (s,v) = eigS' m
+    where (s,v) = eigS' (fmat m)
           s' = fromList . reverse . toList $  s
 
 eigS' m
-    | r == 1 = (fromList [cdat m `at` 0], singleton 1)
+    | r == 1 = (fromList [flatten 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]
+        dsyev // matf m // vec l // matf v // check "eigS" [fdat m]
         return (l,v)
   where r = rows m
 
 -----------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
 
 -- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix:
 --
@@ -171,165 +158,120 @@ foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
 -- The eigenvalues are sorted in descending order (use eigH' for ascending order).
 eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
 eigH m = (s', fliprl v)
-    where (s,v) = eigH' m
+    where (s,v) = eigH' (fmat m)
           s' = fromList . reverse . toList $  s
 
 eigH' m
-    | r == 1 = (fromList [realPart (cdat m `at` 0)], singleton 1)
+    | r == 1 = (fromList [realPart (flatten 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]
+        zheev // matf m // vec l // matf v // check "eigH" [fdat m]
         return (l,v)
   where r = rows m
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
+foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
 
--- | 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 b
+linearSolveSQAux f st 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]
+        f // matf a // matf b // matf s // check st [fdat a, fdat b]
         return s
-    | otherwise = error "linearSolveR of nonsquare matrix"
+    | otherwise = error $ st ++ " 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 /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 b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b)
 
 -- | 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 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
+linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b)
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM
+foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM
+foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM
+foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
 
--- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'.
-linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
-linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSR_l a b
-
-linearSolveLSR_l a b = unsafePerformIO $ do
+linearSolveAux f st 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]
+    f // matf a // matf b // matf r // check st [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
+-- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'.
+linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
+linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $
+                     linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b)
 
 -- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'.
 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 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
+linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $
+                     linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b)
 
 -- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
 linearSolveSVDR :: Maybe Double   -- ^ rcond
                 -> Matrix Double  -- ^ coefficient matrix
                 -> Matrix Double  -- ^ right hand sides (as columns)
                 -> Matrix Double  -- ^ solution vectors (as columns)
-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 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
+linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
+                                   linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b)
+linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b)
 
 -- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
 linearSolveSVDC :: Maybe Double            -- ^ rcond
                 -> Matrix (Complex Double) -- ^ coefficient matrix
                 -> Matrix (Complex Double) -- ^ right hand sides (as columns)
                 -> Matrix (Complex Double) -- ^ solution vectors (as columns)
-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 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
+linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
+                                   linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b)
+linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM
+foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
 
 -- | Wrapper for LAPACK's /zpotrf/, which computes the Cholesky factorization of a
 -- complex Hermitian positive definite matrix.
 cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
-cholH a = unsafePerformIO $ do
-    r <- createMatrix ColumnMajor n n
-    zpotrf // mat fdat a // mat dat r // check "cholH" [fdat a]
-    return r
-  where n = rows a
-
------------------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
+cholH = cholAux zpotrf "cholH" . fmat
 
 -- | Wrapper for LAPACK's /dpotrf/, which computes the Cholesky factorization of a
 -- real symmetric positive definite matrix.
 cholS :: Matrix Double -> Matrix Double
-cholS a = unsafePerformIO $ do
+cholS = cholAux dpotrf "cholS" . fmat
+
+cholAux f st a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor n n
-    dpotrf // mat fdat a // mat dat r // check "cholS" [fdat a]
+    f // matf a // matf r // check st [fdat a]
     return r
   where n = rows a
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM
+foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM
 
 -- | Wrapper for LAPACK's /dgeqr2/, which computes a QR factorization of a real matrix.
 qrR :: Matrix Double -> (Matrix Double, Vector Double)
-qrR a = unsafePerformIO $ do
-    r <- createMatrix ColumnMajor m n
-    tau <- createVector mn
-    dgeqr2 // mat fdat a // vec tau // mat dat r // check "qrR" [fdat a]
-    return (r,tau)
-  where m = rows a
-        n = cols a
-        mn = min m n
-
------------------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM
+qrR = qrAux dgeqr2 "qrR" . fmat
 
 -- | Wrapper for LAPACK's /zgeqr2/, which computes a QR factorization of a complex matrix.
 qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
-qrC a = unsafePerformIO $ do
+qrC = qrAux zgeqr2 "qrC" . fmat
+
+qrAux f st a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor m n
     tau <- createVector mn
-    zgeqr2 // mat fdat a // vec tau // mat dat r // check "qrC" [fdat a]
+    withForeignPtr (fptr $ fdat $ a) $ \p ->
+        f m n p // vec tau // matf r // check st [fdat a]
     return (r,tau)
   where m = rows a
         n = cols a
@@ -337,52 +279,42 @@ qrC a = unsafePerformIO $ do
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM
+foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM
 
 -- | Wrapper for LAPACK's /dgehrd/, which computes a Hessenberg factorization of a square real matrix.
 hessR :: Matrix Double -> (Matrix Double, Vector Double)
-hessR a = unsafePerformIO $ do
-    r <- createMatrix ColumnMajor m n
-    tau <- createVector (mn-1)
-    dgehrd // mat fdat a // vec tau // mat dat r // check "hessR" [fdat a]
-    return (r,tau)
-  where m = rows a
-        n = cols a
-        mn = min m n
-
------------------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM
+hessR = hessAux dgehrd "hessR" . fmat
 
 -- | Wrapper for LAPACK's /zgehrd/, which computes a Hessenberg factorization of a square complex matrix.
 hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
-hessC a = unsafePerformIO $ do
+hessC = hessAux zgehrd "hessC" . fmat
+
+hessAux f st a = unsafePerformIO $ do
     r <- createMatrix ColumnMajor m n
     tau <- createVector (mn-1)
-    zgehrd // mat fdat a // vec tau // mat dat r // check "hessC" [fdat a]
+    f // matf a // vec tau // matf r // check st [fdat a]
     return (r,tau)
   where m = rows a
         n = cols a
         mn = min m n
 
 -----------------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
+foreign import ccall safe "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
+foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM
 
 -- | Wrapper for LAPACK's /dgees/, which computes a Schur factorization of a square real matrix.
 schurR :: Matrix Double -> (Matrix Double, Matrix Double)
-schurR a = unsafePerformIO $ do
-    u <- createMatrix ColumnMajor n n
-    s <- createMatrix ColumnMajor n n
-    dgees // mat fdat a // mat dat u // mat dat s // check "schurR" [fdat a]
-    return (u,s)
-  where n = rows a
-
------------------------------------------------------------------------------------
-foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM
+schurR = schurAux dgees "schurR" . fmat
 
 -- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix.
 schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
-schurC a = unsafePerformIO $ do
+schurC = schurAux zgees "schurC" . fmat
+
+schurAux f st a = unsafePerformIO $ do
     u <- createMatrix ColumnMajor n n
     s <- createMatrix ColumnMajor n n
-    zgees // mat fdat a // mat dat u // mat dat s // check "schurC" [fdat a]
+    f // matf a // matf u // matf s // check st [fdat a]
     return (u,s)
   where n = rows a
+
+-----------------------------------------------------------------------------------
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 058232c..8392feb 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -709,6 +709,11 @@ int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
     integer n = ac;
     REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
     DEBUGMSG("schur_l_R");
+    int k;
+    //printf("---------------------------\n");
+    //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
+    //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
+    //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
     memcpy(sp,ap,n*n*sizeof(double));
     integer lwork = 6*n; // fixme
     double *WORK = (double*)malloc(lwork*sizeof(double));
@@ -719,6 +724,12 @@ int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
     integer res;
     integer sdim;
     dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);
+    //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
+    //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
+    //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
+    if(res>0) {
+        return NOCONVER;
+    }
     CHECK(res,res);
     free(WR);
     free(WI);
@@ -747,6 +758,9 @@ int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
     zgees_ ("V","N",NULL,&n,(doublecomplex*)sp,&n,&sdim,W,
                             (doublecomplex*)up,&n,
                             WORK,&lwork,RWORK,BWORK,&res);
+    if(res>0) {
+        return NOCONVER;
+    }
     CHECK(res,res);
     free(W);
     free(BWORK);
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 13d69ab..3403591 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -20,7 +20,7 @@ module Numeric.LinearAlgebra.Linear (
 ) where
 
 
-import Data.Packed.Internal
+import Data.Packed.Internal(multiply,at,partit)
 import Data.Packed
 import Numeric.GSL.Vector
 import Complex
@@ -68,7 +68,7 @@ instance Linear Matrix Double where
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
-    equal a b = cols a == cols b && cdat a `equal` cdat b
+    equal a b = cols a == cols b && flatten a `equal` flatten b
 
 
 instance Linear Matrix (Complex Double) where
@@ -79,13 +79,13 @@ instance Linear Matrix (Complex Double) where
     sub = liftMatrix2 sub
     mul = liftMatrix2 mul
     divide = liftMatrix2 divide
-    equal a b = cols a == cols b && cdat a `equal` cdat b
+    equal a b = cols a == cols b && flatten a `equal` flatten b
 
 --------------------------------------------------
 
 -- | euclidean inner product
 dot :: (Element t) => Vector t -> Vector t -> t
-dot u v = dat (multiply r c) `at` 0
+dot u v = multiply r c  @@> (0,0)
     where r = asRow u
           c = asColumn v