diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,15 @@
+0.8.0.0
+=======
+
+- singularValues, fullSVD, thinSVD, compactSVD, leftSV, rightSV
+  and complete interface to [d|z]gesdd.
+  Algorithms based on the SVD of large matrices can now be
+  significantly faster.
+
+- eigenvalues, eigenvaluesSH
+
+- linearSolveLS, rq
+
 0.7.2.0
 =======
 
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.7.2.1
+Version:            0.8.0.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -7,11 +7,11 @@
 Stability:          provisional
 Homepage:           http://www.hmatrix.googlepages.com
 Synopsis:           Linear algebra and numerical computations
-Description:        This library provides a purely functional interface to basic linear algebra
+Description:        Purely functional interface to basic linear algebra
                     and other numerical computations, internally implemented using
                     GSL, BLAS and LAPACK.
 Category:           Math
-tested-with:        GHC ==6.12.1
+tested-with:        GHC ==6.10.4, GHC ==6.12.1
 
 cabal-version:      >=1.2
 build-type:         Custom
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -10,10 +10,9 @@
 Stability   :  provisional
 Portability :  uses ffi
 
-A generic interface for some common functions. Using it we can write higher level algorithms
-and testing properties both for real and complex matrices.
+Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices.
 
-In any case, the specific functions for particular base types can also be explicitly
+Specific functions for particular base types can also be explicitly
 imported from "Numeric.LinearAlgebra.LAPACK".
 
 -}
@@ -28,17 +27,23 @@
 -- * Linear Systems
     linearSolve,
     luSolve,
+    linearSolveLS,
     linearSolveSVD,
     inv, pinv,
     det, rank, rcond,
 -- * Matrix factorizations
 -- ** Singular value decomposition
     svd,
-    full, economy, --thin,
+    fullSVD,
+    thinSVD,
+    compactSVD,
+    singularValues,
+    leftSV, rightSV,
 -- ** Eigensystems
     eig, eigSH, eigSH',
+    eigenvalues, eigenvaluesSH, eigenvaluesSH',
 -- ** QR
-    qr,
+    qr, rq,
 -- ** Cholesky
     chol, cholSH,
 -- ** Hessenberg
@@ -54,6 +59,7 @@
 -- * Nullspace
     nullspacePrec,
     nullVector,
+    nullspaceSVD,
 -- * Norms
     Normed(..), NormType(..),
 -- * Misc
@@ -63,8 +69,8 @@
     haussholder,
     unpackQR, unpackHess,
     pinvTol,
-    rankSVD, ranksv,
-    nullspaceSVD
+    ranksv,
+    full, economy
 ) where
 
 
@@ -80,12 +86,17 @@
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     svd'         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
+    thinSVD'     :: Matrix t -> (Matrix t, Vector Double, Matrix t)
+    sv'          :: Matrix t -> Vector Double
     luPacked'    :: Matrix t -> (Matrix t, [Int])
     luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix t
     linearSolve' :: Matrix t -> Matrix t -> Matrix t
     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
+    linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t
     eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
     eigSH''      :: Matrix t -> (Vector Double, Matrix t)
+    eigOnly      :: Matrix t -> Vector (Complex Double)
+    eigOnlySH    :: Matrix t -> Vector Double
     cholSH'      :: Matrix t -> Matrix t
     qr'          :: Matrix t -> (Matrix t, Matrix t)
     hess'        :: Matrix t -> (Matrix t, Matrix t)
@@ -94,55 +105,191 @@
     multiply'    :: Matrix t -> Matrix t -> Matrix t
 
 
--- | Singular value decomposition using lapack's dgesvd or zgesvd.
-svd         :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
+instance Field Double where
+    svd' = svdRd
+    thinSVD' = thinSVDRd
+    sv' = svR
+    luPacked' = luR
+    luSolve' (l_u,perm) = lusR l_u perm
+    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
+    linearSolveLS' = linearSolveLSR
+    linearSolveSVD' = linearSolveSVDR Nothing
+    ctrans' = trans
+    eig' = eigR
+    eigSH'' = eigS
+    eigOnly = eigOnlyR
+    eigOnlySH = eigOnlyS
+    cholSH' = cholS
+    qr' = unpackQR . qrR
+    hess' = unpackHess hessR
+    schur' = schurR
+    multiply' = multiplyR
+
+instance Field (Complex Double) where
+    svd' = svdCd
+    thinSVD' = thinSVDCd
+    sv' = svC
+    luPacked' = luC
+    luSolve' (l_u,perm) = lusC l_u perm
+    linearSolve' = linearSolveC
+    linearSolveLS' = linearSolveLSC
+    linearSolveSVD' = linearSolveSVDC Nothing
+    ctrans' = conj . trans
+    eig' = eigC
+    eigOnly = eigOnlyC
+    eigSH'' = eigH
+    eigOnlySH = eigOnlyH
+    cholSH' = cholH
+    qr' = unpackQR . qrC
+    hess' = unpackHess hessC
+    schur' = schurC
+    multiply' = multiplyC
+
+--------------------------------------------------------------
+
+-- | Full singular value decomposition.
+svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
 svd = svd'
 
