Conversion, LAPACK -> base

4 files changed, 656 insertions, 2 deletions

diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs
index 777b6c5..9350acb 100644
--- a/packages/base/src/Data/Packed/Development.hs
+++ b/packages/base/src/Data/Packed/Development.hs
@@ -24,7 +24,7 @@ module Data.Packed.Development (
     unsafeFromForeignPtr,
     unsafeToForeignPtr,
     check, (//),
-    at', atM'
+    at', atM', fi, table
 ) where
 
 import Data.Packed.Internal
diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs
index b5a4318..653a257 100644
--- a/packages/base/src/Data/Packed/Vector.hs
+++ b/packages/base/src/Data/Packed/Vector.hs
@@ -18,7 +18,7 @@
 
 module Data.Packed.Vector (
     Vector,
-    fromList, (|>), toList, buildVector,
+    fromList, (|>), toList, buildVector, constant,
     dim, (@>),
     subVector, takesV, vjoin, join,
     mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith,
@@ -27,6 +27,7 @@ module Data.Packed.Vector (
 ) where
 
 import Data.Packed.Internal.Vector
+import Data.Packed.Internal.Matrix
 import Foreign.Storable
 
 -------------------------------------------------------------------
@@ -94,3 +95,13 @@ unzipVector = unzipVectorWith id
 join ::  Storable t => [Vector t] -> Vector t
 join = vjoin
 
+{- | creates a vector with a given number of equal components:
+
+@> constant 2 7
+7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
+-}
+constant :: Element a => a -> Int -> Vector a
+-- constant x n = runSTVector (newVector x n)
+constant = constantD-- about 2x faster
+
+
diff --git a/packages/base/src/Numeric/Conversion.hs b/packages/base/src/Numeric/Conversion.hs
new file mode 100644
index 0000000..a1f9385
--- /dev/null
+++ b/packages/base/src/Numeric/Conversion.hs
@@ -0,0 +1,91 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Conversion
+-- Copyright   :  (c) Alberto Ruiz 2010
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Conversion routines
+--
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+
+module Numeric.Conversion (
+    Complexable(..), RealElement,
+    module Data.Complex
+) where
+
+import Data.Packed.Internal.Vector
+import Data.Packed.Internal.Matrix
+import Data.Complex
+import Control.Arrow((***))
+
+-------------------------------------------------------------------
+
+-- | Supported single-double precision type pairs
+class (Element s, Element d) => Precision s d | s -> d, d -> s where
+    double2FloatG :: Vector d -> Vector s
+    float2DoubleG :: Vector s -> Vector d
+
+instance Precision Float Double where
+    double2FloatG = double2FloatV
+    float2DoubleG = float2DoubleV
+
+instance Precision (Complex Float) (Complex Double) where
+    double2FloatG = asComplex . double2FloatV . asReal
+    float2DoubleG = asComplex . float2DoubleV . asReal
+
+-- | Supported real types
+class (Element t, Element (Complex t), RealFloat t
+--       , RealOf t ~ t, RealOf (Complex t) ~ t
+       )
+    => RealElement t
+
+instance RealElement Double
+instance RealElement Float
+
+
+-- | Structures that may contain complex numbers
+class Complexable c where
+    toComplex'   :: (RealElement e) => (c e, c e) -> c (Complex e)
+    fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e)
+    comp'        :: (RealElement e) => c e -> c (Complex e)
+    single'      :: Precision a b => c b -> c a
+    double'      :: Precision a b => c a -> c b
+
+
+instance Complexable Vector where
+    toComplex' = toComplexV
+    fromComplex' = fromComplexV
+    comp' v = toComplex' (v,constantD 0 (dim v))
+    single' = double2FloatG
+    double' = float2DoubleG
+
+
+-- | creates a complex vector from vectors with real and imaginary parts
+toComplexV :: (RealElement a) => (Vector a, Vector a) ->  Vector (Complex a)
+toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
+
+-- | the inverse of 'toComplex'
+fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a)
+fromComplexV z = (r,i) where
+    [r,i] = toColumns $ reshape 2 $ asReal z
+
+
+instance Complexable Matrix where
+    toComplex' = uncurry $ liftMatrix2 $ curry toComplex'
+    fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z
+        where c = cols z
+    comp' = liftMatrix comp'
+    single' = liftMatrix single'
+    double' = liftMatrix double'
+
diff --git a/packages/base/src/Numeric/LAPACK.hs b/packages/base/src/Numeric/LAPACK.hs
new file mode 100644
index 0000000..08cd759
--- /dev/null
+++ b/packages/base/src/Numeric/LAPACK.hs
@@ -0,0 +1,552 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.LAPACK
+-- Copyright   :  (c) Alberto Ruiz 2006-14
+-- License     :  BSD3
+-- 
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).
+--
+-----------------------------------------------------------------------------
+
+module Numeric.LAPACK (
+    -- * Matrix product
+    multiplyR, multiplyC, multiplyF, multiplyQ,
+    -- * Linear systems
+    linearSolveR, linearSolveC,
+    lusR, lusC,
+    cholSolveR, cholSolveC,
+    linearSolveLSR, linearSolveLSC,
+    linearSolveSVDR, linearSolveSVDC,
+    -- * SVD
+    svR, svRd, svC, svCd,
+    svdR, svdRd, svdC, svdCd,
+    thinSVDR, thinSVDRd, thinSVDC, thinSVDCd,
+    rightSVR, rightSVC, leftSVR, leftSVC,
+    -- * Eigensystems
+    eigR, eigC, eigS, eigS', eigH, eigH',
+    eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH,
+    -- * LU
+    luR, luC,
+    -- * Cholesky
+    cholS, cholH, mbCholS, mbCholH,
+    -- * QR
+    qrR, qrC, qrgrR, qrgrC,
+    -- * Hessenberg
+    hessR, hessC,
+    -- * Schur
+    schurR, schurC
+) where
+
+import Data.Packed.Development
+import Data.Packed
+import Data.Packed.Internal
+import Numeric.Conversion
+
+import Foreign.Ptr(nullPtr)
+import Foreign.C.Types
+import Control.Monad(when)
+import System.IO.Unsafe(unsafePerformIO)
+
+-----------------------------------------------------------------------------------
+
+foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM
+foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM
+foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM
+foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM
+
+isT Matrix{order = ColumnMajor} = 0
+isT Matrix{order = RowMajor} = 1
+
+tt x@Matrix{order = ColumnMajor} = x
+tt x@Matrix{order = RowMajor} = trans x
+
+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)
+    app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st
+    return s
+
+-- | Matrix product based on BLAS's /dgemm/.
+multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
+multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b
+
+-- | Matrix product based on BLAS's /zgemm/.
+multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+multiplyC a b = multiplyAux zgemmc "zgemmc" a b
+
+-- | Matrix product based on BLAS's /sgemm/.
+multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
+multiplyF a b = multiplyAux sgemmc "sgemmc" a b
+
+-- | Matrix product based on BLAS's /cgemm/.
+multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)
+multiplyQ a b = multiplyAux cgemmc "cgemmc" a b
+
+-----------------------------------------------------------------------------
+foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM
+foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM
+foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM
+foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM
+
+-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
+svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+svdR = svdAux dgesvd "svdR" . fmat
+
+-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
+svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+svdRd = svdAux dgesdd "svdRdd" . fmat
+
+-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
+svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
+svdC = svdAux zgesvd "svdC" . fmat
+
+-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
+svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
+svdCd = svdAux zgesdd "svdCdd" . fmat
+
+svdAux f st x = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor r r
+    s <- createVector (min r c)
+    v <- createMatrix ColumnMajor c c
+    app4 f mat x mat u vec s mat v st
+    return (u,s,trans v)
+  where r = rows x
+        c = cols x
+
+
+-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
+thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat
+
+-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
+thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
+thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat
+
+-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
+thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat
+
+-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
+thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
+thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat
+
+thinSVDAux f st x = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor r q
+    s <- createVector q
+    v <- createMatrix ColumnMajor q c
+    app4 f mat x mat u vec s mat v st
+    return (u,s,trans v)
+  where r = rows x
+        c = cols x
+        q = min r c
+
+
+-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
+svR :: Matrix Double -> Vector Double
+svR = svAux dgesvd "svR" . fmat
+
+-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
+svC :: Matrix (Complex Double) -> Vector Double
+svC = svAux zgesvd "svC" . fmat
+
+-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
+svRd :: Matrix Double -> Vector Double
+svRd = svAux dgesdd "svRd" . fmat
+
+-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
+svCd :: Matrix (Complex Double) -> Vector Double
+svCd = svAux zgesdd "svCd" . fmat
+
+svAux f st x = unsafePerformIO $ do
+    s <- createVector q
+    app2 g mat x vec s st
+    return s
+  where r = rows x
+        c = cols x
+        q = min r c
+        g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr
+
+
+-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
+rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
+rightSVR = rightSVAux dgesvd "rightSVR" . fmat
+
+-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
+rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
+rightSVC = rightSVAux zgesvd "rightSVC" . fmat
+
+rightSVAux f st x = unsafePerformIO $ do
+    s <- createVector q
+    v <- createMatrix ColumnMajor c c
+    app3 g mat x vec s mat v st
+    return (s,trans v)
+  where r = rows x
+        c = cols x
+        q = min r c
+        g ra ca pa = f ra ca pa 0 0 nullPtr
+
+
+-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/  with jobu == \'A\' and jobvt == \'N\'.
+leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
+leftSVR = leftSVAux dgesvd "leftSVR" . fmat
+
+-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
+leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
+leftSVC = leftSVAux zgesvd "leftSVC" . fmat
+
+leftSVAux f st x = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor r r
+    s <- createVector q
+    app3 g mat x mat u vec s st
+    return (u,s)
+  where r = rows x
+        c = cols x
+        q = min r c
+        g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr
+
+-----------------------------------------------------------------------------
+
+foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM
+foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM
+foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM
+foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM
+
+eigAux f st m = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix ColumnMajor r r
+        app3 g mat m vec l mat v st
+        return (l,v)
+  where r = rows m
+        g ra ca pa = f ra ca pa 0 0 nullPtr
+
+
+-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
+-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
+eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
+eigC = eigAux zgeev "eigC" . fmat
+
+eigOnlyAux f st m = unsafePerformIO $ do
+        l <- createVector r
+        app2 g mat m vec l st
+        return l
+  where r = rows m
+        g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr
+
+-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
+-- The eigenvalues are not sorted.
+eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
+eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat
+
+-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.
+-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
+eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
+eigR m = (s', v'')
+    where (s,v) = eigRaux (fmat m)
+          s' = fixeig1 s
+          v' = toRows $ trans v
+          v'' = fromColumns $ fixeig (toList s') v'
+
+eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
+eigRaux m = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix ColumnMajor r r
+        app3 g mat m vec l mat v "eigR"
+        return (l,v)
+  where r = rows m
+        g ra ca pa = dgeev ra ca pa 0 0 nullPtr
+
+fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s))
+    where r = dim s
+
+fixeig  []  _ =  []
+fixeig [_] [v] = [comp' v]
+fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
+    | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1, mapVector negate v2) : fixeig r vs
+    | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs)
+fixeig _ _ = error "fixeig with impossible inputs"
+
+
+-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
+-- The eigenvalues are not sorted.
+eigOnlyR :: Matrix Double -> Vector (Complex Double)
+eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
+
+
+-----------------------------------------------------------------------------
+
+eigSHAux f st m = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix ColumnMajor r r
+        app3 f mat m vec l mat v st
+        return (l,v)
+  where r = rows m
+
+-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.
+-- The eigenvectors are the columns of v.
+-- 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' (fmat m)
+          s' = fromList . reverse . toList $  s
+
+-- | 'eigS' in ascending order
+eigS' :: Matrix Double -> (Vector Double, Matrix Double)
+eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
+
+-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.
+-- The eigenvectors are the columns of v.
+-- 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' (fmat m)
+          s' = fromList . reverse . toList $  s
+
+-- | 'eigH' in ascending order
+eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
+eigH' = eigSHAux (zheev 1) "eigH'" . fmat
+
+
+-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.
+-- The eigenvalues are sorted in descending order.
+eigOnlyS :: Matrix Double -> Vector Double
+eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat
+
+-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.
+-- The eigenvalues are sorted in descending order.
+eigOnlyH :: Matrix (Complex Double) -> Vector Double
+eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat
+
+vrev = flatten . flipud . reshape 1
+
+-----------------------------------------------------------------------------
+foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM
+foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM
+foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM
+foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM
+
+linearSolveSQAux f st a b
+    | n1==n2 && n1==r = unsafePerformIO $ do
+        s <- createMatrix ColumnMajor r c
+        app3 f mat a mat b mat s st
+        return s
+    | otherwise = error $ st ++ " of nonsquare matrix"
+  where n1 = rows a
+        n2 = cols a
+        r  = rows b
+        c  = cols b
+
+-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.
+linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
+linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b)
+
+-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'.
+linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b)
+
+
+-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.
+cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
+cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b)
+
+-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'.
+cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
+cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b)
+
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM
+foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM
+foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM
+foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
+
+linearSolveAux f st a b = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor (max m n) nrhs
+    app3 f mat a mat b mat r st
+    return r
+  where m = rows a
+        n = cols a
+        nrhs = cols b
+
+-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. 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)
+
+-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. 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) $
+                     linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b)
+
+-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. 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) $
+                                   linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b)
+linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b)
+
+-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. 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) $
+                                   linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b)
+linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
+
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM
+foreign import ccall unsafe "chol_l_S" dpotrf :: TMM
+
+cholAux f st a = do
+    r <- createMatrix ColumnMajor n n
+    app2 f mat a mat r st
+    return r
+  where n = rows a
+
+-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
+cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
+cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat
+
+-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
+cholS :: Matrix Double -> Matrix Double
+cholS =  unsafePerformIO . cholAux dpotrf "cholS" . fmat
+
+-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version).
+mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
+mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat
+
+-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/  ('Maybe' version).
+mbCholS :: Matrix Double -> Maybe (Matrix Double)
+mbCholS =  unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat
+
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM
+foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM
+
+-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
+qrR :: Matrix Double -> (Matrix Double, Vector Double)
+qrR = qrAux dgeqr2 "qrR" . fmat
+
+-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.
+qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
+qrC = qrAux zgeqr2 "qrC" . fmat
+
+qrAux f st a = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor m n
+    tau <- createVector mn
+    app3 f mat a vec tau mat r st
+    return (r,tau)
+  where
+    m = rows a
+    n = cols a
+    mn = min m n
+
+foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM
+foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM
+
+-- | build rotation from reflectors
+qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
+qrgrR = qrgrAux dorgqr "qrgrR"
+-- | build rotation from reflectors
+qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
+qrgrC = qrgrAux zungqr "qrgrC"
+
+qrgrAux f st n (a, tau) = unsafePerformIO $ do
+    res <- createMatrix ColumnMajor (rows a) n
+    app3 f mat (fmat a) vec (subVector 0 n tau') mat res st
+    return res
+  where
+    tau' = vjoin [tau, constantD 0 n]
+
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM
+foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM
+
+-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
+hessR :: Matrix Double -> (Matrix Double, Vector Double)
+hessR = hessAux dgehrd "hessR" . fmat
+
+-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.
+hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
+hessC = hessAux zgehrd "hessC" . fmat
+
+hessAux f st a = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor m n
+    tau <- createVector (mn-1)
+    app3 f mat a vec tau mat r st
+    return (r,tau)
+  where m = rows a
+        n = cols a
+        mn = min m n
+
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "schur_l_R" dgees :: TMMM
+foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM
+
+-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
+schurR :: Matrix Double -> (Matrix Double, Matrix Double)
+schurR = schurAux dgees "schurR" . fmat
+
+-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.
+schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
+schurC = schurAux zgees "schurC" . fmat
+
+schurAux f st a = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor n n
+    s <- createMatrix ColumnMajor n n
+    app3 f mat a mat u mat s st
+    return (u,s)
+  where n = rows a
+
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM
+foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM
+
+-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
+luR :: Matrix Double -> (Matrix Double, [Int])
+luR = luAux dgetrf "luR" . fmat
+
+-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
+luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
+luC = luAux zgetrf "luC" . fmat
+
+luAux f st a = unsafePerformIO $ do
+    lu <- createMatrix ColumnMajor n m
+    piv <- createVector (min n m)
+    app3 f mat a vec piv mat lu st
+    return (lu, map (pred.round) (toList piv))
+  where n = rows a
+        m = cols a
+
+-----------------------------------------------------------------------------------
+type TW a = CInt -> PD -> a
+type TQ a = CInt -> CInt -> PC -> a
+
+foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM
+foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt))))
+
+-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/.
+lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
+lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b)
+
+-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.
+lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
+lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b)
+
+lusAux f st a piv b
+    | n1==n2 && n2==n =unsafePerformIO $ do
+         x <- createMatrix ColumnMajor n m
+         app4 f mat a vec piv' mat b mat x st
+         return x
+    | otherwise = error $ st ++ " on LU factorization of nonsquare matrix"
+  where n1 = rows a
+        n2 = cols a
+        n = rows b
+        m = cols b
+        piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double
+