+-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
+--
+-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.
+thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
+thinSVD = thinSVD'
+
+-- | Singular values only.
+singularValues :: Field t => Matrix t -> Vector Double
+singularValues = sv'
+
+-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
+--
+-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
+fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
+fullSVD m = (u,d,v) where
+    (u,s,v) = svd m
+    d = diagRect s r c
+    r = rows m
+    c = cols m
+
+-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
+compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
+compactSVD m = (u', subVector 0 d s, v') where
+    (u,s,v) = thinSVD m
+    d = rankSVD (1*eps) m s `max` 1
+    u' = takeColumns d u
+    v' = takeColumns d v
+
+vertical m = rows m >= cols m
+
+-- | Singular values and all right singular vectors.
+rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
+rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
+          | otherwise  = let (_,s,v) = svd m     in (s,v)
+
+-- | Singular values and all right singular vectors.
+leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
+leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
+          | otherwise  = let (u,s,_) = thinSVD m in (u,s)
+
+
+{-# DEPRECATED full "use fullSVD instead" #-}
+full svdFun m = (u, d ,v) where
+    (u,s,v) = svdFun m
+    d = diagRect s r c
+    r = rows m
+    c = cols m
+
+{-# DEPRECATED economy "use compactSVD instead" #-}
+economy svdFun m = (u', subVector 0 d s, v') where
+    (u,s,v) = svdFun m
+    d = rankSVD (1*eps) m s `max` 1
+    u' = takeColumns d u
+    v' = takeColumns d v
+
+
+--------------------------------------------------------------
+
 -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
 luPacked    :: Field t => Matrix t -> (Matrix t, [Int])
 luPacked = luPacked'
 
--- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization
---   obtained by 'luPacked'.
+-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
 luSolve     :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
 luSolve = luSolve'
 
--- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.
+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
 -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
---  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
 linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolve = linearSolve'
+
+-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value.
 linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolveSVD = linearSolveSVD'
 
--- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev.
+
+-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
+linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
+linearSolveLS = linearSolveLS'
+
+--------------------------------------------------------------
+
+-- | Eigenvalues and eigenvectors of a general square matrix.
 --
 -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
 eig         :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
 eig = eig'
 
+-- | Eigenvalues of a general square matrix.
+eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
+eigenvalues = eigOnly
+
 -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric.
 eigSH'      :: Field t => Matrix t -> (Vector Double, Matrix t)
 eigSH' = eigSH''
 
--- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
-cholSH      :: Field t => Matrix t -> Matrix t
-cholSH = cholSH'
+-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric.
+eigenvaluesSH' :: Field t => Matrix t -> Vector Double
+eigenvaluesSH' = eigOnlySH
 
--- | QR factorization using lapack's dgeqr2 or zgeqr2.
+-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
 --
+-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
+eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
+eigSH m | m `equal` ctrans m = eigSH' m
+        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
+
+-- | Eigenvalues of a complex hermitian or real symmetric matrix.
+eigenvaluesSH :: Field t => Matrix t -> Vector Double
+eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
+                | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
+
+--------------------------------------------------------------
+
+-- | QR factorization.
+--
 -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
 qr          :: Field t => Matrix t -> (Matrix t, Matrix t)
 qr = qr'
 
--- | Hessenberg factorization using lapack's dgehrd or zgehrd.
+-- | RQ factorization.
 --
+-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
+rq :: Field t => Matrix t -> (Matrix t, Matrix t)
+rq m = (r,q) where
+    (q',r') = qr $ trans $ rev1 m
+    r = rev2 (trans r')
+    q = rev2 (trans q')
+    rev1 = flipud . fliprl
+    rev2 = fliprl . flipud
+
+-- | Hessenberg factorization.
+--
 -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
--- and h is in upper Hessenberg form.
+-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
 hess        :: Field t => Matrix t -> (Matrix t, Matrix t)
 hess = hess'
 
--- | Schur factorization using lapack's dgees or zgees.
+-- | Schur factorization.
 --
 -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
 -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
@@ -162,106 +309,52 @@
 multiply = multiply'
 
 
-instance Field Double where
-    svd' = svdR
-    luPacked' = luR
-    luSolve' (l_u,perm) = lusR l_u perm
-    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
-    linearSolveSVD' = linearSolveSVDR Nothing
-    ctrans' = trans
-    eig' = eigR
-    eigSH'' = eigS
-    cholSH' = cholS
-    qr' = unpackQR . qrR
-    hess' = unpackHess hessR
-    schur' = schurR
-    multiply' = multiplyR
-
-instance Field (Complex Double) where
-    svd' = svdC
-    luPacked' = luC
-    luSolve' (l_u,perm) = lusC l_u perm
-    linearSolve' = linearSolveC
-    linearSolveSVD' = linearSolveSVDC Nothing
-    ctrans' = conj . trans
-    eig' = eigC
-    eigSH'' = eigH
-    cholSH' = cholH
-    qr' = unpackQR . qrC
-    hess' = unpackHess hessC
-    schur' = schurC
-    multiply' = multiplyC
-
-
--- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
---
--- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
-eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
-eigSH m | m `equal` ctrans m = eigSH' m
-        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
+-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
+cholSH      :: Field t => Matrix t -> Matrix t
+cholSH = cholSH'
 
--- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
+-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
 --
 -- If @c = chol m@ then @m == ctrans c \<> c@.
 chol :: Field t => Matrix t ->  Matrix t
 chol m | m `equal` ctrans m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
 
+
+
+
 square m = rows m == cols m
 
--- | determinant of a square matrix, computed from the LU decomposition.
+-- | Determinant of a square matrix.
 det :: Field t => Matrix t -> t
 det m | square m = s * (product $ toList $ takeDiag $ lup)
       | otherwise = error "det of nonsquare matrix"
     where (lup,perm) = luPacked m
           s = signlp (rows m) perm
 
--- | Explicit LU factorization of a general matrix using lapack's dgetrf or zgetrf.
+-- | Explicit LU factorization of a general matrix.
 --
 -- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
 -- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
 lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
 lu = luFact . luPacked
 
--- | Inverse of a square matrix using lapacks' dgesv and zgesv.
+-- | Inverse of a square matrix.
 inv :: Field t => Matrix t -> Matrix t
 inv m | square m = m `linearSolve` ident (rows m)
       | otherwise = error "inv of nonsquare matrix"
 
--- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss.
+-- | Pseudoinverse of a general matrix.
 pinv :: Field t => Matrix t -> Matrix t
 pinv m = linearSolveSVD m (ident (rows m))
 
--- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
---
--- If @(u,d,v) = full svd m@ then @m == u \<> d \<> trans v@.
-full :: Element t 
-     => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
-full svdFun m = (u, d ,v) where
-    (u,s,v) = svdFun m
-    d = diagRect s r c
-    r = rows m
-    c = cols m
-
--- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations.
---
--- If @(u,s,v) = economy svd m@ then @m == u \<> diag s \<> trans v@.
-economy :: Element t 
-        => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
-economy svdFun m = (u', subVector 0 d s, v') where
-    r@(u,s,v) = svdFun m
-    d = rankSVD (1*eps) m r  `max` 1
-    u' = takeColumns d u
-    v' = takeColumns d v
-
-
 -- | Numeric rank of a matrix from the SVD decomposition.
 rankSVD :: Element t
         => Double   -- ^ numeric zero (e.g. 1*'eps')
         -> Matrix t -- ^ input matrix m
-        -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m
+        -> Vector Double -- ^ 'sv' of m
         -> Int      -- ^ rank of m
-rankSVD teps m (_,s,_) = ranksv teps (max (rows m) (cols m)) (toList s)
+rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)
 
 -- | Numeric rank of a matrix from its singular values.
 ranksv ::  Double   -- ^ numeric zero (e.g. 1*'eps')
@@ -316,12 +409,12 @@
 pnormCV Infinity = vectorMax . mapVector magnitude
 --pnormCV _ = error "pnormCV not yet defined"
 
-pnormRM PNorm2 m = head (toList s) where (_,s,_) = svdR m
+pnormRM PNorm2 m = singularValues m @> 0
 pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m
 pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)
 --pnormRM _ _ = error "p norm not yet defined"
 
-pnormCM PNorm2 m = head (toList s) where (_,s,_) = svdC m
+pnormCM PNorm2 m = singularValues m @> 0
 pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m
 pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m)
 --pnormCM _ _ = error "p norm not yet defined"
@@ -354,9 +447,9 @@
              => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
                                   --   or Right \"theoretical\" matrix rank.
              -> Matrix t          -- ^ input matrix m
-             -> (Matrix t, Vector Double, Matrix t) -- ^ 'svd' of m
+             -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
              -> [Vector t]        -- ^ list of unitary vectors spanning the nullspace
-nullspaceSVD hint a z@(_,_,v) = vs where
+nullspaceSVD hint a (s,v) = vs where
     r = rows a
     c = cols a
     tol = case hint of
@@ -364,22 +457,18 @@
         _      -> eps
     k = case hint of
         Right t -> t
-        _       -> rankSVD tol a z
-    nsol = c - k
-    vs = drop k (toColumns v)
+        _       -> rankSVD tol a s
+    vs = drop k $ toRows $ ctrans v
 
 
 -- | The nullspace of a matrix. See also 'nullspaceSVD'.
 nullspacePrec :: Field t
-              => Double     -- ^ relative tolerance in 'eps' units
+              => Double     -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
               -> Matrix t   -- ^ input matrix
               -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
-nullspacePrec t m = ns where
-    r@(_,_,v) = svd m
-    k = rankSVD (t*eps) m r
-    ns = drop k $ toRows $ ctrans v
+nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
 
--- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'.
+-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
 nullVector :: Field t => Matrix t -> Vector t
 nullVector = last . nullspacePrec 1
 
@@ -405,7 +494,7 @@
 -}
 --pinvTol :: Double -> Matrix Double -> Matrix Double
 pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
-    (u,s,v) = svdR m
+    (u,s,v) = thinSVDRd m
     sl@(g:_) = toList s
     s' = fromList . map rec $ sl
     rec x = if x < g*tol then 1 else 1/x
@@ -460,15 +549,14 @@
 
 --------------------------------------------------------------------------
 
--- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD.
+-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
 rcond :: Field t => Matrix t -> Double
 rcond m = last s / head s
-    where (_,s',_) = svd m
-          s = toList s'
+    where s = toList (singularValues m)
 
 -- | Number of linearly independent rows or columns.
 rank :: Field t => Matrix t -> Int
-rank m = rankSVD eps m (svd m)
+rank m = rankSVD eps m (singularValues m)
 
 {-
 expm' m = case diagonalize (complex m) of
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -190,8 +190,8 @@
 
 ---------------------------------------------------------------
 
-instance (Storable a) => Monoid (Vector a) where
-    mempty = V { dim = 0, fptr = undefined }
+instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
+    mempty = 0 { dim = 0 }
     mappend a b = mconcat [a,b]
     mconcat = j . filter ((>0).dim)
         where j [] = mempty
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -1,4 +1,3 @@
-{-# OPTIONS_GHC #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.LinearAlgebra.LAPACK
@@ -9,21 +8,35 @@
 -- Stability   :  provisional
 -- Portability :  portable (uses FFI)
 --
--- Wrappers for a few LAPACK functions (<http://www.netlib.org/lapack>).
+-- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).
 --
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.LAPACK (
+    -- * Matrix product
     multiplyR, multiplyC,
-    svdR, svdRdd, svdC,
-    eigC, eigR, eigS, eigH, eigS', eigH',
+    -- * Linear systems
     linearSolveR, linearSolveC,
+    lusR, lusC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
-    luR, luC, lusR, lusC,
+    -- * 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,
+    -- * QR
     qrR, qrC,
+    -- * Hessenberg
     hessR, hessC,
+    -- * Schur
     schurR, schurC
 ) where
 
@@ -61,28 +74,27 @@
 multiplyC a b = multiplyAux zgemmc "zgemmc" a b
 
 -----------------------------------------------------------------------------
-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
+foreign import ccall "svd_l_R" dgesvd :: TMMVM
+foreign import ccall "svd_l_C" zgesvd :: TCMCMVCM
+foreign import ccall "svd_l_Rdd" dgesdd :: TMMVM
+foreign import ccall "svd_l_Cdd" zgesdd :: TCMCMVCM
 
--- | 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@.
+-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
 svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
 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 = svdAux dgesdd "svdRdd" . 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
 
--- | 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@.
+-- | 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)
@@ -92,59 +104,152 @@
   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
+
 -----------------------------------------------------------------------------
-eigAux f st m
-    | r == 1 = (fromList [flatten m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
+
+foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: CInt -> TMVM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: CInt -> TCMVCM
+
+eigAux f st m = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
-        dummy <- createMatrix ColumnMajor 1 1
-        app4 f mat m mat dummy vec l mat v st
+        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
 
 
-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:
---
--- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@.
---
--- The eigenvectors are the columns of v.
--- The eigenvalues are not sorted.
+-- | 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
 
--- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
---
--- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@.
---
--- The eigenvectors are the columns of v.
+-- | 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' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
+          s' = fixeig1 s
           v' = toRows $ trans v
           v'' = fromColumns $ fixeig (toList s') v'
           r = rows m
 
 eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
-eigRaux m
-    | r == 1 = (fromList [(flatten m `at` 0):+0], singleton 1)
-    | otherwise = unsafePerformIO $ do
+eigRaux m = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
-        dummy <- createMatrix ColumnMajor 1 1
-        app4 dgeev mat m mat dummy vec l mat v "eigR"
+        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)
@@ -153,50 +258,59 @@
   where scale = vectorMapValR Scale
 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
+
+
 -----------------------------------------------------------------------------
 
--- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
---
--- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@.
---
+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).
+-- 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' m
-    | r == 1 = (fromList [flatten m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        app3 dsyev mat m vec l mat v "eigS"
-        return (l,v)
-  where r = rows m
-
------------------------------------------------------------------------------
+-- | 'eigS' in ascending order
+eigS' :: Matrix Double -> (Vector Double, Matrix Double)
+eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
 
--- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix:
---
--- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@.
---
+-- | 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).
+-- 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' m
-    | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        app3 zheev mat m vec l mat v "eigH"
-        return (l,v)
-  where r = rows m
+-- | '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 "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
 foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
@@ -212,11 +326,11 @@
         r  = rows b
         c  = cols b
 
--- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition.
+-- | 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)
 
--- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition.
+-- | 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)
 
@@ -234,17 +348,17 @@
         n = cols a
         nrhs = cols b
 
--- | 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'.
+-- | 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)
 
--- | 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'.
+-- | 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)
 
--- | 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.
+-- | 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)
@@ -253,7 +367,7 @@
                                    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.
+-- | 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)
@@ -266,13 +380,11 @@
 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.
+-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
 cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
 cholH = cholAux zpotrf "cholH" . fmat
 
--- | Wrapper for LAPACK's /dpotrf/, which computes the Cholesky factorization of a
--- real symmetric positive definite matrix.
+-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
 cholS :: Matrix Double -> Matrix Double
 cholS = cholAux dpotrf "cholS" . fmat
 
@@ -286,11 +398,11 @@
 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.
+-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
 qrR :: Matrix Double -> (Matrix Double, Vector Double)
 qrR = qrAux dgeqr2 "qrR" . fmat
 
--- | Wrapper for LAPACK's /zgeqr2/, which computes a QR factorization of a complex matrix.
+-- | 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
 
@@ -307,11 +419,11 @@
 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.
+-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
 hessR :: Matrix Double -> (Matrix Double, Vector Double)
 hessR = hessAux dgehrd "hessR" . fmat
 
--- | Wrapper for LAPACK's /zgehrd/, which computes a Hessenberg factorization of a square complex matrix.
+-- | 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
 
@@ -328,11 +440,11 @@
 foreign import ccall "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.
+-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
 schurR :: Matrix Double -> (Matrix Double, Matrix Double)
 schurR = schurAux dgees "schurR" . fmat
 
--- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix.
+-- | 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
 
@@ -347,11 +459,11 @@
 foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM
 foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM
 
--- | Wrapper for LAPACK's /dgetrf/, which computes a LU factorization of a general real matrix.
+-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
 luR :: Matrix Double -> (Matrix Double, [Int])
 luR = luAux dgetrf "luR" . fmat
 
--- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix.
+-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
 luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
 luC = luAux zgetrf "luC" . fmat
 
@@ -370,11 +482,11 @@
 foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM
 foreign import ccall "LAPACK/lapack-aux.h luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt))))
 
--- | Wrapper for LAPACK's /dgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition.
+-- | 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)
 
--- | Wrapper for LAPACK's /zgetrs/, which solves a general real linear system (for several right-hand sides) from a precomputed LU decomposition.
+-- | 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)
 
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -48,7 +48,27 @@
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
-    REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES(up==NULL || ur==m && (uc==m || uc==q),BAD_SIZE);
+    char* jobu  = "A";
+    if (up==NULL) {
+        jobu = "N";
+    } else {
+        if (uc==q) {
+            jobu = "S";
+        }
+    }
+    REQUIRES(vp==NULL || vc==n && (vr==n || vr==q),BAD_SIZE);
+    char* jobvt  = "A";
+    integer ldvt = n;
+    if (vp==NULL) {
+        jobvt = "N";
+    } else {
+        if (vr==q) {
+            jobvt = "S";
+            ldvt = q;
+        }
+    }
     DEBUGMSG("svd_l_R");
     double *B = (double*)malloc(m*n*sizeof(double));
     CHECK(!B,MEM);
@@ -57,25 +77,21 @@
     integer res;
     // ask for optimal lwork
     double ans;
-    //printf("ask zgesvd\n");
-    char* job = "A";
-    dgesvd_ (job,job,
+    dgesvd_ (jobu,jobvt,
              &m,&n,B,&m,
              sp,
              up,&m,
-             vp,&n,
+             vp,&ldvt,
              &ans, &lwork,
              &res);
     lwork = ceil(ans);
-    //printf("ans = %d\n",lwork);
     double * work = (double*)malloc(lwork*sizeof(double));
     CHECK(!work,MEM);
-    //printf("dgesdd\n");
-    dgesvd_ (job,job,
+    dgesvd_ (jobu,jobvt,
              &m,&n,B,&m,
              sp,
              up,&m,
-             vp,&n,
+             vp,&ldvt,
              work, &lwork,
              &res);
     CHECK(res,res);
@@ -90,25 +106,36 @@
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
-    REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES(up == NULL && vp == NULL
+             || ur==m && vc==n
+                &&   (uc == q && vr == q
+                   || uc == m && vc==n),BAD_SIZE);
+    char* jobz  = "A";
+    integer ldvt = n;
+    if (up==NULL) {
+        jobz = "N";
+    } else {
+        if (uc==q && vr == q) {
+            jobz = "S";
+            ldvt = q;
+        }
+    }
     DEBUGMSG("svd_l_Rdd");
     double *B = (double*)malloc(m*n*sizeof(double));
     CHECK(!B,MEM);
     memcpy(B,ap,m*n*sizeof(double));
-    integer* iwk = (integer*) malloc(8*q*sizeof(int));
+    integer* iwk = (integer*) malloc(8*q*sizeof(integer));
     CHECK(!iwk,MEM);
     integer lwk = -1;
     integer res;
     // ask for optimal lwk
     double ans;
-    //printf("ask dgesdd\n");
-    dgesdd_ ("A",&m,&n,B,&m,sp,up,&m,vp,&n,&ans,&lwk,iwk,&res);
-    lwk = 2*ceil(ans); // ????? otherwise 50x100 rejects lwk
-    //printf("lwk = %d\n",lwk);
+    dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
+    lwk = ans;
     double * workv = (double*)malloc(lwk*sizeof(double));
     CHECK(!workv,MEM);
-    //printf("dgesdd\n");
-    dgesdd_ ("A",&m,&n,B,&m,sp,up,&m,vp,&n,workv,&lwk,iwk,&res);
+    dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
     CHECK(res,res);
     free(iwk);
     free(workv);
@@ -120,17 +147,36 @@
 
 // not in clapack.h
 
-int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, 
-    doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, 
-    integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, 
+int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
+    doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
+    integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
     integer *lwork, doublereal *rwork, integer *info);
 
 int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
-    REQUIRES(ur==m && uc==m && sn==q && vr==n && vc==n,BAD_SIZE);
-    DEBUGMSG("svd_l_C");
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES(up==NULL || ur==m && (uc==m || uc==q),BAD_SIZE);
+    char* jobu  = "A";
+    if (up==NULL) {
+        jobu = "N";
+    } else {
+        if (uc==q) {
+            jobu = "S";
+        }
+    }
+    REQUIRES(vp==NULL || vc==n && (vr==n || vr==q),BAD_SIZE);
+    char* jobvt  = "A";
+    integer ldvt = n;
+    if (vp==NULL) {
+        jobvt = "N";
+    } else {
+        if (vr==q) {
+            jobvt = "S";
+            ldvt = q;
+        }
+    }DEBUGMSG("svd_l_C");
     double *B = (double*)malloc(2*m*n*sizeof(double));
     CHECK(!B,MEM);
     memcpy(B,ap,m*n*2*sizeof(double));
@@ -141,26 +187,22 @@
     integer res;
     // ask for optimal lwork
     doublecomplex ans;
-    //printf("ask zgesvd\n");
-    char* job = "A";
-    zgesvd_ (job,job,
+    zgesvd_ (jobu,jobvt,
              &m,&n,(doublecomplex*)B,&m,
              sp,
              (doublecomplex*)up,&m,
-             (doublecomplex*)vp,&n,
+             (doublecomplex*)vp,&ldvt,
              &ans, &lwork,
              rwork,
              &res);
     lwork = ceil(ans.r);
-    //printf("ans = %d\n",lwork);
     doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));
     CHECK(!work,MEM);
-    //printf("zgesvd\n");
-    zgesvd_ (job,job,
+    zgesvd_ (jobu,jobvt,
              &m,&n,(doublecomplex*)B,&m,
              sp,
              (doublecomplex*)up,&m,
-             (doublecomplex*)vp,&n,
+             (doublecomplex*)vp,&ldvt,
              work, &lwork,
              rwork,
              &res);
@@ -171,23 +213,81 @@
     OK
 }
 
+int zgesdd_ (char *jobz, integer *m, integer *n,
+    doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
+    integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
+    integer *lwork, doublereal *rwork, integer* iwork, integer *info);
 
+int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
+    //printf("entro\n");
+    integer m = ar;
+    integer n = ac;
+    integer q = MIN(m,n);
+    REQUIRES(sn==q,BAD_SIZE);
+    REQUIRES(up == NULL && vp == NULL
+             || ur==m && vc==n
+                &&   (uc == q && vr == q
+                   || uc == m && vc==n),BAD_SIZE);
+    char* jobz  = "A";
+    integer ldvt = n;
+    if (up==NULL) {
+        jobz = "N";
+    } else {
+        if (uc==q && vr == q) {
+            jobz = "S";
+            ldvt = q;
+        }
+    }
+    DEBUGMSG("svd_l_Cdd");
+    doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
+    CHECK(!B,MEM);
+    memcpy(B,ap,m*n*sizeof(doublecomplex));
+    integer* iwk = (integer*) malloc(8*q*sizeof(integer));
+    CHECK(!iwk,MEM);
+    int lrwk;
+    if (0 && *jobz == 'N') {
+        lrwk = 5*q; // does not work, crash at free below
+    } else {
+        lrwk = 5*q*q + 7*q;
+    }
+    double *rwk = (double*)malloc(lrwk*sizeof(double));;
+    CHECK(!rwk,MEM);
+    //printf("%s %ld %d\n",jobz,q,lrwk);
+    integer lwk = -1;
+    integer res;
+    // ask for optimal lwk
+    doublecomplex ans;
+    zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
+    lwk = ans.r;
+    //printf("lwk = %ld\n",lwk);
+    doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
+    CHECK(!workv,MEM);
+    zgesdd_ (jobz,&m,&n,B,&m,sp,(doublecomplex*)up,&m,(doublecomplex*)vp,&ldvt,workv,&lwk,rwk,iwk,&res);
+    //printf("res = %ld\n",res);
+    CHECK(res,res);
+    free(workv); // printf("freed workv\n");
+    free(rwk);   // printf("freed rwk\n");
+    free(iwk);   // printf("freed iwk\n");
+    free(B);     // printf("freed B, salgo\n");
+    OK
+}
 
 //////////////////// general complex eigensystem ////////////
 
-int eig_l_C(KCMAT(a),CMAT(u), CVEC(s),CMAT(v)) {
+int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac==n && (ur==1 || (ur==n && uc==n)) && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(up==NULL || ur==n && uc==n, BAD_SIZE);
+    char jobvl = up==NULL?'N':'V';
+    REQUIRES(vp==NULL || vr==n && vc==n, BAD_SIZE);
+    char jobvr = vp==NULL?'N':'V';
     DEBUGMSG("eig_l_C");
-    double *B = (double*)malloc(2*n*n*sizeof(double));
+    doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
     CHECK(!B,MEM);
-    memcpy(B,ap,n*n*2*sizeof(double));
-
+    memcpy(B,ap,n*n*sizeof(doublecomplex));
     double *rwork = (double*) malloc(2*n*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
-    char jobvl = ur==1?'N':'V';
-    char jobvr = vr==1?'N':'V';
     integer res;
     // ask for optimal lwork
     doublecomplex ans;
@@ -202,7 +302,7 @@
              &res);
     lwork = ceil(ans.r);
     //printf("ans = %d\n",lwork);
-    doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
     //printf("zgeev\n");
     zgeev_  (&jobvl,&jobvr,
@@ -226,14 +326,16 @@
 
 int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac == n && (ur==1 || (ur==n && uc==n)) && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(up==NULL || ur==n && uc==n, BAD_SIZE);
+    char jobvl = up==NULL?'N':'V';
+    REQUIRES(vp==NULL || vr==n && vc==n, BAD_SIZE);
+    char jobvr = vp==NULL?'N':'V';
     DEBUGMSG("eig_l_R");
     double *B = (double*)malloc(n*n*sizeof(double));
     CHECK(!B,MEM);
     memcpy(B,ap,n*n*sizeof(double));
     integer lwork = -1;
-    char jobvl = ur==1?'N':'V';
-    char jobvr = vr==1?'N':'V';
     integer res;
     // ask for optimal lwork
     double ans;
@@ -267,13 +369,14 @@
 //////////////////// symmetric real eigensystem ////////////
 
 
-int eig_l_S(KDMAT(a),DVEC(s),DMAT(v)) {
+int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac == n && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(vr==n && vc==n, BAD_SIZE);
+    char jobz = wantV?'V':'N';
     DEBUGMSG("eig_l_S");
     memcpy(vp,ap,n*n*sizeof(double));
     integer lwork = -1;
-    char jobz = vr==1?'N':'V';
     char uplo = 'U';
     integer res;
     // ask for optimal lwork
@@ -300,15 +403,16 @@
 
 //////////////////// hermitian complex eigensystem ////////////
 
-int eig_l_H(KCMAT(a),DVEC(s),CMAT(v)) {
+int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac==n && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(vr==n && vc==n, BAD_SIZE);
+    char jobz = wantV?'V':'N';
     DEBUGMSG("eig_l_H");
     memcpy(vp,ap,2*n*n*sizeof(double));
     double *rwork = (double*) malloc((3*n-2)*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
-    char jobz = vr==1?'N':'V';
     char uplo = 'U';
     integer res;
     // ask for optimal lwork
@@ -322,7 +426,7 @@
              &res);
     lwork = ceil(ans.r);
     //printf("ans = %d\n",lwork);
-    doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
     zheev_  (&jobz,&uplo,
              &n,(doublecomplex*)vp,&n,
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -68,8 +68,8 @@
 int eig_l_C(KCMAT(a),CMAT(u),CVEC(s),CMAT(v));
 int eig_l_R(KDMAT(a),DMAT(u),CVEC(s),DMAT(v));
 
-int eig_l_S(KDMAT(a),DVEC(s),DMAT(v));
-int eig_l_H(KCMAT(a),DVEC(s),CMAT(v));
+int eig_l_S(int,KDMAT(a),DVEC(s),DMAT(v));
+int eig_l_H(int,KCMAT(a),DVEC(s),CMAT(v));
 
 int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x));
 int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x));
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -22,6 +22,7 @@
 ) where
 
 import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.LAPACK
 import Numeric.LinearAlgebra.Tests.Instances
 import Numeric.LinearAlgebra.Tests.Properties
 import Test.HUnit hiding ((~:),test,Testable)
@@ -177,6 +178,11 @@
     putStrLn "------ luSolve"
     test (linearSolveProp (luSolve.luPacked) . rSqWC)
     test (linearSolveProp (luSolve.luPacked) . cSqWC)
+    putStrLn "------ luSolveLS"
+    test (linearSolveProp linearSolveLS . rSqWC)
+    test (linearSolveProp linearSolveLS . cSqWC)
+    test (linearSolveProp2 linearSolveLS . rConsist)
+    test (linearSolveProp2 linearSolveLS . cConsist)
     putStrLn "------ pinv (linearSolveSVD)"
     test (pinvProp  . rM)
     test (pinvProp  . cM)
@@ -186,19 +192,41 @@
     putStrLn "------ svd"
     test (svdProp1  . rM)
     test (svdProp1  . cM)
-    test (svdProp2  . rM)
-    test (svdProp2  . cM)
+    test (svdProp1a svdR)
+    test (svdProp1a svdC)
+    test (svdProp1a svdRd)
+    test (svdProp1a svdCd)
+    test (svdProp2 thinSVDR)
+    test (svdProp2 thinSVDC)
+    test (svdProp2 thinSVDRd)
+    test (svdProp2 thinSVDCd)
+    test (svdProp3  . rM)
+    test (svdProp3  . cM)
+    test (svdProp4  . rM)
+    test (svdProp4  . cM)
+    test (svdProp5a)
+    test (svdProp5b)
+    test (svdProp6a)
+    test (svdProp6b)
+    test (svdProp7  . rM)
+    test (svdProp7  . cM)
     putStrLn "------ eig"
     test (eigSHProp . rHer)
     test (eigSHProp . cHer)
     test (eigProp   . rSq)
     test (eigProp   . cSq)
+    test (eigSHProp2 . rHer)
+    test (eigSHProp2 . cHer)
+    test (eigProp2   . rSq)
+    test (eigProp2   . cSq)
     putStrLn "------ nullSpace"
     test (nullspaceProp . rM)
     test (nullspaceProp . cM)
     putStrLn "------ qr"
     test (qrProp     . rM)
     test (qrProp     . cM)
+    test (rqProp     . rM)
+    test (rqProp     . cM)
     putStrLn "------ hess"
     test (hessProp   . rSq)
     test (hessProp   . cSq)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -143,8 +143,8 @@
             r = rows m
             c = cols m
             n = min r c
-        sv <- replicateM n (choose (1,100))
-        let s = diagRect (fromList sv) r c
+        sv' <- replicateM n (choose (1,100))
+        let s = diagRect (fromList sv') r c
         return $ WC (u <> real s <> trans v)
 
 #if MIN_VERSION_QuickCheck(2,0,0)
@@ -160,8 +160,8 @@
         Sq m <- arbitrary
         let (u,_,v) = svd m
             n = rows m
-        sv <- replicateM n (choose (1,100))
-        let s = diag (fromList sv)
+        sv' <- replicateM n (choose (1,100))
+        let s = diag (fromList sv')
         return $ SqWC (u <> real s <> trans v)
 
 #if MIN_VERSION_QuickCheck(2,0,0)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -29,23 +29,26 @@
     pinvProp,
     detProp,
     nullspaceProp,
-    svdProp1, svdProp2,
-    eigProp, eigSHProp,
-    qrProp,
+    svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4,
+    svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
+    eigProp, eigSHProp, eigProp2, eigSHProp2,
+    qrProp, rqProp,
     hessProp,
     schurProp1, schurProp2,
     cholProp,
     expmDiagProp,
     multProp1, multProp2,
-    linearSolveProp
+    linearSolveProp, linearSolveProp2
 ) where
 
 import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra.LAPACK
+import Debug.Trace
 #include "quickCheckCompat.h"
--- import Debug.Trace
 
--- debug x = trace (show x) x
 
+debug x = trace (show x) x
+
 -- relative error
 dist :: (Normed t, Num t) => t -> t -> Double
 dist a b = r
@@ -71,8 +74,11 @@
 
 square m = rows m == cols m
 
-unitary m = square m && m <> ctrans m |~| ident (rows m)
+-- orthonormal columns
+orthonormal m = ctrans m <> m |~| ident (cols m)
 
+unitary m = square m && orthonormal m
+
 hermitian m = square m && m |~| ctrans m
 
 wellCond m = rcond m > 1/100
@@ -119,13 +125,66 @@
           r = rows m
           c = cols m - rank m
 
-svdProp1 m = u <> real d <> trans v |~| m
-          && unitary u && unitary v
-    where (u,d,v) = full svd m
+------------------------------------------------------------------
 
-svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m)
-    where (u,s,v) = economy svd m
+-- fullSVD
+svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v
+    where (u,d,v) = fullSVD m
 
+svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where
+    (u,s,v) = svdfun m
+    d = diagRect s (rows m) (cols m)
+
+-- thinSVD
+svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m)
+    where (u,s,v) = thinSVDfun m
+
+-- compactSVD
+svdProp3 m = (m |~| u <> real (diag s) <> trans v
+             && orthonormal u && orthonormal v)
+    where (u,s,v) = compactSVD m
+
+svdProp4 m' = m |~| u <> real (diag s) <> trans v
+           && orthonormal u && orthonormal v
+           && (dim s == r || r == 0 && dim s == 1)
+    where (u,s,v) = compactSVD m
+          m = m' <-> m'
+          r = rank m'
+
+svdProp5a m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where
+    s1       = svR  m
+    s2       = svRd m
+    (_,s3,_) = svdR m
+    (_,s4,_) = svdRd m
+    (_,s5,_) = thinSVDR m
+    (_,s6,_) = thinSVDRd m
+
+svdProp5b m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where
+    s1       = svC  m
+    s2       = svCd m
+    (_,s3,_) = svdC m
+    (_,s4,_) = svdCd m
+    (_,s5,_) = thinSVDC m
+    (_,s6,_) = thinSVDCd m
+
+svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
+    where (u,s,v) = svdR m
+          (s',v') = rightSVR m
+          (u',s'') = leftSVR m
+
+svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
+    where (u,s,v) = svdC m
+          (s',v') = rightSVC m
+          (u',s'') = leftSVC m
+
+svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
+    where (u,s,v) = svd m
+          (s',v') = rightSV m
+          (u',s'') = leftSV m
+          s''' = singularValues m
+
+------------------------------------------------------------------
+
 eigProp m = complex m <> v |~| v <> diag s
     where (s, v) = eig m
 
@@ -134,9 +193,25 @@
               && m |~| v <> real (diag s) <> ctrans v
     where (s, v) = eigSH m
 
+eigProp2 m = fst (eig m) |~| eigenvalues m
+
+eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
+
+------------------------------------------------------------------
+
 qrProp m = q <> r |~| m && unitary q && upperTriang r
     where (q,r) = qr m
 
+rqProp m = r <> q |~| m && unitary q && utr
+    where (r,q) = rq m
+          upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
+              where t = f-c
+          utr = upptr (rows r) (cols r) * r |~| r
+
+upperTriang' m = rows m == 1 || down |~| z
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
 hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
     where (p,h) = hess m
 
@@ -162,3 +237,8 @@
 multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a
 
 linearSolveProp f m = f m m |~| ident (rows m)
+
+linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
+    where q = min (rows a) (cols a)
+          b = a <> x
+          wc = rank a == q