diff --git a/lapack.cabal b/lapack.cabal
--- a/lapack.cabal
+++ b/lapack.cabal
@@ -1,5 +1,5 @@
 Name:             lapack
-Version:          0.0
+Version:          0.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -9,6 +9,9 @@
 Synopsis:         Numerical Linear Algebra using LAPACK
 Description:
   This is a high-level interface to LAPACK.
+  It provides solvers for simultaneous linear equations,
+  linear least-squares problems, eigenvalue and singular value problems
+  for matrices with certain kinds of structures.
   .
   Features:
   .
@@ -34,7 +37,7 @@
 Build-Type:       Simple
 
 Source-Repository this
-  Tag:         0.0
+  Tag:         0.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/lapack/
 
@@ -46,8 +49,8 @@
   Build-Depends:
     lapack-ffi >=0.0.1 && <0.1,
     blas-ffi >=0.0 && <0.1,
-    netlib-ffi >=0.0.1 && <0.1,
-    comfort-array >=0.0 && <0.1,
+    netlib-ffi >=0.1 && <0.2,
+    comfort-array >=0.0.1 && <0.1,
     transformers >=0.3 && <0.6,
     non-empty >=0.3 && <0.4,
     utility-ht >=0.0.10 && <0.1,
@@ -58,9 +61,23 @@
   Exposed-Modules:
     Numeric.LAPACK.Matrix
     Numeric.LAPACK.Matrix.Shape
+    Numeric.LAPACK.Matrix.Square
+    Numeric.LAPACK.Matrix.Hermitian
+    Numeric.LAPACK.Matrix.Triangular
     Numeric.LAPACK.Vector
-    Numeric.LAPACK.LinearSystem
+    Numeric.LAPACK.Orthogonal
+    Numeric.LAPACK.Linear.General
+    Numeric.LAPACK.Linear.HermitianPositiveDefinite
+    Numeric.LAPACK.Linear.Hermitian
+    Numeric.LAPACK.Linear.Triangular
+    Numeric.LAPACK.Eigen.General
+    Numeric.LAPACK.Eigen.Hermitian
+    Numeric.LAPACK.Eigen.Triangular
+    Numeric.LAPACK.Singular
   Other-Modules:
+    Numeric.LAPACK.Matrix.Triangular.Private
     Numeric.LAPACK.Matrix.Shape.Private
+    Numeric.LAPACK.Matrix.Multiply
+    Numeric.LAPACK.Matrix.Private
     Numeric.LAPACK.Private
     Numeric.LAPACK.Format
diff --git a/src/Numeric/LAPACK/Eigen/General.hs b/src/Numeric/LAPACK/Eigen/General.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Eigen/General.hs
@@ -0,0 +1,319 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Eigen.General (
+   values,
+   schur,
+   decompose,
+   ComplexOf,
+   ) where
+
+import Numeric.LAPACK.Matrix.Square (Square)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private
+         (ComplexOf, RealOf, zero, withAutoWorkspaceInfo,
+          copyToTemp, copyToColumnMajor, allocArray)
+
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
+import qualified Numeric.LAPACK.FFI.Real as LapackReal
+import qualified Numeric.BLAS.FFI.Complex as BlasComplex
+import qualified Numeric.BLAS.FFI.Real as BlasReal
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Foreign.Marshal.Array (advancePtr, peekArray)
+import Foreign.C.Types (CInt, CChar)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr, nullPtr, nullFunPtr, castPtr)
+import Foreign.Storable (Storable)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+import Data.Complex (Complex)
+
+
+values ::
+   (Shape.C sh, Class.Floating a) =>
+   Square sh a -> Vector sh (ComplexOf a)
+values =
+   getValues $
+   Class.switchFloating
+      (Values valuesAux) (Values valuesAux)
+      (Values valuesAux) (Values valuesAux)
+
+type Values_ sh a = Square sh a -> Vector sh (ComplexOf a)
+
+newtype Values sh a = Values {getValues :: Values_ sh a}
+
+valuesAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Values_ sh a
+valuesAux (Array (MatrixShape.Square _order size) a) =
+      Array.unsafeCreateWithSize size $ \n wPtr -> do
+   let lda = n
+   evalContT $ do
+      jobvsPtr <- Call.char 'N'
+      sortPtr <- Call.char 'N'
+      aPtr <- copyToTemp (n*n) a
+      ldaPtr <- Call.cint lda
+      sdimPtr <- Call.alloca
+      let vsPtr = nullPtr
+      ldvsPtr <- Call.cint n
+      let bworkPtr = nullPtr
+      liftIO $ withAutoWorkspaceInfo "gees" $ \workPtr lworkPtr infoPtr ->
+         gees
+            jobvsPtr sortPtr n aPtr ldaPtr sdimPtr
+            wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr
+
+
+{- |
+If @(q,r) = schur a@, then @a = q \<#\> r \<#\> adjoint q@,
+where @q@ is unitary (orthogonal)
+and @r@ is a right-upper triangular matrix for complex @a@
+and a 1x1-or-2x2-block upper triangular matrix for real @a@.
+With @getDiagonal r@ you get all eigenvalues of @a@ if @a@ is complex
+and the real parts of the eigenvalues if @a@ is real.
+Complex conjugated eigenvalues of a real matrix @a@
+are encoded as 2x2 blocks along the diagonal.
+-}
+schur ::
+   (Shape.C sh, Class.Floating a) =>
+   Square sh a -> (Square sh a, Square sh a)
+schur =
+   getSchur $
+   Class.switchFloating
+      (Schur schurAux) (Schur schurAux)
+      (Schur schurAux) (Schur schurAux)
+
+type Schur_ sh a = Square sh a -> (Square sh a, Square sh a)
+
+newtype Schur sh a = Schur {getSchur :: Schur_ sh a}
+
+schurAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Schur_ sh a
+schurAux (Array (MatrixShape.Square order size) a) = unsafePerformIO $ do
+   let n = Shape.size size
+   let lda = n
+   let sh = MatrixShape.Square ColumnMajor size
+   evalContT $ do
+      jobvsPtr <- Call.char 'V'
+      sortPtr <- Call.char 'N'
+      aPtr <- ContT $ withForeignPtr a
+      (s,sPtr) <- allocArray sh
+      liftIO $ copyToColumnMajor order n n aPtr sPtr
+      ldaPtr <- Call.cint lda
+      sdimPtr <- Call.alloca
+      wPtr <- Call.allocaArray n
+      (vs,vsPtr) <- allocArray sh
+      ldvsPtr <- Call.cint n
+      let bworkPtr = nullPtr
+      liftIO $ withAutoWorkspaceInfo "gees" $ \workPtr lworkPtr infoPtr ->
+         gees
+            jobvsPtr sortPtr n sPtr ldaPtr sdimPtr
+            wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr
+      return (vs, s)
+
+
+
+type GEES_ ar a =
+   Ptr CChar -> Ptr CChar -> Int -> Ptr a -> Ptr CInt ->
+   Ptr CInt -> Ptr (Complex ar) -> Ptr a -> Ptr CInt ->
+   Ptr a -> Ptr CInt -> Ptr Bool -> Ptr CInt -> IO ()
+
+newtype GEES a = GEES {getGEES :: GEES_ (RealOf a) a}
+
+gees :: Class.Floating a => GEES_ (RealOf a) a
+gees =
+   getGEES $
+   Class.switchFloating
+      (GEES geesReal) (GEES geesReal) (GEES geesComplex) (GEES geesComplex)
+
+geesReal :: Class.Real a => GEES_ a a
+geesReal
+      jobvsPtr sortPtr n aPtr ldaPtr sdimPtr
+      wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr =
+   evalContT $ do
+      let selectPtr = nullFunPtr
+      nPtr <- Call.cint n
+      wrPtr <- Call.allocaArray n
+      wiPtr <- Call.allocaArray n
+      liftIO $
+         LapackReal.gees
+            jobvsPtr sortPtr selectPtr nPtr aPtr ldaPtr sdimPtr
+            wrPtr wiPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr
+      liftIO $ zipComplex n wrPtr wiPtr wPtr
+
+geesComplex :: Class.Real a => GEES_ a (Complex a)
+geesComplex
+      jobvsPtr sortPtr n aPtr ldaPtr sdimPtr
+      wPtr vsPtr ldvsPtr workPtr lworkPtr bworkPtr infoPtr =
+   evalContT $ do
+      let selectPtr = nullFunPtr
+      nPtr <- Call.cint n
+      rworkPtr <- Call.allocaArray n
+      liftIO $
+         LapackComplex.gees
+            jobvsPtr sortPtr selectPtr nPtr aPtr ldaPtr sdimPtr
+            wPtr vsPtr ldvsPtr workPtr lworkPtr rworkPtr bworkPtr infoPtr
+
+
+
+{- |
+@(vr,d,vl) = Eigen.decompose a@
+
+Counterintuitively, @vr@ contains the right eigenvectors
+and @vl@ contains the left eigenvectors as columns.
+The idea is to provide a decomposition of @a@.
+If @a@ is diagonalizable, then @vr@ and @vl@ are almost inverse to each other.
+More precisely, @adjoint vl \<#\> vr@ is a diagonal matrix.
+This is because all eigenvectors are normalized to Euclidean norm 1.
+With the following scaling, the decomposition becomes perfect:
+
+> let scal = Array.map recip $ getDiagonal $ adjoint vl <#> vr
+> a == vr <#> diagonal d <#> diagonal scal <#> adjoint vl
+
+If @a@ is non-diagonalizable then some columns of @vr@ and @vl@ are left zero
+and the above property does not hold.
+-}
+decompose ::
+   (Shape.C sh, Class.Floating a) =>
+   Square sh a ->
+   (Square sh (ComplexOf a),
+    Vector sh (ComplexOf a),
+    Square sh (ComplexOf a))
+decompose =
+   getDecompose $
+   Class.switchFloating
+      (Decompose decomposeReal)
+      (Decompose decomposeReal)
+      (Decompose decomposeComplex)
+      (Decompose decomposeComplex)
+
+newtype Decompose sh a =
+   Decompose {
+      getDecompose ::
+         Square sh a ->
+         (Square sh (ComplexOf a),
+          Vector sh (ComplexOf a),
+          Square sh (ComplexOf a))
+   }
+
+decomposeReal ::
+   (Shape.C sh, Class.Real a) =>
+   Square sh a ->
+   (Square sh (Complex a), Vector sh (Complex a), Square sh (Complex a))
+decomposeReal (Array (MatrixShape.Square order size) a) =
+      unsafePerformIO $ do
+   let n = Shape.size size
+   let lda = n
+   evalContT $ do
+      jobvlPtr <- Call.char 'V'
+      jobvrPtr <- Call.char 'V'
+      nPtr <- Call.cint n
+      aPtr <- copyToTemp (n*n) a
+      ldaPtr <- Call.cint lda
+      wrPtr <- Call.allocaArray n
+      wiPtr <- Call.allocaArray n
+      vlPtr <- Call.allocaArray (n*n)
+      ldvlPtr <- Call.cint n
+      vrPtr <- Call.allocaArray (n*n)
+      ldvrPtr <- Call.cint n
+      liftIO $ withAutoWorkspaceInfo "geev" $
+         LapackReal.geev
+            jobvlPtr jobvrPtr nPtr aPtr ldaPtr
+            wrPtr wiPtr vlPtr ldvlPtr vrPtr ldvrPtr
+      (w,wPtr) <- allocArray size
+      liftIO $ zipComplex n wrPtr wiPtr wPtr
+      let sh = MatrixShape.Square ColumnMajor size
+      (vlc,vlcPtr) <- allocArray sh
+      (vrc,vrcPtr) <- allocArray sh
+      liftIO $ eigenvectorsToComplex n wiPtr vlPtr vlcPtr
+      liftIO $ eigenvectorsToComplex n wiPtr vrPtr vrcPtr
+      return $
+         case order of
+            RowMajor -> (vlc, w, vrc)
+            ColumnMajor -> (vrc, w, vlc)
+
+eigenvectorsToComplex ::
+   (Eq a, Class.Real a) =>
+   Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
+eigenvectorsToComplex n wiPtr vPtr vcPtr = evalContT $ do
+   nPtr <- Call.cint n
+   zeroPtr <- Call.real zero
+   inc0Ptr <- Call.cint 0
+   inc1Ptr <- Call.cint 1
+   inc2Ptr <- Call.cint 2
+   liftIO $ do
+      let go _ _ [] = return ()
+          go xPtr yPtr (False:wi) = do
+            let yrPtr = castPtr yPtr
+            let yiPtr = advancePtr yrPtr 1
+            BlasReal.copy nPtr xPtr    inc1Ptr yrPtr inc2Ptr
+            BlasReal.copy nPtr zeroPtr inc0Ptr yiPtr inc2Ptr
+            go (advancePtr xPtr n) (advancePtr yPtr n) wi
+          go xPtr yPtr (True:True:wi) = do
+            let xrPtr = xPtr
+            let xiPtr = advancePtr xPtr n
+            let yrPtr = castPtr yPtr
+            let yiPtr = advancePtr yrPtr 1
+            let y1Ptr = advancePtr yPtr n
+            BlasReal.copy nPtr xrPtr inc1Ptr yrPtr inc2Ptr
+            BlasReal.copy nPtr xiPtr inc1Ptr yiPtr inc2Ptr
+            BlasComplex.copy nPtr yPtr inc1Ptr y1Ptr inc1Ptr
+            LapackComplex.lacgv nPtr y1Ptr inc1Ptr
+            go (advancePtr xPtr (2*n)) (advancePtr yPtr (2*n)) wi
+          go _xPtr _yPtr wi =
+            error $ "eigenvectorToComplex: invalid non-real pattern " ++ show wi
+      go vPtr vcPtr . map (zero/=) =<< peekArray n wiPtr
+
+decomposeComplex ::
+   (Shape.C sh, Class.Real a) =>
+   Square sh (Complex a) ->
+   (Square sh (Complex a), Vector sh (Complex a), Square sh (Complex a))
+decomposeComplex (Array (MatrixShape.Square order size) a) =
+      unsafePerformIO $ do
+   let n = Shape.size size
+   let lda = n
+   evalContT $ do
+      jobvlPtr <- Call.char 'V'
+      jobvrPtr <- Call.char 'V'
+      nPtr <- Call.cint n
+      aPtr <- copyToTemp (n*n) a
+      ldaPtr <- Call.cint lda
+      (w,wPtr) <- allocArray size
+      let sh = MatrixShape.Square ColumnMajor size
+      (vl,vlPtr) <- allocArray sh
+      ldvlPtr <- Call.cint n
+      (vr,vrPtr) <- allocArray sh
+      ldvrPtr <- Call.cint n
+      rworkPtr <- Call.allocaArray (2*n)
+      liftIO $ withAutoWorkspaceInfo "geev" $ \workPtr lworkPtr infoPtr ->
+         LapackComplex.geev
+            jobvlPtr jobvrPtr nPtr aPtr ldaPtr
+            wPtr vlPtr ldvlPtr vrPtr ldvrPtr
+            workPtr lworkPtr rworkPtr infoPtr
+      return $
+         case order of
+            RowMajor -> (vl, w, vr)
+            ColumnMajor -> (vr, w, vl)
+
+
+zipComplex ::
+   (Class.Real a) => Int -> Ptr a -> Ptr a -> Ptr (Complex a) -> IO ()
+zipComplex n vr vi vc =
+   evalContT $ do
+      nPtr <- Call.cint n
+      incxPtr <- Call.cint 1
+      incyPtr <- Call.cint 2
+      let yPtr = castPtr vc
+      liftIO $ BlasReal.copy nPtr vr incxPtr yPtr incyPtr
+      liftIO $ BlasReal.copy nPtr vi incxPtr (advancePtr yPtr 1) incyPtr
diff --git a/src/Numeric/LAPACK/Eigen/Hermitian.hs b/src/Numeric/LAPACK/Eigen/Hermitian.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Eigen/Hermitian.hs
@@ -0,0 +1,142 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Eigen.Hermitian (
+   values,
+   decompose,
+   ) where
+
+import Numeric.LAPACK.Matrix.Hermitian (Hermitian)
+import Numeric.LAPACK.Matrix.Square (Square)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(ColumnMajor), uploFromOrder, triangleSize)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private (RealOf, copyToTemp, allocArray)
+
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
+import qualified Numeric.LAPACK.FFI.Real as LapackReal
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt, CChar)
+import Foreign.Ptr (Ptr, nullPtr)
+import Foreign.Storable (Storable, peek)
+
+import Control.Monad.Trans.Cont (evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+import Data.Complex (Complex)
+
+
+values ::
+   (Shape.C sh, Class.Floating a) =>
+   Hermitian sh a -> Vector sh (RealOf a)
+values =
+   getValues $
+   Class.switchFloating
+      (Values valuesAux) (Values valuesAux)
+      (Values valuesAux) (Values valuesAux)
+
+newtype Values sh a =
+   Values {getValues :: Hermitian sh a -> Vector sh (RealOf a)}
+
+valuesAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Hermitian sh a -> Vector sh ar
+valuesAux (Array (MatrixShape.Hermitian order size) a) =
+      Array.unsafeCreateWithSize size $ \n wPtr -> do
+   evalContT $ do
+      jobzPtr <- Call.char 'N'
+      uploPtr <- Call.char $ uploFromOrder order
+      aPtr <- copyToTemp (triangleSize n) a
+      let zPtr = nullPtr
+      ldzPtr <- Call.cint (max 1 n)
+      liftIO $ withInfo "hpev" $
+         hpev jobzPtr uploPtr n aPtr wPtr zPtr ldzPtr
+
+
+{- |
+For symmetric eigenvalue problems, @Eigen.decompose@ and @schur@ coincide.
+-}
+decompose ::
+   (Shape.C sh, Class.Floating a) =>
+   Hermitian sh a -> (Square sh a, Vector sh (RealOf a))
+decompose =
+   getDecompose $
+   Class.switchFloating
+      (Decompose decomposeAux) (Decompose decomposeAux)
+      (Decompose decomposeAux) (Decompose decomposeAux)
+
+type Decompose_ sh a = Hermitian sh a -> (Square sh a, Vector sh (RealOf a))
+
+newtype Decompose sh a = Decompose {getDecompose :: Decompose_ sh a}
+
+decomposeAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Decompose_ sh a
+decomposeAux (Array (MatrixShape.Hermitian order size) a) = unsafePerformIO $ do
+   let n = Shape.size size
+   let shZ = MatrixShape.Square ColumnMajor size
+   evalContT $ do
+      jobzPtr <- Call.char 'V'
+      uploPtr <- Call.char $ uploFromOrder order
+      aPtr <- copyToTemp (triangleSize n) a
+      (w,wPtr) <- allocArray size
+      (z,zPtr) <- allocArray shZ
+      ldzPtr <- Call.cint n
+      liftIO $ withInfo "hpev" $
+         hpev jobzPtr uploPtr n aPtr wPtr zPtr ldzPtr
+      return (z, w)
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $
+         printf "%s: %d off-diagonal elements not converging" name info
+
+
+type HPEV_ ar a =
+   Ptr CChar -> Ptr CChar -> Int -> Ptr a -> Ptr ar ->
+   Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
+
+newtype HPEV a = HPEV {getHPEV :: HPEV_ (RealOf a) a}
+
+hpev :: Class.Floating a => HPEV_ (RealOf a) a
+hpev =
+   getHPEV $
+   Class.switchFloating
+      (HPEV spevReal) (HPEV spevReal) (HPEV hpevComplex) (HPEV hpevComplex)
+
+spevReal :: Class.Real a => HPEV_ a a
+spevReal jobzPtr uploPtr n apPtr wPtr zPtr ldzPtr infoPtr =
+   evalContT $ do
+      nPtr <- Call.cint n
+      workPtr <- Call.allocaArray (3*n)
+      liftIO $
+         LapackReal.spev
+            jobzPtr uploPtr nPtr apPtr wPtr zPtr ldzPtr workPtr infoPtr
+
+hpevComplex :: Class.Real a => HPEV_ a (Complex a)
+hpevComplex jobzPtr uploPtr n apPtr wPtr zPtr ldzPtr infoPtr =
+   evalContT $ do
+      nPtr <- Call.cint n
+      workPtr <- Call.allocaArray (max 1 (2*n-1))
+      rworkPtr <- Call.allocaArray (max 1 (3*n-2))
+      liftIO $
+         LapackComplex.hpev
+            jobzPtr uploPtr nPtr apPtr wPtr zPtr ldzPtr workPtr rworkPtr infoPtr
diff --git a/src/Numeric/LAPACK/Eigen/Triangular.hs b/src/Numeric/LAPACK/Eigen/Triangular.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Eigen/Triangular.hs
@@ -0,0 +1,176 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Eigen.Triangular (
+   values,
+   decompose,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Triangular.Private
+         (unpackZero, pack, unpackToTemp, fillTriangle,
+          forPointers, rowMajorPointers)
+import Numeric.LAPACK.Matrix.Triangular (Triangular)
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(ColumnMajor,RowMajor), caseUplo, uploOrder, triangleSize)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private (zero, lacgv, allocArray)
+
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
+import qualified Numeric.LAPACK.FFI.Real as LapackReal
+import qualified Numeric.BLAS.FFI.Generic as BlasGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt, CChar)
+import Foreign.Ptr (Ptr, nullPtr)
+import Foreign.Storable (peek)
+
+import Control.Monad.Trans.Cont (evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+import Data.Complex (Complex)
+
+
+values ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Triangular uplo sh a -> Vector sh a
+values = Triangular.getDiagonal
+
+
+{- |
+@(vr,d,vlAdj) = TriEigen.decompose a@
+
+Counterintuitively, @vr@ contains the right eigenvectors as columns
+and @vlAdj@ contains the left conjugated eigenvectors as rows.
+The idea is to provide a decomposition of @a@.
+If @a@ is diagonalizable, then @vr@ and @vlAdj@
+are almost inverse to each other.
+More precisely, @vlAdj \<#\> vr@ is a diagonal matrix.
+This is because the eigenvectors are not normalized.
+With the following scaling, the decomposition becomes perfect:
+
+> let scal = Array.map recip $ getDiagonal $ vlAdj <#> vr
+> a == vr <#> diagonal d <#> diagonal scal <#> vlAdj
+
+If @a@ is non-diagonalizable
+then some columns of @vr@ and corresponding rows of @vlAdj@ are left zero
+and the above property does not hold.
+-}
+decompose ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Triangular uplo sh a ->
+   (Triangular uplo sh a, Vector sh a, Triangular uplo sh a)
+decompose a =
+   let (vr,vl) = decomposePlain a
+   in  (vr, values a, vl)
+
+decomposePlain ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Triangular uplo sh a -> (Triangular uplo sh a, Triangular uplo sh a)
+decomposePlain (Array (MatrixShape.Triangular uplo order sh) a) =
+      unsafePerformIO $ do
+   let n = Shape.size sh
+   let n2 = n*n
+   let triSize = triangleSize n
+   evalContT $ do
+      sidePtr <- Call.char 'B'
+      howManyPtr <- Call.char 'A'
+      let selectPtr = nullPtr
+      let unpk =
+            case uploOrder uplo order of
+               ColumnMajor -> unpackZero ColumnMajor
+               RowMajor -> unpackZeroRowMajor
+      aPtr <- unpackToTemp unpk n a
+      ldaPtr <- Call.cint n
+      vlPtr <- Call.allocaArray n2
+      vrPtr <- Call.allocaArray n2
+      mmPtr <- Call.cint n
+      mPtr <- Call.alloca
+      liftIO $ withInfo "trevc" $
+         trevc sidePtr howManyPtr selectPtr n
+            aPtr ldaPtr vlPtr ldaPtr vrPtr ldaPtr mmPtr mPtr
+      (vl,vlpPtr) <-
+         allocArray $
+         MatrixShape.Triangular uplo (uploOrder uplo RowMajor) sh
+      (vr,vrpPtr) <-
+         allocArray $
+         MatrixShape.Triangular uplo (uploOrder uplo ColumnMajor) sh
+      sizePtr <- Call.cint triSize
+      incPtr <- Call.cint 1
+      liftIO $ do
+         pack ColumnMajor n vrPtr vrpPtr
+         pack RowMajor n vlPtr vlpPtr
+         lacgv sizePtr vlpPtr incPtr
+      return $ caseUplo uplo (vl,vr) (vr,vl)
+
+
+unpackZeroRowMajor :: Class.Floating a => Int -> Ptr a -> Ptr a -> IO ()
+unpackZeroRowMajor n packedPtr fullPtr = do
+   fillTriangle zero RowMajor n fullPtr
+   unpackRowMajor n packedPtr fullPtr
+
+unpackRowMajor :: Class.Floating a => Int -> Ptr a -> Ptr a -> IO ()
+unpackRowMajor n packedPtr fullPtr = evalContT $ do
+   incxPtr <- Call.cint 1
+   incyPtr <- Call.cint n
+   liftIO $
+      forPointers (rowMajorPointers n fullPtr packedPtr) $
+            \nPtr (dstPtr,srcPtr) ->
+         BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $
+         printf "%s: %d off-diagonal elements not converging" name info
+
+
+type TREVC_ a =
+   Ptr CChar -> Ptr CChar -> Ptr Bool ->
+   Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
+   Ptr CInt -> Ptr CInt -> Ptr CInt -> IO ()
+
+newtype TREVC a = TREVC {getTREVC :: TREVC_ a}
+
+trevc :: Class.Floating a => TREVC_ a
+trevc =
+   getTREVC $
+   Class.switchFloating
+      (TREVC trevcReal) (TREVC trevcReal)
+      (TREVC trevcComplex) (TREVC trevcComplex)
+
+trevcReal :: Class.Real a => TREVC_ a
+trevcReal sidePtr howmnyPtr selectPtr n
+      tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr infoPtr =
+   evalContT $ do
+      nPtr <- Call.cint n
+      workPtr <- Call.allocaArray (3*n)
+      liftIO $
+         LapackReal.trevc sidePtr howmnyPtr selectPtr nPtr
+            tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr workPtr infoPtr
+
+trevcComplex :: Class.Real a => TREVC_ (Complex a)
+trevcComplex sidePtr howmnyPtr selectPtr n
+      tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr infoPtr =
+   evalContT $ do
+      nPtr <- Call.cint n
+      workPtr <- Call.allocaArray (2*n)
+      rworkPtr <- Call.allocaArray n
+      liftIO $
+         LapackComplex.trevc sidePtr howmnyPtr selectPtr nPtr
+            tPtr ldtPtr vlPtr ldvlPtr vrPtr ldvrPtr mmPtr mPtr
+            workPtr rworkPtr infoPtr
diff --git a/src/Numeric/LAPACK/Format.hs b/src/Numeric/LAPACK/Format.hs
--- a/src/Numeric/LAPACK/Format.hs
+++ b/src/Numeric/LAPACK/Format.hs
@@ -1,7 +1,13 @@
-module Numeric.LAPACK.Format where
+module Numeric.LAPACK.Format (
+   (##),
+   Format(format),
+   FormatArray(formatArray),
+   ) where
 
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Square as Square
 import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))
+import Numeric.LAPACK.Matrix.Private (General)
 
 import qualified Numeric.Netlib.Class as Class
 
@@ -9,11 +15,12 @@
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable (Array)
 
-import Foreign.Storable (Storable)
-
 import Text.Printf (PrintfArg, printf)
 
 import qualified Data.List.HT as ListHT
+import qualified Data.Complex as Complex
+import Data.Monoid (Endo(Endo,appEndo))
+import Data.List (mapAccumL, transpose)
 import Data.Complex (Complex((:+)))
 
 
@@ -35,7 +42,7 @@
 instance Format Double where
    format fmt a = [printf fmt a]
 
-instance (PrintfArg a) => Format (Complex a) where
+instance (PrintfArg a, Class.Real a) => Format (Complex a) where
    format fmt a = [printfComplex fmt a]
 
 instance (Format a, Format b) => Format (a,b) where
@@ -45,15 +52,12 @@
    format fmt (a,b,c) =
       format fmt a ++ [""] ++ format fmt b ++ [""] ++ format fmt c
 
-instance
-   (FormatArray sh, Class.Floating a, Storable a) =>
-      Format (Array sh a) where
+instance (FormatArray sh, Class.Floating a) => Format (Array sh a) where
    format = formatArray
 
 
 class (Shape.C sh) => FormatArray sh where
-   formatArray ::
-      (Storable a, Class.Floating a) => String -> Array sh a -> [String]
+   formatArray :: (Class.Floating a) => String -> Array sh a -> [String]
 
 instance (Integral i) => FormatArray (Shape.ZeroBased i) where
    formatArray fmt m = [unwords $ map (printfFloating fmt) $ Array.toList m]
@@ -61,19 +65,20 @@
 instance (Integral i) => FormatArray (Shape.OneBased i) where
    formatArray fmt m = [unwords $ map (printfFloating fmt) $ Array.toList m]
 
+instance (Shape.C sh) => FormatArray (MatrixShape.Square sh) where
+   formatArray fmt = formatGeneral fmt . Square.toGeneral
+
 instance
    (Shape.C height, Shape.C width) =>
       FormatArray (MatrixShape.General height width) where
    formatArray = formatGeneral
 
 formatGeneral ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
-   String -> Array (MatrixShape.General height width) a -> [String]
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   String -> General height width a -> [String]
 formatGeneral fmt m =
    let MatrixShape.General order height width = Array.shape m
-       xss = formatRows fmt order (height,width) $ Array.toList m
-       strWidths = columnWidths xss
-   in  map (unwords . zipWith (ListHT.padLeft ' ') strWidths) xss
+   in  formatAligned $ formatRows fmt order (height,width) $ Array.toList m
 
 instance
    (Shape.C height, Shape.C width) =>
@@ -81,12 +86,94 @@
    formatArray = formatHouseholder
 
 formatHouseholder ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    String -> Array (MatrixShape.Householder height width) a -> [String]
 formatHouseholder fmt m =
    let MatrixShape.Householder order height width = Array.shape m
-       xss = formatRows fmt order (height,width) $ Array.toList m
-       strWidths = columnWidths xss
+   in formatSeparateTriangle $
+      formatRows fmt order (height,width) $ Array.toList m
+
+instance (Shape.C size) => FormatArray (MatrixShape.Hermitian size) where
+   formatArray = formatHermitian
+
+formatHermitian ::
+   (Shape.C size, Class.Floating a) =>
+   String -> Array (MatrixShape.Hermitian size) a -> [String]
+formatHermitian fmt m =
+   let MatrixShape.Hermitian order size = Array.shape m
+   in  formatSeparateTriangle $
+       map (map (printfFloating fmt)) $
+       complementTriangle order (Shape.size size) $ Array.toList m
+
+complementTriangle :: (Class.Floating a) => Order -> Int -> [a] -> [[a]]
+complementTriangle order n xs =
+   let mergeTriangles lower upper =
+         zipWith (++) (map (map conjugate . init) lower) upper
+   in case order of
+         RowMajor ->
+            let tri = slice (take n $ iterate pred n) xs
+                trans = reverse $ transpose $ map reverse tri
+            in  mergeTriangles trans tri
+         ColumnMajor ->
+            let tri = slice (take n [1..]) xs
+            in  mergeTriangles tri (transpose tri)
+
+conjugate :: (Class.Floating a) => a -> a
+conjugate =
+   appEndo $
+   Class.switchFloating
+      (Endo id) (Endo id) (Endo Complex.conjugate) (Endo Complex.conjugate)
+
+instance
+   (MatrixShape.Uplo uplo, Shape.C size) =>
+      FormatArray (MatrixShape.Triangular uplo size) where
+   formatArray = formatTriangular
+
+formatTriangular ::
+   (MatrixShape.Uplo uplo, Shape.C size, Class.Floating a) =>
+   String -> Array (MatrixShape.Triangular uplo size) a -> [String]
+formatTriangular fmt m =
+   let MatrixShape.Triangular uplo order size = Array.shape m
+   in  formatAligned $
+       MatrixShape.caseUplo uplo
+         padLowerTriangle padUpperTriangle order (Shape.size size) $
+       map (printfFloating fmt) $ Array.toList m
+
+padUpperTriangle :: Order -> Int -> [String] -> [[String]]
+padUpperTriangle order n xs =
+   case order of
+      RowMajor ->
+         zipWith (++) (iterate ("":) []) (slice (take n $ iterate pred n) xs)
+      ColumnMajor ->
+         transpose $
+         zipWith (++)
+            (slice (take n [1..]) xs)
+            (reverse $ take n $ iterate ("":) [])
+
+padLowerTriangle :: Order -> Int -> [String] -> [[String]]
+padLowerTriangle order n xs =
+   case order of
+      RowMajor ->
+         map (take n) $ map (++ repeat "") $ slice (take n [1..]) xs
+      ColumnMajor ->
+         transpose $
+         zipWith (++) (iterate ("":) []) (slice (take n $ iterate pred n) xs)
+
+_padLowerTriangle :: Order -> Int -> [a] -> [[a]]
+_padLowerTriangle order n xs =
+   case order of
+      RowMajor -> slice (take n [1..]) xs
+      ColumnMajor ->
+         foldr (\(y:ys) zs -> [y] : zipWith (:) ys zs) []
+            (slice (take n $ iterate pred n) xs)
+
+slice :: [Int] -> [a] -> [[a]]
+slice ns xs =
+   snd $ mapAccumL (\ys n -> let (vs,ws) = splitAt n ys in (ws,vs)) xs ns
+
+formatSeparateTriangle :: [[String]] -> [String]
+formatSeparateTriangle xss =
+   let strWidths = columnWidths xss
    in  zipWith
          (\row xs ->
             concat $
@@ -103,6 +190,11 @@
       ColumnMajor -> ListHT.sliceHorizontal (Shape.size height)) .
    map (printfFloating fmt)
 
+formatAligned :: [[String]] -> [String]
+formatAligned xss =
+   let strWidths = columnWidths xss
+   in  map (unwords . zipWith (ListHT.padLeft ' ') strWidths) xss
+
 columnWidths :: [[[a]]] -> [Int]
 columnWidths xss =
    case map (map length) xss of
@@ -121,5 +213,8 @@
       (Printf printfComplex)
       (Printf printfComplex)
 
-printfComplex :: (PrintfArg a) => String -> Complex a -> String
-printfComplex fmt (r:+i) = printf (fmt ++ "+i" ++ fmt) r i
+printfComplex :: (PrintfArg a, Class.Real a) => String -> Complex a -> String
+printfComplex fmt (r:+i) =
+   if i<0 || isNegativeZero i
+     then printf (fmt ++ "-i" ++ fmt) r (-i)
+     else printf (fmt ++ "+i" ++ fmt) r i
diff --git a/src/Numeric/LAPACK/Linear/General.hs b/src/Numeric/LAPACK/Linear/General.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Linear/General.hs
@@ -0,0 +1,88 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Linear.General (
+   solve,
+   inverse,
+   ) where
+
+import Numeric.LAPACK.Matrix.Square (Square)
+import Numeric.LAPACK.Matrix (General)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor))
+import Numeric.LAPACK.Private (withAutoWorkspace, copyBlock, copyToColumnMajor)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+
+solve ::
+   (Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
+   Square sh a -> General sh nrhs a -> General sh nrhs a
+solve
+   (Array (MatrixShape.Square orderA shA) a)
+   (Array (MatrixShape.General orderB heightB widthB) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $
+         \xPtr -> do
+   Call.assert "Square.solve: height shapes mismatch"
+      (shA == heightB)
+   let n = Shape.size heightB
+   let nrhs = Shape.size widthB
+   let ldb = n
+   evalContT $ do
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      aPtr <- ContT $ withForeignPtr a
+      atmpPtr <- Call.allocaArray (n*n)
+      ldaPtr <- Call.cint ldb
+      ipivPtr <- Call.allocaArray n
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint ldb
+      liftIO $ do
+         copyToColumnMajor orderA n n aPtr atmpPtr
+         copyToColumnMajor orderB n nrhs bPtr xPtr
+         withInfo "gesv" $
+            LapackGen.gesv nPtr nrhsPtr atmpPtr ldaPtr ipivPtr xPtr ldbPtr
+
+
+inverse :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
+inverse (Array shape@(MatrixShape.Square _order sh) a) =
+      Array.unsafeCreateWithSize shape $ \blockSize bPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      nPtr <- Call.cint n
+      aPtr <- ContT $ withForeignPtr a
+      ldbPtr <- Call.cint n
+      ipivPtr <- Call.allocaArray n
+      liftIO $ do
+         copyBlock blockSize aPtr bPtr
+         withInfo "getrf" $ LapackGen.getrf nPtr nPtr bPtr ldbPtr ipivPtr
+         withInfo "getri" $ \infoPtr ->
+            withAutoWorkspace $ \workPtr lworkPtr ->
+               LapackGen.getri nPtr bPtr ldbPtr ipivPtr workPtr lworkPtr infoPtr
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $ printf "%s: %d-th diagonal value is zero" name info
diff --git a/src/Numeric/LAPACK/Linear/Hermitian.hs b/src/Numeric/LAPACK/Linear/Hermitian.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Linear/Hermitian.hs
@@ -0,0 +1,87 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Linear.Hermitian (
+   solve,
+   inverse,
+   ) where
+
+import Numeric.LAPACK.Matrix.Hermitian (Hermitian)
+import Numeric.LAPACK.Matrix (General)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp)
+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)
+import Numeric.LAPACK.Private (copyBlock, copyToColumnMajor)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+
+solve ::
+   (Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
+   Hermitian sh a -> General sh nrhs a -> General sh nrhs a
+solve
+   (Array (MatrixShape.Hermitian orderA shA) a)
+   (Array (MatrixShape.General orderB heightB widthB) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $
+         \xPtr -> do
+   Call.assert "Hermitian.solve: height shapes mismatch"
+      (shA == heightB)
+   let n = Shape.size heightB
+   let nrhs = Shape.size widthB
+   let ldb = n
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder orderA
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      apPtr <- copyTriangleToTemp orderA n a
+      ipivPtr <- Call.allocaArray n
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint ldb
+      liftIO $ do
+         copyToColumnMajor orderB n nrhs bPtr xPtr
+         withInfo "hpsv" $
+            LapackGen.hpsv uploPtr nPtr nrhsPtr apPtr ipivPtr xPtr ldbPtr
+
+
+inverse ::
+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a
+inverse (Array shape@(MatrixShape.Hermitian order sh) a) =
+      Array.unsafeCreateWithSize shape $ \triSize bPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder order
+      nPtr <- Call.cint n
+      aPtr <- ContT $ withForeignPtr a
+      ipivPtr <- Call.allocaArray n
+      workPtr <- Call.allocaArray n
+      liftIO $ do
+         copyBlock triSize aPtr bPtr
+         withInfo "hptrf" $ LapackGen.hptrf uploPtr nPtr bPtr ipivPtr
+         withInfo "hptri" $ LapackGen.hptri uploPtr nPtr bPtr ipivPtr workPtr
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $ printf "%s: %d-th diagonal value is zero" name info
diff --git a/src/Numeric/LAPACK/Linear/HermitianPositiveDefinite.hs b/src/Numeric/LAPACK/Linear/HermitianPositiveDefinite.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Linear/HermitianPositiveDefinite.hs
@@ -0,0 +1,105 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Linear.HermitianPositiveDefinite (
+   solve,
+   inverse,
+   decompose,
+   ) where
+
+import Numeric.LAPACK.Matrix.Hermitian (Hermitian)
+import Numeric.LAPACK.Matrix.Triangular (Upper)
+import Numeric.LAPACK.Matrix (General)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Triangular.Private (copyTriangleToTemp)
+import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)
+import Numeric.LAPACK.Private (copyBlock, copyToColumnMajor)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+
+solve ::
+   (Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
+   Hermitian sh a -> General sh nrhs a -> General sh nrhs a
+solve
+   (Array (MatrixShape.Hermitian orderA shA) a)
+   (Array (MatrixShape.General orderB heightB widthB) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $
+         \xPtr -> do
+   Call.assert "Hermitian.solve: height shapes mismatch"
+      (shA == heightB)
+   let n = Shape.size heightB
+   let nrhs = Shape.size widthB
+   let ldb = n
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder orderA
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      apPtr <- copyTriangleToTemp orderA n a
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint ldb
+      liftIO $ do
+         copyToColumnMajor orderB n nrhs bPtr xPtr
+         withInfo "ppsv" $
+            LapackGen.ppsv uploPtr nPtr nrhsPtr apPtr xPtr ldbPtr
+
+
+inverse ::
+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a
+inverse
+   (Array shape@(MatrixShape.Hermitian order sh) a) =
+      Array.unsafeCreateWithSize shape $ \triSize bPtr -> do
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder order
+      nPtr <- Call.cint $ Shape.size sh
+      aPtr <- ContT $ withForeignPtr a
+      liftIO $ do
+         copyBlock triSize aPtr bPtr
+         withInfo "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr
+         withInfo "pptri" $ LapackGen.pptri uploPtr nPtr bPtr
+
+{- |
+Cholesky decomposition
+-}
+decompose ::
+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Upper sh a
+decompose
+   (Array (MatrixShape.Hermitian order sh) a) =
+      Array.unsafeCreateWithSize
+         (MatrixShape.Triangular MatrixShape.Upper order sh) $
+            \triSize bPtr -> do
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder order
+      nPtr <- Call.cint $ Shape.size sh
+      aPtr <- ContT $ withForeignPtr a
+      liftIO $ do
+         copyBlock triSize aPtr bPtr
+         withInfo "pptrf" $ LapackGen.pptrf uploPtr nPtr bPtr
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $
+         printf "%s: minor of order %d not positive definite" name info
diff --git a/src/Numeric/LAPACK/Linear/Triangular.hs b/src/Numeric/LAPACK/Linear/Triangular.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Linear/Triangular.hs
@@ -0,0 +1,88 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Linear.Triangular (
+   solve,
+   inverse,
+   ) where
+
+import Numeric.LAPACK.Matrix.Triangular (Triangular)
+import Numeric.LAPACK.Matrix (General)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(ColumnMajor),
+          transposeFromOrder, uploFromOrder, uploOrder, triangleSize)
+import Numeric.LAPACK.Private
+         (copyBlock, copyToTemp, copyToColumnMajor)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+
+solve ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
+   Triangular uplo sh a -> General sh nrhs a -> General sh nrhs a
+solve
+   (Array (MatrixShape.Triangular uplo orderA shA) a)
+   (Array (MatrixShape.General orderB heightB widthB) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor heightB widthB) $
+         \xPtr -> do
+   Call.assert "Triangular.solve: height shapes mismatch" (shA == heightB)
+   let n = Shape.size heightB
+   let nrhs = Shape.size widthB
+   let ldb = n
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo orderA
+      transPtr <- Call.char $ transposeFromOrder orderA
+      diagPtr <- Call.char 'N'
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      apPtr <- copyToTemp (triangleSize n) a
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint ldb
+      liftIO $ do
+         copyToColumnMajor orderB n nrhs bPtr xPtr
+         withInfo "tptrs" $
+            LapackGen.tptrs uploPtr transPtr diagPtr
+               nPtr nrhsPtr apPtr xPtr ldbPtr
+
+
+inverse ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Triangular uplo sh a -> Triangular uplo sh a
+inverse (Array shape@(MatrixShape.Triangular uplo order sh) a) =
+      Array.unsafeCreateWithSize shape $ \triSize bPtr -> do
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo order
+      diagPtr <- Call.char 'N'
+      nPtr <- Call.cint $ Shape.size sh
+      aPtr <- ContT $ withForeignPtr a
+      liftIO $ do
+         copyBlock triSize aPtr bPtr
+         withInfo "tptri" $ LapackGen.tptri uploPtr diagPtr nPtr bPtr
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $ printf "%s: %d-th diagonal element zero" name info
diff --git a/src/Numeric/LAPACK/LinearSystem.hs b/src/Numeric/LAPACK/LinearSystem.hs
deleted file mode 100644
--- a/src/Numeric/LAPACK/LinearSystem.hs
+++ /dev/null
@@ -1,478 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-module Numeric.LAPACK.LinearSystem (
-   leastSquares,
-   minimumNorm,
-   leastSquaresMinimumNorm,
-   pseudoInverseRCond,
-
-   Householder,
-   householder,
-   householderDecompose,
-   householderDeterminant,
-   determinant,
-   householderExtractQ,
-   householderExtractR,
-   orthogonalComplement,
-   ) where
-
-import Numeric.LAPACK.Matrix
-         (General, ZeroInt, zeroInt, transpose, identity, dropColumns)
-
-import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
-import Numeric.LAPACK.Matrix.Shape.Private
-         (Order(RowMajor, ColumnMajor), charFromOrder)
-import Numeric.LAPACK.Vector (Vector)
-import Numeric.LAPACK.Private
-         (RealOf, zero, one, fill, pointerSeq,
-          copyTransposed, copySubMatrix, copyBlock)
-
-import qualified Numeric.LAPACK.FFI.Generic as LapackGen
-import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
-import qualified Numeric.LAPACK.FFI.Real as LapackReal
-import qualified Numeric.Netlib.Utility as Call
-import qualified Numeric.Netlib.Class as Class
-
-import qualified Data.Array.Comfort.Storable.Internal as Array
-import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Internal (Array(Array))
-
-import System.IO.Unsafe (unsafePerformIO)
-
-import Foreign.Marshal.Array (allocaArray, advancePtr)
-import Foreign.Marshal.Alloc (alloca)
-import Foreign.C.Types (CInt)
-import Foreign.ForeignPtr (withForeignPtr, mallocForeignPtrArray)
-import Foreign.Ptr (Ptr)
-import Foreign.Storable (Storable, poke, peek)
-
-import Text.Printf (printf)
-
-import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
-import Control.Monad.IO.Class (liftIO)
-import Control.Monad (when, foldM)
-import Control.Applicative ((<$>))
-
-import qualified Data.Complex as Complex
-import Data.Complex (Complex)
-import Data.Tuple.HT (mapSnd)
-
-
-{- |
-If @x = leastSquares a b@
-then @x@ minimizes @Vector.norm2 (multiply a x `sub` b)@.
-
-Precondition: @a@ must have full rank and @height a >= width a@.
--}
-leastSquares ::
-   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,
-    Storable a, Class.Floating a) =>
-   General height width a -> General height nrhs a -> General width nrhs a
-leastSquares
-   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)
-   (Array        (MatrixShape.General orderB heightB widthB) b) =
-      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $
-         \xPtr -> do
-   Call.assert "leastSquares: height shapes mismatch" (heightA == heightB)
-   Call.assert "leastSquares: height of 'a' must be at least the width"
-      (Shape.size heightA >= Shape.size widthA)
-   let (m,n) = MatrixShape.dimensions shapeA
-   let lda = m
-   let nrhs = Shape.size widthB
-   let ldb = Shape.size heightB
-   let ldx = Shape.size widthA
-   evalContT $ do
-      transPtr <- Call.char $ charFromOrder orderA
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      nrhsPtr <- Call.cint nrhs
-      aPtr <- ContT $ withForeignPtr a
-      ldaPtr <- Call.cint lda
-      let aSize = Shape.size (heightA,widthA)
-      atmpPtr <- Call.allocaArray aSize
-      liftIO $ copyBlock aSize aPtr atmpPtr
-      bPtr <- ContT $ withForeignPtr b
-      ldbPtr <- Call.cint ldb
-      let bSize = Shape.size (heightB,widthB)
-      btmpPtr <- Call.allocaArray bSize
-      liftIO $ copyToColumnMajor orderB ldb nrhs bPtr btmpPtr
-      liftIO $ withAutoWorkspaceInfo "gels" $
-         LapackGen.gels transPtr
-            mPtr nPtr nrhsPtr atmpPtr ldaPtr btmpPtr ldbPtr
-      liftIO $ copySubMatrix ldx nrhs ldb btmpPtr ldx xPtr
-
-{- |
-The vector @x@ with @x = minimumNorm a b@
-is the vector with minimal @Vector.norm2 x@
-that satisfies @multiply a x == b@.
-
-Precondition: @a@ must have full rank and @height a <= width a@.
--}
-minimumNorm ::
-   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,
-    Storable a, Class.Floating a) =>
-   General height width a -> General height nrhs a -> General width nrhs a
-minimumNorm
-   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)
-   (Array        (MatrixShape.General orderB heightB widthB) b) =
-      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $
-         \xPtr -> do
-   Call.assert "minimumNorm: height shapes mismatch" (heightA == heightB)
-   Call.assert "minimumNorm: width of 'a' must be at least the height"
-      (Shape.size widthA >= Shape.size heightA)
-   let (m,n) = MatrixShape.dimensions shapeA
-   let lda = m
-   let nrhs = Shape.size widthB
-   let ldb = Shape.size heightB
-   let ldx = Shape.size widthA
-   evalContT $ do
-      transPtr <- Call.char $ charFromOrder orderA
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      nrhsPtr <- Call.cint nrhs
-      aPtr <- ContT $ withForeignPtr a
-      ldaPtr <- Call.cint lda
-      let aSize = Shape.size (heightA,widthA)
-      atmpPtr <- Call.allocaArray aSize
-      liftIO $ copyBlock aSize aPtr atmpPtr
-      bPtr <- ContT $ withForeignPtr b
-      ldxPtr <- Call.cint ldx
-      liftIO $ copyToSubColumnMajor orderB ldb nrhs bPtr ldx xPtr
-      liftIO $ withAutoWorkspaceInfo "gels" $
-         LapackGen.gels transPtr
-            mPtr nPtr nrhsPtr atmpPtr ldaPtr xPtr ldxPtr
-
-{- |
-If @x = leastSquaresMinimumNorm a b@
-then @x@ is the vector with minimum @Vector.norm2 x@
-that minimizes @Vector.norm2 (multiply a x `sub` b)@.
-
-Matrix @a@ can have any rank
-but you must specify the reciprocal condition of the rank-truncated matrix.
--}
-leastSquaresMinimumNorm ::
-   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,
-    Storable a, Class.Floating a) =>
-   RealOf a ->
-   General height width a -> General height nrhs a ->
-   (Int, General width nrhs a)
-leastSquaresMinimumNorm rcond
-   (Array (MatrixShape.General orderA heightA widthA) a)
-   (Array (MatrixShape.General orderB heightB widthB) b) =
-      unsafePerformIO $ do
-   Call.assert "minimumNorm: height shapes mismatch" (heightA == heightB)
-   let shapeX = MatrixShape.General ColumnMajor widthA widthB
-   let m = Shape.size heightA
-   let n = Shape.size widthA
-   let nrhs = Shape.size widthB
-   let aSize = m*n
-   let lda = m
-   let ldtmp = max m n
-   let tmpSize = ldtmp*nrhs
-   evalContT $ do
-      aPtr <- ContT $ withForeignPtr a
-      atmpPtr <- Call.allocaArray aSize
-      liftIO $ copyToColumnMajor orderA m n aPtr atmpPtr
-      ldaPtr <- Call.cint lda
-      bPtr <- ContT $ withForeignPtr b
-      let needTmp = m>n
-      x <- liftIO $ mallocForeignPtrArray $ Shape.size shapeX
-      tmpPtr <-
-         ContT $ if needTmp then allocaArray tmpSize else withForeignPtr x
-      ldtmpPtr <- Call.cint ldtmp
-      liftIO $ copyToSubColumnMajor orderB m nrhs bPtr ldtmp tmpPtr
-      jpvtPtr <- Call.allocaArray n
-      rankPtr <- Call.alloca
-      gelsy m n nrhs atmpPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr
-      when needTmp $ liftIO $
-         withForeignPtr x $ copySubMatrix n nrhs ldtmp tmpPtr n
-      rank <- liftIO $ fromIntegral <$> peek rankPtr
-      return (rank, Array shapeX x)
-
-
-newtype GELSY r a =
-   GELSY {
-      getGELSY ::
-         Int -> Int -> Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
-         Ptr CInt -> RealOf a -> Ptr CInt -> ContT r IO ()
-   }
-
-gelsy ::
-   (Class.Floating a) =>
-   Int -> Int -> Int ->
-   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
-   Ptr CInt -> RealOf a -> Ptr CInt -> ContT r IO ()
-gelsy =
-   getGELSY $
-   Class.switchFloating
-      (GELSY gelsyReal)
-      (GELSY gelsyReal)
-      (GELSY gelsyComplex)
-      (GELSY gelsyComplex)
-
-gelsyReal ::
-   (Class.Real a, Class.Floating a) =>
-   Int -> Int -> Int ->
-   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
-   Ptr CInt -> a -> Ptr CInt -> ContT r IO ()
-gelsyReal m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do
-   mPtr <- Call.cint m
-   nPtr <- Call.cint n
-   nrhsPtr <- Call.cint nrhs
-   rcondPtr <- Call.real rcond
-   liftIO $ withAutoWorkspaceInfo "gelsy" $
-      LapackReal.gelsy mPtr nPtr nrhsPtr
-         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr
-
-gelsyComplex ::
-   (Class.Real a) =>
-   Int -> Int -> Int ->
-   Ptr (Complex a) -> Ptr CInt -> Ptr (Complex a) -> Ptr CInt ->
-   Ptr CInt -> a -> Ptr CInt -> ContT r IO ()
-gelsyComplex m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do
-   mPtr <- Call.cint m
-   nPtr <- Call.cint n
-   nrhsPtr <- Call.cint nrhs
-   rcondPtr <- Call.real rcond
-   rworkPtr <- Call.allocaArray (2*n)
-   liftIO $ withAutoWorkspaceInfo "gelsy" $ \workPtr lworkPtr infoPtr ->
-      LapackComplex.gelsy mPtr nPtr nrhsPtr
-         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr
-         workPtr lworkPtr rworkPtr infoPtr
-
-
-pseudoInverseRCond ::
-   (Shape.C height, Eq height, Shape.C width, Eq width,
-    Storable a, Class.Floating a) =>
-   RealOf a -> General height width a -> (Int, General width height a)
-pseudoInverseRCond rcond a =
-   let (MatrixShape.General _ height width) = Array.shape a
-   in if Shape.size height < Shape.size width
-         then leastSquaresMinimumNorm rcond a $ identity height
-         else mapSnd transpose $
-              leastSquaresMinimumNorm rcond (transpose a) $
-              identity width
-
-
-type Householder height width = Array (MatrixShape.Householder height width)
-
-{-
-@(q,r) = householder a@
-means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.
--}
-householder ::
-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>
-   General height width a ->
-   (General height height a, General height width a)
-householder a =
-   let hh = householderDecompose a
-   in  (householderExtractQ hh, householderExtractR $ snd hh)
-
-householderDecompose ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
-   General height width a -> (Vector width a, Householder height width a)
-householderDecompose (Array (MatrixShape.General order height width) a) =
-   unsafePerformIO $ do
-
-   let (m,n) =
-         case order of
-            RowMajor -> (Shape.size width, Shape.size height)
-            ColumnMajor -> (Shape.size height, Shape.size width)
-   let lda = m
-   let mn = min m n
-   evalContT $ do
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      aPtr <- ContT $ withForeignPtr a
-      ldaPtr <- Call.cint lda
-      qr <- liftIO $ mallocForeignPtrArray (m*n)
-      qrPtr <- ContT $ withForeignPtr qr
-      liftIO $ copyBlock (m*n) aPtr qrPtr
-      tau <- liftIO $ mallocForeignPtrArray n
-      tauPtr <- ContT $ withForeignPtr tau
-      liftIO $ fill zero (n-mn) (advancePtr tauPtr mn)
-      liftIO $
-         case order of
-            RowMajor ->
-               withAutoWorkspaceInfo "gelqf" $
-                  LapackGen.gelqf mPtr nPtr qrPtr ldaPtr tauPtr
-            ColumnMajor ->
-               withAutoWorkspaceInfo "geqrf" $
-                  LapackGen.geqrf mPtr nPtr qrPtr ldaPtr tauPtr
-      return (Array width tau,
-              Array (MatrixShape.Householder order height width) qr)
-
-householderDeterminant ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
-   Householder height width a -> a
-householderDeterminant
-      (Array (MatrixShape.Householder order height width) a) =
-   let m = Shape.size height
-       n = Shape.size width
-       k = case order of RowMajor -> n; ColumnMajor -> m
-   in unsafePerformIO $
-      withForeignPtr a $ \aPtr ->
-         foldM (\x ptr -> do y <- peek ptr; return $! mul x y) one $
-         take (min m n) $ pointerSeq (k+1) aPtr
-
-newtype Mul a = Mul {getMul :: a -> a -> a}
-
-mul :: (Class.Floating a) => a -> a -> a
-mul = getMul $ Class.switchFloating (Mul (*)) (Mul (*)) (Mul (*)) (Mul (*))
-
-
-{-|
-Generalized determinant - works also for non-square matrices.
-In contrast to the square root of the Gramian determinant
-it has the proper sign.
--}
-determinant ::
-   (Shape.C height, Shape.C width, Eq a, Storable a, Class.Floating a) =>
-   General height width a -> a
-determinant a =
-   let (tau,hh) = householderDecompose a
-   in  foldl (\x _ -> neg x)
-         (householderDeterminant hh)
-         (takeWhile (/=zero) $ Array.toList tau)
-
-newtype Neg a = Neg {getNeg :: a -> a}
-
-neg :: (Class.Floating a) => a -> a
-neg =
-   getNeg $
-   Class.switchFloating (Neg negate) (Neg negate) (Neg negate) (Neg negate)
-
-
-householderExtractQ ::
-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>
-   (Vector width a, Householder height width a) -> General height height a
-householderExtractQ
-   (Array widthTau tau,
-    Array (MatrixShape.Householder order height width) qr) =
-
-   Array.unsafeCreate (MatrixShape.General order height height) $ \qPtr -> do
-
-   Call.assert "householderExtractQ: width shapes mismatch" (widthTau == width)
-
-   let m = Shape.size height
-   let k = min m $ Shape.size width
-   let lda = m
-   evalContT $ do
-      mPtr <- Call.cint m
-      kPtr <- Call.cint k
-      qrPtr <- ContT $ withForeignPtr qr
-      ldaPtr <- Call.cint lda
-      tauPtr <- ContT $ withForeignPtr tau
-      liftIO $
-         case order of
-            RowMajor -> do
-               copySubMatrix k m k qrPtr lda qPtr
-               withAutoWorkspaceInfo "unglq" $
-                  LapackGen.unglq mPtr mPtr kPtr qPtr ldaPtr tauPtr
-            ColumnMajor -> do
-               copyBlock (m*k) qrPtr qPtr
-               withAutoWorkspaceInfo "ungqr" $
-                  LapackGen.ungqr mPtr mPtr kPtr qPtr ldaPtr tauPtr
-
-householderExtractR ::
-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>
-   Householder height width a -> General height width a
-householderExtractR
-      (Array (MatrixShape.Householder order height width) qr) =
-
-   Array.unsafeCreate (MatrixShape.General order height width) $
-      \rPtr -> do
-
-   let (uplo, (m,n)) =
-         case order of
-            RowMajor -> ('L', (Shape.size width, Shape.size height))
-            ColumnMajor -> ('U', (Shape.size height, Shape.size width))
-   fill zero (m*n) rPtr
-   evalContT $ do
-      uploPtr <- Call.char uplo
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      qrPtr <- ContT $ withForeignPtr qr
-      ldqrPtr <- Call.cint m
-      ldrPtr <- Call.cint m
-      liftIO $ LapackGen.lacpy uploPtr mPtr nPtr qrPtr ldqrPtr rPtr ldrPtr
-
-{- |
-For an m-by-n-matrix @a@ with m>=n
-this function computes an m-by-(m-n)-matrix @b@
-such that @Matrix.multiply (transpose b) a@ is a zero matrix.
-The function does not try to compensate a rank deficiency of @a@.
-That is, @a|||b@ has full rank if and only if @a@ has full rank.
-
-For full-rank matrices you might also call this @kernel@ or @nullspace@.
--}
-orthogonalComplement ::
-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>
-   General height width a -> General height ZeroInt a
-orthogonalComplement a =
-   dropColumns (Shape.size $ MatrixShape.generalWidth $ Array.shape a) $
-   Array.mapShape zeroIntWidth $ householderExtractQ $ householderDecompose a
-
-zeroIntWidth ::
-   (Shape.C width) =>
-   MatrixShape.General height width -> MatrixShape.General height ZeroInt
-zeroIntWidth (MatrixShape.General order height width) =
-   MatrixShape.General order height (zeroInt $ Shape.size width)
-
-
-
-withAutoWorkspaceInfo ::
-   (Storable a, Class.Floating a) =>
-   String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
-withAutoWorkspaceInfo name computation = evalContT $ do
-   infoPtr <- Call.alloca
-   liftIO $ withAutoWorkspace $ \workPtr lworkPtr ->
-      computation workPtr lworkPtr infoPtr
-   info <- liftIO $ fromIntegral <$> peek infoPtr
-   case compare info (0::Int) of
-      EQ -> return ()
-      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
-      GT -> error $ printf "%s: deficient rank %d" name info
-
-withAutoWorkspace ::
-   (Storable a, Class.Floating a) =>
-   (Ptr a -> Ptr CInt -> IO ()) -> IO ()
-withAutoWorkspace computation = evalContT $ do
-   lworkPtr <- Call.cint (-1)
-   lwork <- liftIO $ alloca $ \workPtr -> do
-      computation workPtr lworkPtr
-      ceilingSize <$> peek workPtr
-   workPtr <- Call.allocaArray lwork
-   liftIO $ poke lworkPtr $ fromIntegral lwork
-   liftIO $ computation workPtr lworkPtr
-
-
-copyToColumnMajor ::
-   (Storable a, Class.Floating a) =>
-   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
-copyToColumnMajor order m n aPtr bPtr =
-   case order of
-      RowMajor -> copyTransposed m n aPtr m bPtr
-      ColumnMajor -> copyBlock (m*n) aPtr bPtr
-
-copyToSubColumnMajor ::
-   (Storable a, Class.Floating a) =>
-   Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
-copyToSubColumnMajor order m n aPtr ldb bPtr =
-   case order of
-      RowMajor -> copyTransposed m n aPtr ldb bPtr
-      ColumnMajor ->
-         if m==ldb
-           then copyBlock (m*n) aPtr bPtr
-           else copySubMatrix m n m aPtr ldb bPtr
-
-
-newtype FuncArg b a = FuncArg {runFuncArg :: a -> b}
-
-ceilingSize :: (Class.Floating a) => a -> Int
-ceilingSize =
-   runFuncArg $
-   Class.switchFloating
-      (FuncArg ceiling)
-      (FuncArg ceiling)
-      (FuncArg $ ceiling . Complex.realPart)
-      (FuncArg $ ceiling . Complex.realPart)
diff --git a/src/Numeric/LAPACK/Matrix.hs b/src/Numeric/LAPACK/Matrix.hs
--- a/src/Numeric/LAPACK/Matrix.hs
+++ b/src/Numeric/LAPACK/Matrix.hs
@@ -2,11 +2,12 @@
 {-# LANGUAGE TypeOperators #-}
 module Numeric.LAPACK.Matrix (
    General,
-   (Format.##),
-   Format.Format,
-   Format.FormatArray,
+   (##),
+   Format,
+   FormatArray,
    ZeroInt, zeroInt,
-   transpose,
+   transpose, adjoint,
+   fromScalar, toScalar,
    fromList,
    identity,
    diagonal, getDiagonal,
@@ -27,15 +28,21 @@
    multiply,
    multiplyVector,
 
-   trace,
+   Multiply, (<#>),
+   MultiplyLeft, (<#),
+   MultiplyRight, (#>),
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
-import qualified Numeric.LAPACK.Private as Private
-import qualified Numeric.LAPACK.Format as Format
+import qualified Numeric.LAPACK.Matrix.Square as Square
 import qualified Numeric.LAPACK.Vector as Vector
-import Numeric.LAPACK.Matrix.Shape.Private
-         (Order(RowMajor, ColumnMajor), charFromOrder)
+import Numeric.LAPACK.Format (Format, FormatArray, (##))
+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))
+import Numeric.LAPACK.Matrix.Multiply
+         (Multiply((<#>)), MultiplyLeft((<#)), MultiplyRight((#>)),
+          transpose, multiplyVector, multiply, multiplyVectorUnchecked)
+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)
+import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Private
          (zero, one, pointerSeq, copyTransposed, copySubMatrix, copyBlock)
 
@@ -64,18 +71,21 @@
 import Data.Bool.HT (if')
 
 
-type General height width = Array (MatrixShape.General height width)
-
-
-transpose :: General height width a -> General width height a
-transpose = Array.mapShape MatrixShape.transpose
+{- |
+conjugate transpose
+-}
+adjoint ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> General width height a
+adjoint = transpose . Vector.conjugate
 
 
-type ZeroInt = Shape.ZeroBased Int
-
-zeroInt :: Int -> ZeroInt
-zeroInt = Shape.ZeroBased
+fromScalar :: (Storable a) => a -> General () () a
+fromScalar = Square.toGeneral . Square.fromScalar
 
+toScalar :: (Storable a) => General () () a -> a
+toScalar (Array (MatrixShape.General _ () ()) a) =
+   unsafePerformIO $ withForeignPtr a peek
 
 fromList ::
    (Shape.C height, Shape.C width, Storable a) =>
@@ -84,70 +94,20 @@
    Array.fromList (MatrixShape.General RowMajor height width)
 
 
-type Vector = Array
-
-
-identity, _identity ::
-   (Shape.C sh, Storable a, Class.Floating a) =>
+identity ::
+   (Shape.C sh, Class.Floating a) =>
    sh -> General sh sh a
-identity sh =
-   Array.unsafeCreate (MatrixShape.General ColumnMajor sh sh) $ \aPtr ->
-   evalContT $ do
-      uploPtr <- Call.char 'A'
-      nPtr <- Call.cint $ Shape.size sh
-      alphaPtr <- Call.number zero
-      betaPtr <- Call.number one
-      liftIO $ LapackGen.laset uploPtr nPtr nPtr alphaPtr betaPtr aPtr nPtr
-
-_identity sh =
-   Array.unsafeCreate (MatrixShape.General ColumnMajor sh sh) $ \yPtr ->
-   evalContT $ do
-      nPtr <- Call.alloca
-      xPtr <- Call.number zero
-      incxPtr <- Call.cint 0
-      incyPtr <- Call.cint 1
-      liftIO $ do
-         let n = fromIntegral $ Shape.size sh
-         poke nPtr $ n*n
-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
-         poke nPtr n
-         poke xPtr one
-         poke incyPtr (n+1)
-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+identity = Square.toGeneral . Square.identity
 
 diagonal ::
-   (Shape.C sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Class.Floating a) =>
    Vector sh a -> General sh sh a
-diagonal (Array sh x) =
-   Array.unsafeCreate (MatrixShape.General ColumnMajor sh sh) $ \yPtr ->
-   evalContT $ do
-      nPtr <- Call.alloca
-      xPtr <- ContT $ withForeignPtr x
-      zPtr <- Call.number zero
-      incxPtr <- Call.cint 1
-      incyPtr <- Call.cint 1
-      inczPtr <- Call.cint 0
-      liftIO $ do
-         let n = fromIntegral $ Shape.size sh
-         poke nPtr $ n*n
-         BlasGen.copy nPtr zPtr inczPtr yPtr incyPtr
-         poke nPtr n
-         poke incyPtr (n+1)
-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+diagonal = Square.toGeneral . Square.diagonal
 
 getDiagonal ::
-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Eq sh, Class.Floating a) =>
    General sh sh a -> Vector sh a
-getDiagonal (Array (MatrixShape.General _ height width) x) =
-      Array.unsafeCreate height $ \yPtr -> do
-   Call.assert "getDiagonal: non-square matrix" (height==width)
-   evalContT $ do
-      let n = Shape.size height
-      nPtr <- Call.cint n
-      xPtr <- ContT $ withForeignPtr x
-      incxPtr <- Call.cint (n+1)
-      incyPtr <- Call.cint 1
-      liftIO $ BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+getDiagonal = Square.getDiagonal . Square.fromGeneral
 
 
 singleRow :: Vector width a -> General () width a
@@ -209,7 +169,7 @@
 
 pickRow ::
    (Shape.C height, Shape.C width, Shape.Index height ~ ix,
-    Storable a, Class.Floating a) =>
+    Class.Floating a) =>
    General height width a -> ix -> Vector width a
 pickRow (Array (MatrixShape.General order height width) x) ix =
    case order of
@@ -218,7 +178,7 @@
 
 pickColumn ::
    (Shape.C height, Shape.C width, Shape.Index width ~ ix,
-    Storable a, Class.Floating a) =>
+    Class.Floating a) =>
    General height width a -> ix -> Vector height a
 pickColumn (Array (MatrixShape.General order height width) x) ix =
    case order of
@@ -227,11 +187,10 @@
 
 pickConsecutive ::
    (Shape.C height, Shape.C width, Shape.Index height ~ ix,
-    Storable a, Class.Floating a) =>
+    Class.Floating a) =>
    height -> width -> ForeignPtr a -> ix -> Vector width a
 pickConsecutive height width x ix =
-   Array.unsafeCreate width $ \yPtr -> evalContT $ do
-      let n = Shape.size width
+   Array.unsafeCreateWithSize width $ \n yPtr -> evalContT $ do
       let offset = Shape.offset height ix
       nPtr <- Call.cint n
       xPtr <- ContT $ withForeignPtr x
@@ -242,11 +201,10 @@
 
 pickScattered ::
    (Shape.C height, Shape.C width, Shape.Index width ~ ix,
-    Storable a, Class.Floating a) =>
+    Class.Floating a) =>
    height -> width -> ForeignPtr a -> ix -> Vector height a
 pickScattered height width x ix =
-   Array.unsafeCreate height $ \yPtr -> evalContT $ do
-      let n = Shape.size height
+   Array.unsafeCreateWithSize height $ \n yPtr -> evalContT $ do
       let offset = Shape.offset width ix
       nPtr <- Call.cint n
       xPtr <- ContT $ withForeignPtr x
@@ -257,18 +215,19 @@
 
 
 takeRows, dropRows ::
-   (Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C width, Class.Floating a) =>
    Int -> General ZeroInt width a -> General ZeroInt width a
 takeRows k
       (Array (MatrixShape.General order (Shape.ZeroBased heightA) width) a) =
    let heightB = min k heightA
        n = Shape.size width
    in if' (k<0) (error "take: negative number") $
-      Array.unsafeCreate
-         (MatrixShape.General order (Shape.ZeroBased heightB) width) $ \bPtr ->
+      Array.unsafeCreateWithSize
+         (MatrixShape.General order (Shape.ZeroBased heightB) width) $
+            \blockSize bPtr ->
       withForeignPtr a $ \aPtr ->
       case order of
-         RowMajor -> copyBlock (heightB*n) aPtr bPtr
+         RowMajor -> copyBlock blockSize aPtr bPtr
          ColumnMajor -> copySubMatrix heightB n heightA aPtr heightB bPtr
 
 dropRows k0
@@ -277,17 +236,18 @@
        heightB = heightA - k
        n = Shape.size width
    in if' (k<0) (error "take: negative number") $
-      Array.unsafeCreate
-         (MatrixShape.General order (Shape.ZeroBased heightB) width) $ \bPtr ->
+      Array.unsafeCreateWithSize
+         (MatrixShape.General order (Shape.ZeroBased heightB) width) $
+            \blockSize bPtr ->
       withForeignPtr a $ \aPtr ->
       case order of
-         RowMajor -> copyBlock (heightB*n) (advancePtr aPtr (k*n)) bPtr
+         RowMajor -> copyBlock blockSize (advancePtr aPtr (k*n)) bPtr
          ColumnMajor ->
             copySubMatrix heightB n heightA (advancePtr aPtr k) heightB bPtr
 
 
 takeColumns, dropColumns ::
-   (Shape.C height, Storable a, Class.Floating a) =>
+   (Shape.C height, Class.Floating a) =>
    Int -> General height ZeroInt a -> General height ZeroInt a
 takeColumns k = transpose . takeRows k . transpose
 dropColumns k = transpose . dropRows k . transpose
@@ -295,10 +255,10 @@
 
 -- alternative: laswp
 reverseRows ::
-   (Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C width, Class.Floating a) =>
    General ZeroInt width a -> General ZeroInt width a
 reverseRows (Array shape@(MatrixShape.General order height width) a) =
-   Array.unsafeCreate shape $ \bPtr -> evalContT $ do
+   Array.unsafeCreateWithSize shape $ \blockSize bPtr -> evalContT $ do
       let n = Shape.size height
       let m = Shape.size width
       fwdPtr <- Call.bool True
@@ -307,26 +267,26 @@
       kPtr <- Call.allocaArray n
       aPtr <- ContT $ withForeignPtr a
       liftIO $ do
-         copyBlock (n*m) aPtr bPtr
+         copyBlock blockSize aPtr bPtr
          pokeArray kPtr $ take n $ iterate (subtract 1) $ fromIntegral n
          case order of
             RowMajor -> LapackGen.lapmt fwdPtr mPtr nPtr bPtr mPtr kPtr
             ColumnMajor -> LapackGen.lapmr fwdPtr nPtr mPtr bPtr nPtr kPtr
 
 reverseColumns ::
-   (Shape.C height, Storable a, Class.Floating a) =>
+   (Shape.C height, Class.Floating a) =>
    General height ZeroInt a -> General height ZeroInt a
 reverseColumns = transpose . reverseRows . transpose
 
 
 fromRowMajor ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    Array (height,width) a -> General height width a
 fromRowMajor (Array (height,width) x) =
    Array (MatrixShape.General RowMajor height width) x
 
 toRowMajor ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    General height width a -> Array (height,width) a
 toRowMajor (Array (MatrixShape.General order height width) x) =
    let shape = (height, width)
@@ -347,12 +307,11 @@
                   (pointerSeq n yPtr)
 
 flatten ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    General height width a -> Vector ZeroInt a
 flatten x =
    case toRowMajor x of
-      Array (height,width) fptr ->
-         Array (zeroInt $ Shape.size height * Shape.size width) fptr
+      Array shape fptr -> Array (zeroInt $ Shape.size shape) fptr
 
 
 infixl 3 |||
@@ -360,7 +319,7 @@
 
 (|||) ::
    (Shape.C height, Eq height, Shape.C widtha, Shape.C widthb,
-    Storable a, Class.Floating a) =>
+    Class.Floating a) =>
    General height widtha a ->
    General height widthb a ->
    General height (widtha:+:widthb) a
@@ -440,7 +399,7 @@
 
 (===) ::
    (Shape.C width, Eq width, Shape.C heighta, Shape.C heightb,
-    Storable a, Class.Floating a) =>
+    Class.Floating a) =>
    General heighta width a ->
    General heightb width a ->
    General (heighta:+:heightb) width a
@@ -448,93 +407,22 @@
 
 
 rowSums ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    General height width a -> Vector height a
 rowSums m =
    let MatrixShape.General _ _ width = Array.shape m
    in  multiplyVectorUnchecked m (Vector.constant width one)
 
 columnSums ::
-   (Shape.C height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Class.Floating a) =>
    General height width a -> Vector width a
 columnSums m =
    let MatrixShape.General _ height _ = Array.shape m
    in  multiplyVectorUnchecked (transpose m) (Vector.constant height one)
 
-multiplyVector ::
-   (Shape.C height, Shape.C width, Eq width,
-    Storable a, Class.Floating a) =>
-   General height width a -> Vector width a -> Vector height a
-multiplyVector a x =
-   let MatrixShape.General _order _height width = Array.shape a
-   in if width == Array.shape x
-         then multiplyVectorUnchecked a x
-         else error "multiplyVector: width shapes mismatch"
 
-multiplyVectorUnchecked ::
-   (Shape.C height, Shape.C width,
-    Storable a, Class.Floating a) =>
-   General height width a -> Vector width a -> Vector height a
-multiplyVectorUnchecked
-   (Array shape@(MatrixShape.General order height _width) a) (Array _ x) =
-      Array.unsafeCreate height $ \yPtr -> do
-   let (m,n) = MatrixShape.dimensions shape
-   let lda = m
-   evalContT $ do
-      transPtr <- Call.char $ charFromOrder order
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      alphaPtr <- Call.number one
-      aPtr <- ContT $ withForeignPtr a
-      ldaPtr <- Call.cint lda
-      xPtr <- ContT $ withForeignPtr x
-      incxPtr <- Call.cint 1
-      betaPtr <- Call.number zero
-      incyPtr <- Call.cint 1
-      liftIO $
-         BlasGen.gemv
-            transPtr mPtr nPtr alphaPtr aPtr ldaPtr
-            xPtr incxPtr betaPtr yPtr incyPtr
-
-multiply ::
-   (Shape.C height,
-    Shape.C fuse, Eq fuse,
-    Shape.C width,
-    Storable a, Class.Floating a) =>
-   General height fuse a -> General fuse width a -> General height width a
-multiply
-   (Array (MatrixShape.General orderA height fuseA) a)
-   (Array (MatrixShape.General orderB fuseB width) b) =
-      Array.unsafeCreate (MatrixShape.General ColumnMajor height width) $
-         \cPtr -> do
-   Call.assert "multiply: fuse shapes mismatch" (fuseA == fuseB)
-   let m = Shape.size height
-   let n = Shape.size width
-   let k = Shape.size fuseA
-   let lda = case orderA of RowMajor -> k; ColumnMajor -> m
-   let ldb = case orderB of RowMajor -> n; ColumnMajor -> k
-   let ldc = m
-   evalContT $ do
-      transaPtr <- Call.char $ charFromOrder orderA
-      transbPtr <- Call.char $ charFromOrder orderB
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      kPtr <- Call.cint k
-      alphaPtr <- Call.number one
-      aPtr <- ContT $ withForeignPtr a
-      ldaPtr <- Call.cint lda
-      bPtr <- ContT $ withForeignPtr b
-      ldbPtr <- Call.cint ldb
-      betaPtr <- Call.number zero
-      ldcPtr <- Call.cint ldc
-      liftIO $
-         BlasGen.gemm
-            transaPtr transbPtr mPtr nPtr kPtr alphaPtr aPtr ldaPtr
-            bPtr ldbPtr betaPtr cPtr ldcPtr
-
-
 scaleRows ::
-   (Shape.C height, Eq height, Shape.C width, Storable a, Class.Floating a) =>
+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
    Vector height a -> General height width a -> General height width a
 scaleRows
    (Array heightX x) (Array shape@(MatrixShape.General order height width) a) =
@@ -583,14 +471,6 @@
                (pointerSeq n bPtr)
 
 scaleColumns ::
-   (Shape.C height, Shape.C width, Eq width, Storable a, Class.Floating a) =>
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
    Vector width a -> General height width a -> General height width a
 scaleColumns x = transpose . scaleRows x . transpose
-
-
-
-trace :: (Shape.C sh, Eq sh, Class.Floating a) => General sh sh a -> a
-trace (Array (MatrixShape.General _ height width) x) = unsafePerformIO $ do
-   Call.assert "trace: non-square matrix" (height==width)
-   let n = Shape.size height
-   withForeignPtr x $ \xPtr -> Private.sum n xPtr (n+1)
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian.hs b/src/Numeric/LAPACK/Matrix/Hermitian.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Hermitian.hs
@@ -0,0 +1,540 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Matrix.Hermitian (
+   Hermitian,
+   fromList,
+   autoFromList,
+   identity,
+   diagonal,
+   getDiagonal,
+
+   multiplyVector,
+   square,
+   multiplySquareLeft,
+   multiplyGeneralLeft,
+   multiplySquareRight,
+   multiplyGeneralRight,
+   outer,
+   sumRank1,
+   sumRank2,
+
+   toSquare,
+   covariance,
+   addTransposed,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Triangular.Private
+         (forPointers, pack, unpack, unpackToTemp,
+          diagonalPointers, rowMajorPointers, columnMajorPointers)
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder)
+import Numeric.LAPACK.Matrix.Square (Square)
+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private
+         (RealOf, fill, zero, one, lacgv, fromReal, realPart, copyToTemp)
+
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
+import qualified Numeric.BLAS.FFI.Generic as BlasGen
+import qualified Numeric.BLAS.FFI.Complex as BlasComplex
+import qualified Numeric.BLAS.FFI.Real as BlasReal
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.C.Types (CInt, CChar)
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (Storable, poke, peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Monad (when)
+
+import qualified Data.NonEmpty as NonEmpty
+import Data.Foldable (forM_)
+import Data.Complex (Complex)
+
+
+type Hermitian sh = Array (MatrixShape.Hermitian sh)
+
+
+fromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Hermitian sh a
+fromList order sh =
+   Array.fromList (MatrixShape.Hermitian order sh)
+
+autoFromList :: (Storable a) => Order -> [a] -> Hermitian ZeroInt a
+autoFromList order xs =
+   fromList order
+      (zeroInt $ MatrixShape.triangleExtent "Hermitian.autoFromList" $
+       length xs)
+      xs
+
+
+identity :: (Shape.C sh, Class.Floating a) => Order -> sh -> Hermitian sh a
+identity order sh =
+   Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $
+      \triSize aPtr -> do
+   fill zero triSize aPtr
+   mapM_ (flip poke one . snd) $
+      diagonalPointers order (Shape.size sh) aPtr aPtr
+
+diagonal ::
+   (Shape.C sh, Class.Floating a) =>
+   Order -> Vector sh (RealOf a) -> Hermitian sh a
+diagonal order =
+   runDiagonal $
+   Class.switchFloating
+      (Diagonal $ diagonalAux order) (Diagonal $ diagonalAux order)
+      (Diagonal $ diagonalAux order) (Diagonal $ diagonalAux order)
+
+newtype Diagonal sh a =
+   Diagonal {runDiagonal :: Vector sh (RealOf a) -> Hermitian sh a}
+
+diagonalAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Order -> Vector sh ar -> Hermitian sh a
+diagonalAux order (Array sh x) =
+   Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $
+      \triSize aPtr -> do
+   fill zero triSize aPtr
+   withForeignPtr x $ \xPtr ->
+      forM_ (diagonalPointers order (Shape.size sh) xPtr aPtr) $
+         \(srcPtr,dstPtr) -> poke dstPtr . fromReal =<< peek srcPtr
+
+
+getDiagonal ::
+   (Shape.C sh, Class.Floating a) =>
+   Hermitian sh a -> Vector sh (RealOf a)
+getDiagonal =
+   runGetDiagonal $
+   Class.switchFloating
+      (GetDiagonal $ getDiagonalAux) (GetDiagonal $ getDiagonalAux)
+      (GetDiagonal $ getDiagonalAux) (GetDiagonal $ getDiagonalAux)
+
+newtype GetDiagonal sh a =
+   GetDiagonal {runGetDiagonal :: Hermitian sh a -> Vector sh (RealOf a)}
+
+getDiagonalAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Hermitian sh a -> Vector sh ar
+getDiagonalAux (Array (MatrixShape.Hermitian order sh) a) =
+   Array.unsafeCreateWithSize sh $ \n xPtr ->
+   withForeignPtr a $ \aPtr ->
+      forM_ (diagonalPointers order n xPtr aPtr) $
+         \(dstPtr,srcPtr) -> poke dstPtr . realPart =<< peek srcPtr
+
+
+multiplyVector ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   Hermitian sh a -> Vector sh a -> Vector sh a
+multiplyVector (Array (MatrixShape.Hermitian order shA) a) (Array shX x) =
+      Array.unsafeCreateWithSize shX $ \n yPtr -> do
+   Call.assert "Hermitian.multiplyVector: width shapes mismatch" (shA == shX)
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder order
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr a
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      betaPtr <- Call.number zero
+      incyPtr <- Call.cint 1
+      liftIO $
+         BlasGen.hpmv
+            uploPtr nPtr alphaPtr aPtr xPtr incxPtr betaPtr yPtr incyPtr
+
+
+square ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   Hermitian sh a -> Hermitian sh a
+square
+   (Array shape@(MatrixShape.Hermitian order sh) a) =
+      Array.unsafeCreate shape $ \cpPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      sidePtr <- Call.char 'L'
+      uploPtr <- Call.char 'U'
+      nPtr <- Call.cint n
+      let ldPtr = nPtr
+      bPtr <- unpackToTemp (unpackFull order) n a
+      cPtr <- Call.allocaArray (n*n)
+      alphaPtr <- Call.number one
+      betaPtr <- Call.number zero
+      liftIO $ do
+         BlasGen.hemm sidePtr uploPtr
+            nPtr nPtr alphaPtr bPtr ldPtr
+            bPtr ldPtr betaPtr cPtr ldPtr
+         pack order n cPtr cpPtr
+
+
+multiplySquareLeft ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   Square sh a -> Hermitian sh a -> Square sh a
+multiplySquareLeft
+   (Array shapeB@(MatrixShape.Square orderB shB) b)
+   (Array        (MatrixShape.Hermitian orderA shA) a) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Hermitian.multiplySquareLeft: shapes mismatch" (shA == shB)
+   let n = Shape.size shB
+   multiplyAux True orderA n a (flipOrder orderB) n b cPtr
+
+multiplyGeneralLeft ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a -> Hermitian width a -> General height width a
+multiplyGeneralLeft
+   (Array shapeB@(MatrixShape.General orderB height width) b)
+   (Array        (MatrixShape.Hermitian orderA shA) a) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Hermitian.multiplyGeneralLeft: shapes mismatch" (shA == width)
+   multiplyAux True
+      orderA (Shape.size width) a (flipOrder orderB) (Shape.size height) b cPtr
+
+multiplySquareRight ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   Hermitian sh a -> Square sh a -> Square sh a
+multiplySquareRight
+   (Array        (MatrixShape.Hermitian orderA shA) a)
+   (Array shapeB@(MatrixShape.Square orderB shB) b) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Hermitian.multiplySquareRight: shapes mismatch" (shA == shB)
+   let n = Shape.size shB
+   multiplyAux False orderA n a orderB n b cPtr
+
+multiplyGeneralRight ::
+   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
+   Hermitian height a -> General height width a -> General height width a
+multiplyGeneralRight
+   (Array        (MatrixShape.Hermitian orderA shA) a)
+   (Array shapeB@(MatrixShape.General orderB height width) b) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Hermitian.multiplyGeneralRight: shapes mismatch" (shA == height)
+   multiplyAux False
+      orderA (Shape.size height) a orderB (Shape.size width) b cPtr
+
+multiplyAux ::
+   Class.Floating a =>
+   Bool ->
+   Order -> Int -> ForeignPtr a ->
+   Order -> Int -> ForeignPtr a -> Ptr a -> IO ()
+multiplyAux extraConjugate orderA m0 a orderB n0 b cPtr = do
+   let size = m0*m0
+   evalContT $ do
+      let (side,(m,n)) =
+            case orderB of
+               ColumnMajor -> ('L',(m0,n0))
+               RowMajor -> ('R',(n0,m0))
+      sidePtr <- Call.char side
+      uploPtr <- Call.char $ uploFromOrder orderA
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number one
+      aPtr <- unpackToTemp (unpack orderA) m0 a
+      ldaPtr <- Call.cint m0
+      incaPtr <- Call.cint 1
+      sizePtr <- Call.cint size
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint m
+      betaPtr <- Call.number zero
+      ldcPtr <- Call.cint m
+      liftIO $ do
+         when ((orderA/=orderB) /= extraConjugate) $
+            lacgv sizePtr aPtr incaPtr
+         BlasGen.hemm sidePtr uploPtr
+            mPtr nPtr alphaPtr aPtr ldaPtr
+            bPtr ldbPtr betaPtr cPtr ldcPtr
+
+
+outer :: (Shape.C sh, Class.Floating a) => Vector sh a -> Hermitian sh a
+outer =
+   getMap $
+   Class.switchFloating
+      (Map outerAux) (Map outerAux)
+      (Map outerAux) (Map outerAux)
+
+outerAux ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   Vector sh a -> Hermitian sh a
+outerAux (Array sh x) =
+   Array.unsafeCreateWithSize (MatrixShape.Hermitian ColumnMajor sh) $
+      \triSize aPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder ColumnMajor
+      nPtr <- Call.cint n
+      alphaPtr <- Call.real one
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      liftIO $ fill zero triSize aPtr
+      liftIO $ hpr uploPtr nPtr alphaPtr xPtr incxPtr aPtr
+
+
+sumRank1 ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   NonEmpty.T [] (RealOf a, Vector sh a) -> Hermitian sh a
+sumRank1 =
+   getSumRank1 $
+   Class.switchFloating
+      (SumRank1 sumRank1Aux) (SumRank1 sumRank1Aux)
+      (SumRank1 sumRank1Aux) (SumRank1 sumRank1Aux)
+
+type SumRank1_ sh a = NonEmpty.T [] (RealOf a, Vector sh a) -> Hermitian sh a
+
+newtype SumRank1 sh a = SumRank1 {getSumRank1 :: SumRank1_ sh a}
+
+sumRank1Aux ::
+   (Shape.C sh, Eq sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   SumRank1_ sh a
+sumRank1Aux xs@(NonEmpty.Cons (_, Array sh _) _) =
+   Array.unsafeCreateWithSize (MatrixShape.Hermitian ColumnMajor sh) $
+      \triSize aPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder ColumnMajor
+      nPtr <- Call.cint n
+      alphaPtr <- Call.alloca
+      incxPtr <- Call.cint 1
+      liftIO $ do
+         fill zero triSize aPtr
+         forM_ xs $ \(alpha, Array shX x) ->
+            withForeignPtr x $ \xPtr -> do
+               Call.assert
+                  "Hermitian.sumRank1: non-matching vector size" (sh==shX)
+               poke alphaPtr alpha
+               hpr uploPtr nPtr alphaPtr xPtr incxPtr aPtr
+
+
+type HPR_ a =
+   Ptr CChar -> Ptr CInt ->
+   Ptr (RealOf a) -> Ptr a -> Ptr CInt -> Ptr a -> IO ()
+
+newtype HPR a = HPR {getHPR :: HPR_ a}
+
+hpr :: Class.Floating a => HPR_ a
+hpr =
+   getHPR $
+   Class.switchFloating
+      (HPR BlasReal.spr) (HPR BlasReal.spr)
+      (HPR BlasComplex.hpr) (HPR BlasComplex.hpr)
+
+
+sumRank2 ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   NonEmpty.T [] (a, (Vector sh a, Vector sh a)) -> Hermitian sh a
+sumRank2 xys@(NonEmpty.Cons (_, (Array sh _, _)) _) =
+   Array.unsafeCreateWithSize (MatrixShape.Hermitian ColumnMajor sh) $
+      \triSize aPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder ColumnMajor
+      nPtr <- Call.cint n
+      alphaPtr <- Call.alloca
+      incPtr <- Call.cint 1
+      liftIO $ do
+         fill zero triSize aPtr
+         forM_ xys $ \(alpha, (Array shX x, Array shY y)) ->
+            withForeignPtr x $ \xPtr ->
+            withForeignPtr y $ \yPtr -> do
+               Call.assert
+                  "Hermitian.sumRank2: non-matching x vector size" (sh==shX)
+               Call.assert
+                  "Hermitian.sumRank2: non-matching y vector size" (sh==shY)
+               poke alphaPtr alpha
+               BlasGen.hpr2 uploPtr nPtr alphaPtr xPtr incPtr yPtr incPtr aPtr
+
+
+{-
+It is not strictly necessary to keep the 'order'.
+It would be neither more complicated nor less efficient
+to change the order via the conversion.
+-}
+toSquare, _toSquare ::
+   (Shape.C sh, Class.Floating a) => Hermitian sh a -> Square sh a
+_toSquare (Array (MatrixShape.Hermitian order sh) a) =
+      Array.unsafeCreate (MatrixShape.Square order sh) $ \bPtr ->
+   evalContT $ do
+      let n = Shape.size sh
+      aPtr <- ContT $ withForeignPtr a
+      conjPtr <- conjugateToTemp (MatrixShape.triangleSize n) a
+      liftIO $ do
+         unpack (flipOrder order) n conjPtr bPtr -- wrong
+         unpack order n aPtr bPtr
+
+toSquare (Array (MatrixShape.Hermitian order sh) a) =
+      Array.unsafeCreate (MatrixShape.Square order sh) $ \bPtr ->
+   withForeignPtr a $ \aPtr ->
+      unpackFull order (Shape.size sh) aPtr bPtr
+
+
+{- |
+Make a temporary copy only for complex matrices.
+-}
+conjugateToTemp ::
+   (Class.Floating a) => Int -> ForeignPtr a -> ContT r IO (Ptr a)
+conjugateToTemp n =
+   runCopyToTemp $
+   Class.switchFloating
+      (CopyToTemp $ ContT . withForeignPtr)
+      (CopyToTemp $ ContT . withForeignPtr)
+      (CopyToTemp $ complexConjugateToTemp n)
+      (CopyToTemp $ complexConjugateToTemp n)
+
+newtype CopyToTemp r a =
+   CopyToTemp {runCopyToTemp :: ForeignPtr a -> ContT r IO (Ptr a)}
+
+complexConjugateToTemp ::
+   Class.Real a =>
+   Int -> ForeignPtr (Complex a) -> ContT r IO (Ptr (Complex a))
+complexConjugateToTemp n x = do
+   nPtr <- Call.cint n
+   xPtr <- copyToTemp n x
+   incxPtr <- Call.cint 1
+   liftIO $ LapackComplex.lacgv nPtr xPtr incxPtr
+   return xPtr
+
+
+{- |
+A^H * A
+-}
+covariance ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> Hermitian width a
+covariance =
+   getMap $
+   Class.switchFloating
+      (Map covarianceAux) (Map covarianceAux)
+      (Map covarianceAux) (Map covarianceAux)
+
+newtype Map f g a = Map {getMap :: f a -> g a}
+
+covarianceAux ::
+   (Shape.C height, Shape.C width,
+    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   General height width a -> Hermitian width a
+covarianceAux (Array (MatrixShape.General order height width) a) =
+      Array.unsafeCreate (MatrixShape.Hermitian order width) $ \bPtr -> do
+   let n = Shape.size width
+   let k = Shape.size height
+   evalContT $ do
+      nPtr <- Call.cint n
+      kPtr <- Call.cint k
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr a
+      betaPtr <- Call.number zero
+      cPtr <- Call.allocaArray (n*n)
+      ldcPtr <- Call.cint n
+
+      case order of
+         ColumnMajor -> do
+            uploPtr <- Call.char 'U'
+            transPtr <- Call.char 'C'
+            ldaPtr <- Call.cint k
+            liftIO $ do
+               herk uploPtr transPtr
+                  nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr
+               pack ColumnMajor n cPtr bPtr
+
+         RowMajor -> do
+            uploPtr <- Call.char 'L'
+            transPtr <- Call.char 'N'
+            ldaPtr <- Call.cint n
+            liftIO $ do
+               herk uploPtr transPtr
+                  nPtr kPtr alphaPtr aPtr ldaPtr betaPtr cPtr ldcPtr
+               pack RowMajor n cPtr bPtr
+
+
+type HERK_ a =
+   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt -> Ptr (RealOf a) -> Ptr a ->
+   Ptr CInt -> Ptr (RealOf a) -> Ptr a -> Ptr CInt -> IO ()
+
+newtype HERK a = HERK {getHERK :: HERK_ a}
+
+herk :: Class.Floating a => HERK_ a
+herk =
+   getHERK $
+   Class.switchFloating
+      (HERK BlasReal.syrk)
+      (HERK BlasReal.syrk)
+      (HERK BlasComplex.herk)
+      (HERK BlasComplex.herk)
+
+
+{- |
+A^H + A
+-}
+addTransposed, _addTransposed ::
+   (Shape.C sh, Class.Floating a) => Square sh a -> Hermitian sh a
+_addTransposed (Array (MatrixShape.Square order sh) a) =
+      Array.unsafeCreateWithSize (MatrixShape.Hermitian order sh) $ \bSize bPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      alphaPtr <- Call.number one
+      incxPtr <- Call.cint 1
+      aPtr <- ContT $ withForeignPtr a
+      sizePtr <- Call.cint bSize
+      conjPtr <- Call.allocaArray bSize
+      liftIO $ do
+         pack order n aPtr bPtr
+         pack (flipOrder order) n aPtr conjPtr -- wrong
+         lacgv sizePtr conjPtr incxPtr
+         BlasGen.axpy sizePtr alphaPtr conjPtr incxPtr bPtr incxPtr
+
+addTransposed (Array (MatrixShape.Square order sh) a) =
+      Array.unsafeCreate (MatrixShape.Hermitian order sh) $ \bPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      alphaPtr <- Call.number one
+      incxPtr <- Call.cint 1
+      incnPtr <- Call.cint n
+      aPtr <- ContT $ withForeignPtr a
+      liftIO $ case order of
+         RowMajor ->
+            forPointers (rowMajorPointers n aPtr bPtr) $
+               \nPtr (srcPtr,dstPtr) -> do
+            BlasGen.copy nPtr srcPtr incnPtr dstPtr incxPtr
+            lacgv nPtr dstPtr incxPtr
+            BlasGen.axpy nPtr alphaPtr srcPtr incxPtr dstPtr incxPtr
+         ColumnMajor ->
+            forPointers (columnMajorPointers n aPtr bPtr) $
+               \nPtr ((srcRowPtr,srcColumnPtr),dstPtr) -> do
+            BlasGen.copy nPtr srcRowPtr incnPtr dstPtr incxPtr
+            lacgv nPtr dstPtr incxPtr
+            BlasGen.axpy nPtr alphaPtr srcColumnPtr incxPtr dstPtr incxPtr
+
+
+unpackFull :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
+unpackFull order n packedPtr fullPtr = evalContT $ do
+   incxPtr <- Call.cint 1
+   incyPtr <- Call.cint n
+   liftIO $ case order of
+      RowMajor ->
+         forPointers (rowMajorPointers n fullPtr packedPtr) $
+               \nPtr (dstPtr,srcPtr) -> do
+            BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr
+            lacgv nPtr dstPtr incyPtr
+            BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr
+      ColumnMajor ->
+         forPointers (columnMajorPointers n fullPtr packedPtr) $
+               \nPtr ((dstRowPtr,dstColumnPtr),srcPtr) -> do
+            BlasGen.copy nPtr srcPtr incxPtr dstRowPtr incyPtr
+            lacgv nPtr dstRowPtr incyPtr
+            BlasGen.copy nPtr srcPtr incxPtr dstColumnPtr incxPtr
+
+_pack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
+_pack order n fullPtr packedPtr =
+   evalContT $ do
+      incxPtr <- Call.cint 1
+      liftIO $
+         case order of
+            ColumnMajor ->
+               forPointers (columnMajorPointers n fullPtr packedPtr) $
+                  \nPtr ((_,srcPtr),dstPtr) ->
+                     BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr
+            RowMajor ->
+               forPointers (rowMajorPointers n fullPtr packedPtr) $
+                  \nPtr (srcPtr,dstPtr) ->
+                     BlasGen.copy nPtr srcPtr incxPtr dstPtr incxPtr
diff --git a/src/Numeric/LAPACK/Matrix/Multiply.hs b/src/Numeric/LAPACK/Matrix/Multiply.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Multiply.hs
@@ -0,0 +1,328 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module Numeric.LAPACK.Matrix.Multiply where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Vector as Vector
+import qualified Numeric.LAPACK.Private as Private
+import Numeric.LAPACK.Matrix.Shape.Private
+         (HeightOf, WidthOf, Order(ColumnMajor), transposeFromOrder)
+import Numeric.LAPACK.Matrix.Triangular (Triangular)
+import Numeric.LAPACK.Matrix.Private (General)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private (zero, one)
+
+import qualified Numeric.BLAS.FFI.Generic as BlasGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.ForeignPtr (withForeignPtr)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+
+transpose :: General height width a -> General width height a
+transpose = Array.mapShape MatrixShape.transpose
+
+multiplyVector ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a -> Vector width a -> Vector height a
+multiplyVector a x =
+   let MatrixShape.General _order _height width = Array.shape a
+   in if width == Array.shape x
+         then multiplyVectorUnchecked a x
+         else error "multiplyVector: width shapes mismatch"
+
+multiplyVectorUnchecked ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> Vector width a -> Vector height a
+multiplyVectorUnchecked
+   (Array shape@(MatrixShape.General order height _width) a) (Array _ x) =
+      Array.unsafeCreate height $ \yPtr -> do
+   let (m,n) = MatrixShape.dimensions shape
+   let lda = m
+   evalContT $ do
+      transPtr <- Call.char $ transposeFromOrder order
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr a
+      ldaPtr <- Call.cint lda
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      betaPtr <- Call.number zero
+      incyPtr <- Call.cint 1
+      liftIO $
+         BlasGen.gemv
+            transPtr mPtr nPtr alphaPtr aPtr ldaPtr
+            xPtr incxPtr betaPtr yPtr incyPtr
+
+multiply ::
+   (Shape.C height,
+    Shape.C fuse, Eq fuse,
+    Shape.C width,
+    Class.Floating a) =>
+   General height fuse a -> General fuse width a -> General height width a
+multiply
+   (Array (MatrixShape.General orderA height fuseA) a)
+   (Array (MatrixShape.General orderB fuseB width) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor height width) $
+         \cPtr -> do
+   Call.assert "multiply: fuse shapes mismatch" (fuseA == fuseB)
+   let m = Shape.size height
+   let n = Shape.size width
+   let k = Shape.size fuseA
+   Private.multiplyMatrix orderA orderB m k n a b cPtr
+
+
+infixl 7 <#, <#>
+infixr 7 #>
+
+class MultiplyRight shape where
+   (#>) ::
+      (Class.Floating a) =>
+      Array shape a -> Array (WidthOf shape) a -> Array (HeightOf shape) a
+
+class MultiplyLeft shape where
+   (<#) ::
+      (Class.Floating a) =>
+      Array (HeightOf shape) a -> Array shape a -> Array (WidthOf shape) a
+
+class Multiply shapeA shapeB where
+   type Multiplied shapeA shapeB
+   (<#>) ::
+      (Class.Floating a) =>
+      Array shapeA a -> Array shapeB a -> Array (Multiplied shapeA shapeB) a
+
+
+instance
+   (Eq width, Shape.C width, Shape.C height) =>
+      MultiplyRight (MatrixShape.General height width) where
+   (#>) = multiplyVector
+
+instance
+   (Eq height, Shape.C width, Shape.C height) =>
+      MultiplyLeft (MatrixShape.General height width) where
+   v <# m = multiplyVector (transpose m) v
+
+
+instance
+   (Eq shape, Shape.C shape) =>
+      MultiplyRight (MatrixShape.Square shape) where
+   m #> v = multiplyVector (Square.toGeneral m) v
+
+instance
+   (Eq shape, Shape.C shape) =>
+      MultiplyLeft (MatrixShape.Square shape) where
+   v <# m = multiplyVector (transpose $ Square.toGeneral m) v
+
+
+instance
+   (Eq shape, Shape.C shape) =>
+      MultiplyRight (MatrixShape.Hermitian shape) where
+   m #> v = Hermitian.multiplyVector m v
+
+instance
+   (Eq shape, Shape.C shape) =>
+      MultiplyLeft (MatrixShape.Hermitian shape) where
+   v <# m = Hermitian.multiplyVector (Vector.conjugate m) v
+
+
+instance
+   (MatrixShape.Uplo uplo, Eq shape, Shape.C shape) =>
+      MultiplyRight (MatrixShape.Triangular uplo shape) where
+   m #> v = Triangular.multiplyVectorRight m v
+
+instance
+   (MatrixShape.Uplo uplo, Eq shape, Shape.C shape) =>
+      MultiplyLeft (MatrixShape.Triangular uplo shape) where
+   v <# m = Triangular.multiplyVectorLeft m v
+
+
+instance
+   (Shape.C heightA, Shape.C widthA, Shape.C widthB,
+    widthA ~ heightB, Eq heightB) =>
+      Multiply
+         (MatrixShape.General heightA widthA)
+         (MatrixShape.General heightB widthB) where
+   type Multiplied
+         (MatrixShape.General heightA widthA)
+         (MatrixShape.General heightB widthB) =
+            MatrixShape.General heightA widthB
+   (<#>) = multiply
+
+instance
+   (Shape.C shapeA, Shape.C widthB, shapeA ~ heightB, Eq heightB) =>
+      Multiply
+         (MatrixShape.Square shapeA)
+         (MatrixShape.General heightB widthB) where
+   type Multiplied
+         (MatrixShape.Square shapeA)
+         (MatrixShape.General heightB widthB) =
+            MatrixShape.General heightB widthB
+   a <#> b = multiply (Square.toGeneral a) b
+
+instance
+   (Shape.C heightA, Shape.C widthA, widthA ~ shapeB, Eq shapeB) =>
+      Multiply
+         (MatrixShape.General heightA widthA)
+         (MatrixShape.Square shapeB) where
+   type Multiplied
+         (MatrixShape.General heightA widthA)
+         (MatrixShape.Square shapeB) =
+            MatrixShape.General heightA widthA
+   a <#> b = multiply a (Square.toGeneral b)
+
+instance
+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>
+      Multiply (MatrixShape.Square shapeA) (MatrixShape.Square shapeB) where
+   type Multiplied (MatrixShape.Square shapeA) (MatrixShape.Square shapeB) =
+            MatrixShape.Square shapeA
+   (<#>) = Square.multiply
+
+
+instance
+   (Shape.C shapeA, shapeA ~ width, Eq width, Shape.C height) =>
+      Multiply
+         (MatrixShape.General height width)
+         (MatrixShape.Hermitian shapeA)
+            where
+   type Multiplied
+         (MatrixShape.General height width) (MatrixShape.Hermitian shapeA) =
+            MatrixShape.General height width
+   (<#>) = Hermitian.multiplyGeneralLeft
+
+instance
+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>
+      Multiply (MatrixShape.Square shapeB) (MatrixShape.Hermitian shapeA)
+         where
+   type Multiplied
+         (MatrixShape.Square shapeB) (MatrixShape.Hermitian shapeA) =
+            MatrixShape.Square shapeA
+   (<#>) = Hermitian.multiplySquareLeft
+
+instance
+   (Shape.C shapeA, shapeA ~ height, Eq height, Shape.C width) =>
+      Multiply
+         (MatrixShape.Hermitian shapeA)
+         (MatrixShape.General height width)
+            where
+   type Multiplied
+         (MatrixShape.Hermitian shapeA) (MatrixShape.General height width) =
+            MatrixShape.General height width
+   (<#>) = Hermitian.multiplyGeneralRight
+
+instance
+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>
+      Multiply (MatrixShape.Hermitian shapeA) (MatrixShape.Square shapeB)
+         where
+   type Multiplied
+         (MatrixShape.Hermitian shapeA) (MatrixShape.Square shapeB) =
+            MatrixShape.Square shapeA
+   (<#>) = Hermitian.multiplySquareRight
+
+instance
+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>
+      Multiply (MatrixShape.Hermitian shapeA) (MatrixShape.Hermitian shapeB)
+         where
+   type Multiplied
+         (MatrixShape.Hermitian shapeA) (MatrixShape.Hermitian shapeB) =
+            MatrixShape.Square shapeA
+   a <#> b = Hermitian.multiplySquareRight a (Hermitian.toSquare b)
+
+
+
+
+instance
+   (MatrixShape.Uplo uplo,
+    Shape.C shapeA, shapeA ~ width, Eq width, Shape.C height) =>
+      Multiply
+         (MatrixShape.General height width)
+         (MatrixShape.Triangular uplo shapeA)
+            where
+   type Multiplied
+         (MatrixShape.General height width)
+         (MatrixShape.Triangular uplo shapeA) =
+            MatrixShape.General height width
+   (<#>) = Triangular.multiplyGeneralLeft
+
+instance
+   (MatrixShape.Uplo uplo, Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>
+      Multiply (MatrixShape.Square shapeB) (MatrixShape.Triangular uplo shapeA)
+         where
+   type Multiplied
+         (MatrixShape.Square shapeB) (MatrixShape.Triangular uplo shapeA) =
+            MatrixShape.Square shapeA
+   (<#>) = Triangular.multiplySquareLeft
+
+instance
+   (MatrixShape.Uplo uplo,
+    Shape.C shapeA, shapeA ~ height, Eq height, Shape.C width) =>
+      Multiply
+         (MatrixShape.Triangular uplo shapeA)
+         (MatrixShape.General height width)
+            where
+   type Multiplied
+         (MatrixShape.Triangular uplo shapeA)
+         (MatrixShape.General height width) =
+            MatrixShape.General height width
+   (<#>) = Triangular.multiplyGeneralRight
+
+instance
+   (MatrixShape.Uplo uplo, Shape.C shapeA, shapeA ~ shapeB, Eq shapeB) =>
+      Multiply (MatrixShape.Triangular uplo shapeA) (MatrixShape.Square shapeB)
+         where
+   type Multiplied
+         (MatrixShape.Triangular uplo shapeA) (MatrixShape.Square shapeB) =
+            MatrixShape.Square shapeA
+   (<#>) = Triangular.multiplySquareRight
+
+instance
+   (Shape.C shapeA, shapeA ~ shapeB, Eq shapeB,
+    MultiplyTriangular uploA uploB) =>
+      Multiply
+         (MatrixShape.Triangular uploA shapeA)
+         (MatrixShape.Triangular uploB shapeB) where
+   type Multiplied
+         (MatrixShape.Triangular uploA shapeA)
+         (MatrixShape.Triangular uploB shapeB) =
+            MultipliedTriangular uploA uploB shapeB
+   (<#>) = multiplyTriangular
+
+class MultiplyTriangular uploA uploB where
+   type MultipliedTriangular uploA uploB :: * -> *
+   multiplyTriangular ::
+      (Class.Floating a, Shape.C shape, Eq shape) =>
+      Triangular uploA shape a ->
+      Triangular uploB shape a ->
+      Array (MultipliedTriangular uploA uploB shape) a
+
+instance MultiplyTriangular MatrixShape.Lower MatrixShape.Lower where
+   type MultipliedTriangular MatrixShape.Lower MatrixShape.Lower =
+         MatrixShape.Triangular MatrixShape.Lower
+   multiplyTriangular = Triangular.multiply
+
+instance MultiplyTriangular MatrixShape.Upper MatrixShape.Upper where
+   type MultipliedTriangular MatrixShape.Upper MatrixShape.Upper =
+         MatrixShape.Triangular MatrixShape.Upper
+   multiplyTriangular = Triangular.multiply
+
+instance MultiplyTriangular MatrixShape.Lower MatrixShape.Upper where
+   type MultipliedTriangular MatrixShape.Lower MatrixShape.Upper =
+         MatrixShape.Square
+   multiplyTriangular a b =
+      Square.multiply (Triangular.toSquare a) (Triangular.toSquare b)
+
+instance MultiplyTriangular MatrixShape.Upper MatrixShape.Lower where
+   type MultipliedTriangular MatrixShape.Upper MatrixShape.Lower =
+         MatrixShape.Square
+   multiplyTriangular a b =
+      Square.multiply (Triangular.toSquare a) (Triangular.toSquare b)
diff --git a/src/Numeric/LAPACK/Matrix/Private.hs b/src/Numeric/LAPACK/Matrix/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Private.hs
@@ -0,0 +1,15 @@
+module Numeric.LAPACK.Matrix.Private where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Data.Array.Comfort.Storable (Array)
+
+import qualified Data.Array.Comfort.Shape as Shape
+
+
+type General height width = Array (MatrixShape.General height width)
+
+
+type ZeroInt = Shape.ZeroBased Int
+
+zeroInt :: Int -> ZeroInt
+zeroInt = Shape.ZeroBased
diff --git a/src/Numeric/LAPACK/Matrix/Shape/Private.hs b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Shape/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
@@ -3,6 +3,10 @@
 
 import qualified Data.Array.Comfort.Shape as Shape
 
+import Control.Applicative (Const(Const, getConst))
+
+import Data.Functor.Identity (Identity(Identity), runIdentity)
+import Data.List (tails)
 import Data.Tuple.HT (swap)
 
 
@@ -13,11 +17,15 @@
 flipOrder RowMajor = ColumnMajor
 flipOrder ColumnMajor = RowMajor
 
-charFromOrder :: Order -> Char
-charFromOrder RowMajor = 'T'
-charFromOrder ColumnMajor = 'N'
+transposeFromOrder :: Order -> Char
+transposeFromOrder RowMajor = 'T'
+transposeFromOrder ColumnMajor = 'N'
 
 
+type family HeightOf shape
+type family WidthOf shape
+
+
 data General height width =
    General {
       generalOrder :: Order,
@@ -25,6 +33,9 @@
       generalWidth :: width
    } deriving (Eq, Show)
 
+type instance HeightOf (General height width) = height
+type instance WidthOf (General height width) = width
+
 instance (Shape.C height, Shape.C width) => Shape.C (General height width) where
    type Index (General height width) = (Shape.Index height, Shape.Index width)
    indices (General _ height width) = Shape.indices (height,width)
@@ -63,6 +74,48 @@
       ColumnMajor -> (Shape.size height, Shape.size width)
 
 
+data Square size =
+   Square {
+      squareOrder :: Order,
+      squareSize :: size
+   } deriving (Eq, Show)
+
+type instance HeightOf (Square size) = size
+type instance WidthOf (Square size) = size
+
+generalFromSquare :: Square size -> General size size
+generalFromSquare (Square order sh) = General order sh sh
+
+transposeSquare :: Square sh -> Square sh
+transposeSquare (Square order size) = Square (flipOrder order) size
+
+instance (Shape.C size) => Shape.C (Square size) where
+   type Index (Square size) = (Shape.Index size, Shape.Index size)
+   indices (Square _ size) = Shape.indices (size,size)
+
+   offset (Square RowMajor size) =
+      Shape.offset (size,size)
+   offset (Square ColumnMajor size) =
+      Shape.offset (size,size) . swap
+   uncheckedOffset (Square RowMajor size) =
+      Shape.uncheckedOffset (size,size)
+   uncheckedOffset (Square ColumnMajor size) =
+      Shape.uncheckedOffset (size,size) . swap
+
+   sizeOffset (Square RowMajor size) =
+      Shape.sizeOffset (size,size)
+   sizeOffset (Square ColumnMajor size) =
+      Shape.sizeOffset (size,size) . swap
+   uncheckedSizeOffset (Square RowMajor size) =
+      Shape.uncheckedSizeOffset (size,size)
+   uncheckedSizeOffset (Square ColumnMajor size) =
+      Shape.uncheckedSizeOffset (size,size) . swap
+
+   inBounds (Square _ size) = Shape.inBounds (size,size)
+   size (Square _ size) = Shape.size (size,size)
+   uncheckedSize (Square _ size) = Shape.uncheckedSize (size,size)
+
+
 data Householder height width =
    Householder {
       householderOrder :: Order,
@@ -70,24 +123,27 @@
       householderWidth :: width
    } deriving (Eq, Show)
 
+type instance HeightOf (Householder height width) = height
+type instance WidthOf (Householder height width) = width
+
 data Reflector = Reflector deriving (Eq)
-data Triangular = Triangular deriving (Eq)
+data Triangle = Triangle deriving (Eq)
 
 householderPart ::
    (Shape.C height, Shape.C width) =>
    Householder height width ->
-   (Shape.Index height, Shape.Index width) -> Either Reflector Triangular
+   (Shape.Index height, Shape.Index width) -> Either Reflector Triangle
 householderPart (Householder _ height width) (r,c) =
    if Shape.offset height r > Shape.offset width c
      then Left Reflector
-     else Right Triangular
+     else Right Triangle
 
 instance
    (Shape.C height, Shape.C width) =>
       Shape.C (Householder height width) where
 
    type Index (Householder height width) =
-            (Either Reflector Triangular,
+            (Either Reflector Triangle,
              (Shape.Index height, Shape.Index width))
 
    indices sh@(Householder _ height width) =
@@ -125,3 +181,147 @@
       Shape.inBounds (height,width) ix
       &&
       part == householderPart sh ix
+
+
+{- |
+Store the upper triangular half of a real symmetric or complex Hermitian matrix.
+-}
+data Hermitian size =
+   Hermitian {
+      hermitianOrder :: Order,
+      hermitianSize :: size
+   } deriving (Eq, Show)
+
+type instance HeightOf (Hermitian size) = size
+type instance WidthOf (Hermitian size) = size
+
+uploFromOrder :: Order -> Char
+uploFromOrder RowMajor = 'L'
+uploFromOrder ColumnMajor = 'U'
+
+instance (Shape.C size) => Shape.C (Hermitian size) where
+   type Index (Hermitian size) = (Shape.Index size, Shape.Index size)
+
+   indices (Hermitian _ size) =
+      let ixs = Shape.indices size
+      in  concat $ zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs
+
+   uncheckedOffset sh ix =
+      snd $ Shape.uncheckedSizeOffset sh ix
+
+   sizeOffset sh ix =
+      if Shape.inBounds sh ix
+        then Shape.uncheckedSizeOffset sh ix
+        else error "Shape.Hermitian.sizeOffset: wrong matrix part"
+
+   uncheckedSizeOffset (Hermitian RowMajor size) (rs,cs) =
+      let (s,r) = Shape.uncheckedSizeOffset size rs
+          c = Shape.uncheckedOffset size cs
+      in  (s, triangleSize s - triangleSize (s-r) + c-r)
+   uncheckedSizeOffset (Hermitian ColumnMajor size) (rs,cs) =
+      let (s,r) = Shape.uncheckedSizeOffset size rs
+          c = Shape.uncheckedOffset size cs
+      in  (s, triangleSize c + r)
+
+   size (Hermitian _ size) = triangleSize $ Shape.size size
+   uncheckedSize (Hermitian _ size) = triangleSize $ Shape.uncheckedSize size
+   inBounds (Hermitian _ size) ix@(r,c) =
+      Shape.inBounds (size,size) ix
+      &&
+      Shape.offset size r <= Shape.offset size c
+
+
+data Triangular uplo size =
+   Triangular {
+      triangularUplo :: uplo,
+      triangularOrder :: Order,
+      triangularSize :: size
+   } deriving (Eq, Show)
+
+type instance HeightOf (Triangular uplo size) = size
+type instance WidthOf (Triangular uplo size) = size
+
+data Lower = Lower deriving (Eq, Show)
+data Upper = Upper deriving (Eq, Show)
+
+type LowerTriangular = Triangular Lower
+type UpperTriangular = Triangular Upper
+
+triangularTransposeUp :: LowerTriangular sh -> UpperTriangular sh
+triangularTransposeUp (Triangular Lower order size) =
+   Triangular Upper (flipOrder order) size
+
+triangularTransposeDown :: UpperTriangular sh -> LowerTriangular sh
+triangularTransposeDown (Triangular Upper order size) =
+   Triangular Lower (flipOrder order) size
+
+
+class Uplo uplo where switchUplo :: f Lower -> f Upper -> f uplo
+instance Uplo Lower where switchUplo f _ = f
+instance Uplo Upper where switchUplo _ f = f
+
+autoUplo :: Uplo uplo => uplo
+autoUplo = runIdentity $ switchUplo (Identity Lower) (Identity Upper)
+
+uploOrder :: Uplo uplo => uplo -> Order -> Order
+uploOrder uplo order = caseUplo uplo (flipOrder order) order
+
+getUploConst :: uplo -> Const a uplo -> a
+getUploConst _ = getConst
+
+caseUplo :: Uplo uplo => uplo -> a -> a -> a
+caseUplo uplo lo up =
+   getUploConst uplo $ switchUplo (Const lo) (Const up)
+
+instance (Uplo uplo, Shape.C size) => Shape.C (Triangular uplo size) where
+   type Index (Triangular uplo size) = (Shape.Index size, Shape.Index size)
+
+   indices (Triangular uplo _ size) =
+      let ixs = Shape.indices size
+          rcs = concat $ zipWith (\r cs -> map ((,) r) cs) ixs $ tails ixs
+      in  caseUplo uplo (map swap rcs) rcs
+
+   uncheckedOffset sh ix =
+      snd $ Shape.uncheckedSizeOffset sh ix
+
+   sizeOffset sh ix =
+      if Shape.inBounds sh ix
+        then Shape.uncheckedSizeOffset sh ix
+        else error "Shape.Triangular.sizeOffset: wrong matrix part"
+
+   uncheckedSizeOffset (Triangular uplo RowMajor size) (rs,cs) =
+      let (s,r) = Shape.uncheckedSizeOffset size rs
+          c = Shape.uncheckedOffset size cs
+      in  (s,
+           caseUplo uplo
+               (triangleSize r + c)
+               (triangleSize s - triangleSize (s-r) + c-r))
+   uncheckedSizeOffset (Triangular uplo ColumnMajor size) (rs,cs) =
+      let (s,r) = Shape.uncheckedSizeOffset size rs
+          c = Shape.uncheckedOffset size cs
+      in  (s,
+           caseUplo uplo
+               (triangleSize s - triangleSize (s-c) + r-c)
+               (triangleSize c + r))
+
+   size (Triangular _ _ size) = triangleSize $ Shape.size size
+   uncheckedSize (Triangular _ _ size) = triangleSize $ Shape.uncheckedSize size
+   inBounds (Triangular uplo _ size) ix@(r,c) =
+      Shape.inBounds (size,size) ix
+      &&
+      caseUplo uplo
+         (Shape.offset size r >= Shape.offset size c)
+         (Shape.offset size r <= Shape.offset size c)
+
+triangleSize :: Int -> Int
+triangleSize n = div (n*(n+1)) 2
+
+triangleRoot :: Floating a => a -> a
+triangleRoot size = (sqrt (8*size+1)-1)/2
+
+triangleExtent :: String -> Int -> Int
+triangleExtent name size =
+   let n = round (triangleRoot (fromIntegral size :: Double))
+   in if size == triangleSize n
+        then n
+        else error (name ++ ": no triangular number of elements")
diff --git a/src/Numeric/LAPACK/Matrix/Square.hs b/src/Numeric/LAPACK/Matrix/Square.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Square.hs
@@ -0,0 +1,202 @@
+module Numeric.LAPACK.Matrix.Square (
+   Square,
+   size,
+   toGeneral,
+   fromGeneral,
+   fromScalar,
+   toScalar,
+   fromList,
+   autoFromList,
+
+   transpose,
+   adjoint,
+
+   identity,
+   identityFrom,
+   diagonal,
+   getDiagonal,
+   trace,
+
+   multiply,
+   square,
+   power,
+   ) where
+
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Vector as Vector
+import qualified Numeric.LAPACK.Private as Private
+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))
+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private (zero, one)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.BLAS.FFI.Generic as BlasGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Storable (Storable, peek, poke)
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+import Data.Function.HT (powerAssociative)
+
+
+type Square sh = Array (MatrixShape.Square sh)
+
+size :: Square sh a -> sh
+size = MatrixShape.squareSize . Array.shape
+
+toGeneral :: Square sh a -> General sh sh a
+toGeneral (Array sh a) = Array (MatrixShape.generalFromSquare sh) a
+
+fromGeneral :: (Eq sh) => General sh sh a -> Square sh a
+fromGeneral (Array (MatrixShape.General order height width) a) =
+   if height==width
+     then Array (MatrixShape.Square order height) a
+     else error "Square.fromGeneral: no square shape"
+
+
+fromScalar :: (Storable a) => a -> Square () a
+fromScalar a =
+   Array.unsafeCreate (MatrixShape.Square RowMajor ()) $ flip poke a
+
+toScalar :: (Storable a) => Square () a -> a
+toScalar (Array (MatrixShape.Square _ ()) a) =
+   unsafePerformIO $ withForeignPtr a peek
+
+fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Square sh a
+fromList sh =
+   Array.fromList (MatrixShape.Square RowMajor sh)
+
+autoFromList :: (Storable a) => [a] -> Square ZeroInt a
+autoFromList xs =
+   let n = length xs
+       m = round $ sqrt (fromIntegral n :: Double)
+   in if n == m*m
+        then fromList (zeroInt m) xs
+        else error "Square.autoFromList: no quadratic number of elements"
+
+
+transpose :: Square sh a -> Square sh a
+transpose = Array.mapShape MatrixShape.transposeSquare
+
+{- |
+conjugate transpose
+-}
+adjoint :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
+adjoint = transpose . Vector.conjugate
+
+
+identity :: (Shape.C sh, Class.Floating a) => sh -> Square sh a
+identity = identityOrder ColumnMajor
+
+identityFrom :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
+identityFrom (Array (MatrixShape.Square order sh) _) = identityOrder order sh
+
+identityOrder, _identityOrder ::
+   (Shape.C sh, Class.Floating a) => Order -> sh -> Square sh a
+identityOrder order sh =
+   Array.unsafeCreate (MatrixShape.Square order sh) $ \aPtr ->
+   evalContT $ do
+      uploPtr <- Call.char 'A'
+      nPtr <- Call.cint $ Shape.size sh
+      alphaPtr <- Call.number zero
+      betaPtr <- Call.number one
+      liftIO $ LapackGen.laset uploPtr nPtr nPtr alphaPtr betaPtr aPtr nPtr
+
+_identityOrder order sh =
+   Array.unsafeCreateWithSize (MatrixShape.Square order sh) $ \blockSize yPtr ->
+   evalContT $ do
+      nPtr <- Call.alloca
+      xPtr <- Call.number zero
+      incxPtr <- Call.cint 0
+      incyPtr <- Call.cint 1
+      liftIO $ do
+         poke nPtr $ fromIntegral blockSize
+         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+         let n = fromIntegral $ Shape.size sh
+         poke nPtr n
+         poke xPtr one
+         poke incyPtr (n+1)
+         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+
+diagonal :: (Shape.C sh, Class.Floating a) => Vector sh a -> Square sh a
+diagonal (Array sh x) =
+   Array.unsafeCreateWithSize (MatrixShape.Square ColumnMajor sh) $
+      \blockSize yPtr ->
+   evalContT $ do
+      nPtr <- Call.alloca
+      xPtr <- ContT $ withForeignPtr x
+      zPtr <- Call.number zero
+      incxPtr <- Call.cint 1
+      incyPtr <- Call.cint 1
+      inczPtr <- Call.cint 0
+      liftIO $ do
+         poke nPtr $ fromIntegral blockSize
+         BlasGen.copy nPtr zPtr inczPtr yPtr incyPtr
+         let n = fromIntegral $ Shape.size sh
+         poke nPtr n
+         poke incyPtr (n+1)
+         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+
+getDiagonal :: (Shape.C sh, Class.Floating a) => Square sh a -> Vector sh a
+getDiagonal (Array (MatrixShape.Square _ sh) x) =
+   Array.unsafeCreateWithSize sh $ \n yPtr -> evalContT $ do
+      nPtr <- Call.cint n
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint (n+1)
+      incyPtr <- Call.cint 1
+      liftIO $ BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
+
+trace :: (Shape.C sh, Class.Floating a) => Square sh a -> a
+trace (Array (MatrixShape.Square _ sh) x) = unsafePerformIO $ do
+   let n = Shape.size sh
+   withForeignPtr x $ \xPtr -> Private.sum n xPtr (n+1)
+
+
+multiply ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   Square sh a -> Square sh a -> Square sh a
+multiply
+   (Array (MatrixShape.Square orderA shA) a)
+   (Array (MatrixShape.Square orderB shB) b) =
+      Array.unsafeCreate (MatrixShape.Square ColumnMajor shA) $ \cPtr -> do
+   Call.assert "Square.multiply: shapes mismatch" (shA == shB)
+   let n = Shape.size shA
+   Private.multiplyMatrix orderA orderB n n n a b cPtr
+
+square :: (Shape.C sh, Class.Floating a) => Square sh a -> Square sh a
+square a = multiplyCommutativeUnchecked a a
+
+power ::
+   (Shape.C sh, Class.Floating a) =>
+   Integer -> Square sh a -> Square sh a
+power n a =
+   powerAssociative multiplyCommutativeUnchecked (identityFrom a) a n
+
+{-
+orderA and orderB must be equal but this is not checked.
+-}
+multiplyCommutativeUnchecked ::
+   (Shape.C sh, Class.Floating a) =>
+   Square sh a -> Square sh a -> Square sh a
+multiplyCommutativeUnchecked
+   (Array shape@(MatrixShape.Square  order  sh) a)
+   (Array       (MatrixShape.Square _order _sh) b) =
+      Array.unsafeCreate shape $ \cPtr ->
+   let n = Shape.size sh
+       (at,bt) =
+         case order of
+            ColumnMajor -> (a,b)
+            RowMajor -> (b,a)
+   in  Private.multiplyMatrix ColumnMajor ColumnMajor n n n at bt cPtr
diff --git a/src/Numeric/LAPACK/Matrix/Triangular.hs b/src/Numeric/LAPACK/Matrix/Triangular.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Triangular.hs
@@ -0,0 +1,315 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Matrix.Triangular (
+   Triangular, MatrixShape.Uplo(..),
+   Upper, Lower,
+   fromList, autoFromList,
+   lowerFromList, autoLowerFromList,
+   upperFromList, autoUpperFromList,
+   identity,
+   diagonal,
+   getDiagonal,
+   transposeUp, transposeDown,
+   adjointUp, adjointDown,
+
+   toSquare,
+
+   multiplyVectorLeft,
+   multiplyVectorRight,
+   square,
+   multiply,
+   multiplySquareLeft,
+   multiplyGeneralLeft,
+   multiplySquareRight,
+   multiplyGeneralRight,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix.Triangular.Private
+         (diagonalPointers, pack, unpack, unpackZero, unpackToTemp)
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor,ColumnMajor),
+          flipOrder, transposeFromOrder, uploFromOrder, uploOrder)
+import Numeric.LAPACK.Matrix.Square (Square)
+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private (fill, zero, one, copyBlock)
+
+import qualified Numeric.BLAS.FFI.Generic as BlasGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (Storable, poke, peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+import Data.Foldable (forM_)
+
+
+type Triangular uplo sh = Array (MatrixShape.Triangular uplo sh)
+
+type Lower sh = Array (MatrixShape.LowerTriangular sh)
+type Upper sh = Array (MatrixShape.UpperTriangular sh)
+
+transposeUp :: Lower sh a -> Upper sh a
+transposeUp (Array sh a) =
+   Array (MatrixShape.triangularTransposeUp sh) a
+
+transposeDown :: Upper sh a -> Lower sh a
+transposeDown (Array sh a) =
+   Array (MatrixShape.triangularTransposeDown sh) a
+
+adjointUp :: (Shape.C sh, Class.Floating a) => Lower sh a -> Upper sh a
+adjointUp = Vector.conjugate . transposeUp
+
+adjointDown :: (Shape.C sh, Class.Floating a) => Upper sh a -> Lower sh a
+adjointDown = Vector.conjugate . transposeDown
+
+
+fromList ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Storable a) =>
+   Order -> sh -> [a] -> Triangular uplo sh a
+fromList order sh =
+   Array.fromList (MatrixShape.Triangular MatrixShape.autoUplo order sh)
+
+lowerFromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Lower sh a
+lowerFromList = fromList
+
+upperFromList :: (Shape.C sh, Storable a) => Order -> sh -> [a] -> Upper sh a
+upperFromList = fromList
+
+
+autoFromList ::
+   (MatrixShape.Uplo uplo, Storable a) =>
+   Order -> [a] -> Triangular uplo ZeroInt a
+autoFromList order xs =
+   fromList order
+      (zeroInt $ MatrixShape.triangleExtent "Triangular.autoFromList" $
+       length xs)
+      xs
+
+autoLowerFromList :: (Storable a) => Order -> [a] -> Lower ZeroInt a
+autoLowerFromList = autoFromList
+
+autoUpperFromList :: (Storable a) => Order -> [a] -> Upper ZeroInt a
+autoUpperFromList = autoFromList
+
+
+toSquare ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Triangular uplo sh a -> Square sh a
+toSquare (Array (MatrixShape.Triangular uplo order sh) a) =
+   Array.unsafeCreate (MatrixShape.Square order sh) $ \bPtr ->
+      withForeignPtr a $ \aPtr ->
+         unpackZero (uploOrder uplo order) (Shape.size sh) aPtr bPtr
+
+
+identity ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Order -> sh -> Triangular uplo sh a
+identity order sh =
+   let (realOrder, uplo) = autoUploOrder order
+   in Array.unsafeCreate (MatrixShape.Triangular uplo order sh) $ \aPtr -> do
+      let n = Shape.size sh
+      fill zero (MatrixShape.triangleSize n) aPtr
+      forM_ (diagonalPointers realOrder n aPtr aPtr) $ flip poke one . snd
+
+diagonal ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Order -> Vector sh a -> Triangular uplo sh a
+diagonal order (Array sh x) =
+   let (realOrder, uplo) = autoUploOrder order
+   in Array.unsafeCreate (MatrixShape.Triangular uplo order sh) $ \aPtr -> do
+      let n = Shape.size sh
+      fill zero (MatrixShape.triangleSize n) aPtr
+      withForeignPtr x $ \xPtr ->
+         forM_ (diagonalPointers realOrder n xPtr aPtr) $
+            \(srcPtr,dstPtr) -> poke dstPtr =<< peek srcPtr
+
+getDiagonal ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Class.Floating a) =>
+   Triangular uplo sh a -> Vector sh a
+getDiagonal (Array (MatrixShape.Triangular uplo order sh) a) =
+      Array.unsafeCreate sh $ \xPtr -> do
+   withForeignPtr a $ \aPtr ->
+      mapM_
+         (\(dstPtr,srcPtr) -> poke dstPtr =<< peek srcPtr)
+         (diagonalPointers (uploOrder uplo order) (Shape.size sh) xPtr aPtr)
+
+
+multiplyVectorLeft, multiplyVectorRight ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>
+   Triangular uplo sh a -> Vector sh a -> Vector sh a
+multiplyVectorLeft = multiplyVector True
+multiplyVectorRight = multiplyVector False
+
+multiplyVector ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>
+   Bool -> Triangular uplo sh a -> Vector sh a -> Vector sh a
+multiplyVector transp
+   (Array (MatrixShape.Triangular uplo order shA) a) (Array shX x) =
+      Array.unsafeCreate shX $ \yPtr -> do
+   Call.assert "Triangular.multiplyVector: width shapes mismatch" (shA == shX)
+   let n = Shape.size shA
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder $ uploOrder uplo order
+      transPtr <-
+         Call.char $ transposeFromOrder $
+         (if transp then flipOrder else id) order
+      diagPtr <- Call.char 'N'
+      nPtr <- Call.cint n
+      aPtr <- ContT $ withForeignPtr a
+      xPtr <- ContT $ withForeignPtr x
+      incyPtr <- Call.cint 1
+      liftIO $ do
+         copyBlock n xPtr yPtr
+         BlasGen.tpmv uploPtr transPtr diagPtr nPtr aPtr yPtr incyPtr
+
+
+square ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>
+   Triangular uplo sh a -> Triangular uplo sh a
+square
+   (Array shape@(MatrixShape.Triangular uplo order sh) a) =
+      Array.unsafeCreate shape $ \bpPtr -> do
+   let n = Shape.size sh
+   evalContT $ do
+      sidePtr <- Call.char 'L'
+      let realOrder = uploOrder uplo order
+      uploPtr <- Call.char $ uploFromOrder realOrder
+      transPtr <- Call.char 'N'
+      diagPtr <- Call.char 'N'
+      nPtr <- Call.cint n
+      let ldPtr = nPtr
+      aPtr <- unpackToTemp (unpack realOrder) n a
+      bPtr <- unpackToTemp (unpackZero realOrder) n a
+      alphaPtr <- Call.number one
+      liftIO $ do
+         BlasGen.trmm sidePtr uploPtr transPtr diagPtr
+            nPtr nPtr alphaPtr aPtr ldPtr bPtr ldPtr
+         pack realOrder n bPtr bpPtr
+
+multiply ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>
+   Triangular uplo sh a -> Triangular uplo sh a -> Triangular uplo sh a
+multiply
+   (Array        (MatrixShape.Triangular uploA orderA shA) a)
+   (Array shapeB@(MatrixShape.Triangular uploB orderB shB) b) =
+      Array.unsafeCreate shapeB $ \cpPtr -> do
+   Call.assert "Triangular.multiply: width shapes mismatch" (shA == shB)
+   let n = Shape.size shA
+   evalContT $ do
+      let (side,trans) =
+            case orderB of
+               ColumnMajor -> ('L', orderA)
+               RowMajor -> ('R', flipOrder orderA)
+      sidePtr <- Call.char side
+      let realOrderA = uploOrder uploA orderA
+      let realOrderB = uploOrder uploB orderB
+      uploPtr <- Call.char $ uploFromOrder realOrderA
+      transPtr <- Call.char $ transposeFromOrder trans
+      diagPtr <- Call.char 'N'
+      nPtr <- Call.cint n
+      let ldPtr = nPtr
+      aPtr <- unpackToTemp (unpack realOrderA) n a
+      bPtr <- unpackToTemp (unpackZero realOrderB) n b
+      alphaPtr <- Call.number one
+      liftIO $ do
+         BlasGen.trmm sidePtr uploPtr transPtr diagPtr
+            nPtr nPtr alphaPtr aPtr ldPtr bPtr ldPtr
+         pack realOrderB n bPtr cpPtr
+
+
+multiplySquareLeft ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>
+   Square sh a -> Triangular uplo sh a -> Square sh a
+multiplySquareLeft
+   (Array shapeB@(MatrixShape.Square orderB shB) b)
+   (Array        (MatrixShape.Triangular uploA orderA shA) a) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Triangular.multiplySquareLeft: shapes mismatch" (shA == shB)
+   let n = Shape.size shB
+   MatrixShape.caseUplo uploA
+      (multiplyAux MatrixShape.Upper)
+      (multiplyAux MatrixShape.Lower)
+      (flipOrder orderA) n a (flipOrder orderB) n b cPtr
+
+multiplyGeneralLeft ::
+   (MatrixShape.Uplo uplo,
+    Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a -> Triangular uplo width a -> General height width a
+multiplyGeneralLeft
+   (Array shapeB@(MatrixShape.General orderB height width) b)
+   (Array        (MatrixShape.Triangular uploA orderA shA) a) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Triangular.multiplyGeneralLeft: shapes mismatch" (shA == width)
+   MatrixShape.caseUplo uploA
+      (multiplyAux MatrixShape.Upper)
+      (multiplyAux MatrixShape.Lower)
+      (flipOrder orderA) (Shape.size width) a
+      (flipOrder orderB) (Shape.size height) b cPtr
+
+multiplySquareRight ::
+   (MatrixShape.Uplo uplo, Shape.C sh, Eq sh, Class.Floating a) =>
+   Triangular uplo sh a -> Square sh a -> Square sh a
+multiplySquareRight
+   (Array        (MatrixShape.Triangular uploA orderA shA) a)
+   (Array shapeB@(MatrixShape.Square orderB shB) b) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Triangular.multiplySquareRight: shapes mismatch" (shA == shB)
+   let n = Shape.size shB
+   multiplyAux uploA orderA n a orderB n b cPtr
+
+multiplyGeneralRight ::
+   (MatrixShape.Uplo uplo,
+    Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
+   Triangular uplo height a -> General height width a -> General height width a
+multiplyGeneralRight
+   (Array        (MatrixShape.Triangular uploA orderA shA) a)
+   (Array shapeB@(MatrixShape.General orderB height width) b) =
+      Array.unsafeCreate shapeB $ \cPtr -> do
+   Call.assert "Triangular.multiplyGeneralRight: shapes mismatch"
+      (shA == height)
+   multiplyAux
+      uploA orderA (Shape.size height) a orderB (Shape.size width) b cPtr
+
+multiplyAux ::
+   (MatrixShape.Uplo uplo, Class.Floating a) =>
+   uplo ->
+   Order -> Int -> ForeignPtr a ->
+   Order -> Int -> ForeignPtr a -> Ptr a -> IO ()
+multiplyAux uploA orderA m0 a orderB n0 b cPtr =
+   evalContT $ do
+      let (side,trans,(m,n)) =
+            case orderB of
+               ColumnMajor -> ('L', orderA, (m0,n0))
+               RowMajor -> ('R', flipOrder orderA, (n0,m0))
+      sidePtr <- Call.char side
+      let realOrderA = uploOrder uploA orderA
+      uploPtr <- Call.char $ uploFromOrder realOrderA
+      transPtr <- Call.char $ transposeFromOrder trans
+      diagPtr <- Call.char 'N'
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number one
+      aPtr <- unpackToTemp (unpack realOrderA) m0 a
+      ldaPtr <- Call.cint m0
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint m
+      liftIO $ do
+         copyBlock (m0*n0) bPtr cPtr
+         BlasGen.trmm sidePtr uploPtr transPtr diagPtr
+            mPtr nPtr alphaPtr aPtr ldaPtr cPtr ldbPtr
+
+
+autoUploOrder :: MatrixShape.Uplo uplo => Order -> (Order, uplo)
+autoUploOrder order =
+   case MatrixShape.autoUplo of
+      uplo -> (uploOrder uplo order, uplo)
diff --git a/src/Numeric/LAPACK/Matrix/Triangular/Private.hs b/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Triangular/Private.hs
@@ -0,0 +1,121 @@
+module Numeric.LAPACK.Matrix.Triangular.Private where
+
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder, triangleSize)
+import Numeric.LAPACK.Private (pointerSeq, copyToTemp, lacgv, fill, zero)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.Marshal.Array (advancePtr)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (Storable, poke)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+import Data.Foldable (forM_)
+
+
+diagonalPointers ::
+   (Storable a, Storable ar) =>
+   Order -> Int -> Ptr ar -> Ptr a -> [(Ptr ar, Ptr a)]
+diagonalPointers order n xPtr aPtr =
+   take n $ zip (pointerSeq 1 xPtr) $ scanl advancePtr aPtr $
+   case order of
+      RowMajor -> iterate pred n
+      ColumnMajor -> iterate succ 2
+
+
+columnMajorPointers ::
+   (Storable a) => Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
+columnMajorPointers n fullPtr packedPtr =
+   let ds = iterate succ 1
+   in  take n $ zip ds $
+       zip
+         (zip (pointerSeq 1 fullPtr) (pointerSeq n fullPtr))
+         (scanl advancePtr packedPtr ds)
+
+rowMajorPointers ::
+   (Storable a) => Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
+rowMajorPointers n fullPtr packedPtr =
+   let ds = iterate pred n
+   in  take n $ zip ds $
+       zip (pointerSeq (n+1) fullPtr) (scanl advancePtr packedPtr ds)
+
+
+forPointers :: [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()
+forPointers xs act =
+   alloca $ \nPtr ->
+   forM_ xs $ \(d,ptrs) -> do
+      poke nPtr $ fromIntegral d
+      act nPtr ptrs
+
+
+copyTriangleToTemp ::
+   Class.Floating a =>
+   Order -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
+copyTriangleToTemp order n a = do
+   let aSize = triangleSize n
+   apPtr <- copyToTemp aSize a
+   liftIO $ evalContT $ do
+      aSizePtr <- Call.cint aSize
+      incPtr <- Call.cint 1
+      case order of
+         RowMajor -> liftIO $ lacgv aSizePtr apPtr incPtr
+         ColumnMajor -> return ()
+   return apPtr
+
+
+unpackToTemp ::
+   Storable a =>
+   (Int -> Ptr a -> Ptr a -> IO ()) ->
+   Int -> ForeignPtr a -> ContT r IO (Ptr a)
+unpackToTemp f n a = do
+   apPtr <- ContT $ withForeignPtr a
+   aPtr <- Call.allocaArray (n*n)
+   liftIO $ f n apPtr aPtr
+   return aPtr
+
+
+unpack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
+unpack order n packedPtr fullPtr =
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder order
+      nPtr <- Call.cint n
+      ldaPtr <- Call.cint n
+      liftIO $ withInfo $ LapackGen.tpttr uploPtr nPtr packedPtr fullPtr ldaPtr
+
+pack :: Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
+pack order n fullPtr packedPtr =
+   evalContT $ do
+      uploPtr <- Call.char $ uploFromOrder order
+      nPtr <- Call.cint n
+      ldaPtr <- Call.cint n
+      liftIO $ withInfo $ LapackGen.trttp uploPtr nPtr fullPtr ldaPtr packedPtr
+
+
+unpackZero, _unpackZero ::
+   Class.Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
+_unpackZero order n packedPtr fullPtr = do
+   fill zero (n*n) fullPtr
+   unpack order n packedPtr fullPtr
+
+unpackZero order n packedPtr fullPtr = do
+   fillTriangle zero (flipOrder order) n fullPtr
+   unpack order n packedPtr fullPtr
+
+fillTriangle :: Class.Floating a => a -> Order -> Int -> Ptr a -> IO ()
+fillTriangle z order n aPtr = evalContT $ do
+   uploPtr <- Call.char $ uploFromOrder order
+   nPtr <- Call.cint n
+   zPtr <- Call.number z
+   liftIO $ LapackGen.laset uploPtr nPtr nPtr zPtr zPtr aPtr nPtr
+
+
+withInfo :: (Ptr CInt -> IO ()) -> IO ()
+withInfo = alloca
diff --git a/src/Numeric/LAPACK/Orthogonal.hs b/src/Numeric/LAPACK/Orthogonal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Orthogonal.hs
@@ -0,0 +1,372 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Orthogonal (
+   leastSquares,
+   minimumNorm,
+   leastSquaresMinimumNorm,
+   pseudoInverseRCond,
+
+   Householder,
+   householder,
+   householderDecompose,
+   householderDeterminant,
+   determinant,
+   householderExtractQ,
+   householderExtractR,
+   orthogonalComplement,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import Numeric.LAPACK.Matrix
+         (General, ZeroInt, zeroInt, transpose, identity, dropColumns)
+import Numeric.LAPACK.Matrix.Square (Square)
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Private as Private
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor, ColumnMajor), transposeFromOrder)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private
+         (RealOf, zero, fill,
+          copySubMatrix, copyBlock, copyToTemp,
+          copyToColumnMajor, copyToSubColumnMajor,
+          withAutoWorkspaceInfo, allocArray, allocHigherArray)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
+import qualified Numeric.LAPACK.FFI.Real as LapackReal
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Foreign.Marshal.Array (advancePtr)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (peek)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Data.Complex (Complex)
+import Data.Tuple.HT (mapSnd)
+
+
+{- |
+If @x = leastSquares a b@
+then @x@ minimizes @Vector.norm2 (multiply a x `sub` b)@.
+
+Precondition: @a@ must have full rank and @height a >= width a@.
+-}
+leastSquares ::
+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
+   General height width a -> General height nrhs a -> General width nrhs a
+leastSquares
+   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)
+   (Array        (MatrixShape.General orderB heightB widthB) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $
+         \xPtr -> do
+   Call.assert "leastSquares: height shapes mismatch" (heightA == heightB)
+   Call.assert "leastSquares: height of 'a' must be at least the width"
+      (Shape.size heightA >= Shape.size widthA)
+   let (m,n) = MatrixShape.dimensions shapeA
+   let lda = m
+   let nrhs = Shape.size widthB
+   let ldb = Shape.size heightB
+   let ldx = Shape.size widthA
+   evalContT $ do
+      transPtr <- Call.char $ transposeFromOrder orderA
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      aPtr <- copyToTemp (Shape.size shapeA) a
+      ldaPtr <- Call.cint lda
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint ldb
+      let bSize = Shape.size (heightB,widthB)
+      btmpPtr <- Call.allocaArray bSize
+      liftIO $ copyToColumnMajor orderB ldb nrhs bPtr btmpPtr
+      liftIO $ withAutoWorkspaceInfo "gels" $
+         LapackGen.gels transPtr
+            mPtr nPtr nrhsPtr aPtr ldaPtr btmpPtr ldbPtr
+      liftIO $ copySubMatrix ldx nrhs ldb btmpPtr ldx xPtr
+
+{- |
+The vector @x@ with @x = minimumNorm a b@
+is the vector with minimal @Vector.norm2 x@
+that satisfies @multiply a x == b@.
+
+Precondition: @a@ must have full rank and @height a <= width a@.
+-}
+minimumNorm ::
+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
+   General height width a -> General height nrhs a -> General width nrhs a
+minimumNorm
+   (Array shapeA@(MatrixShape.General orderA heightA widthA) a)
+   (Array        (MatrixShape.General orderB heightB widthB) b) =
+      Array.unsafeCreate (MatrixShape.General ColumnMajor widthA widthB) $
+         \xPtr -> do
+   Call.assert "minimumNorm: height shapes mismatch" (heightA == heightB)
+   Call.assert "minimumNorm: width of 'a' must be at least the height"
+      (Shape.size widthA >= Shape.size heightA)
+   let (m,n) = MatrixShape.dimensions shapeA
+   let lda = m
+   let nrhs = Shape.size widthB
+   let ldb = Shape.size heightB
+   let ldx = Shape.size widthA
+   evalContT $ do
+      transPtr <- Call.char $ transposeFromOrder orderA
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      aPtr <- copyToTemp (Shape.size shapeA) a
+      ldaPtr <- Call.cint lda
+      bPtr <- ContT $ withForeignPtr b
+      ldxPtr <- Call.cint ldx
+      liftIO $ copyToSubColumnMajor orderB ldb nrhs bPtr ldx xPtr
+      liftIO $ withAutoWorkspaceInfo "gels" $
+         LapackGen.gels transPtr
+            mPtr nPtr nrhsPtr aPtr ldaPtr xPtr ldxPtr
+
+{- |
+If @x = leastSquaresMinimumNorm a b@
+then @x@ is the vector with minimum @Vector.norm2 x@
+that minimizes @Vector.norm2 (multiply a x `sub` b)@.
+
+Matrix @a@ can have any rank
+but you must specify the reciprocal condition of the rank-truncated matrix.
+-}
+leastSquaresMinimumNorm ::
+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
+   RealOf a ->
+   General height width a -> General height nrhs a ->
+   (Int, General width nrhs a)
+leastSquaresMinimumNorm rcond
+   (Array (MatrixShape.General orderA heightA widthA) a)
+   (Array (MatrixShape.General orderB heightB widthB) b) =
+      unsafePerformIO $ do
+   Call.assert "leastSquaresMinimumNorm: height shapes mismatch"
+      (heightA == heightB)
+   let shapeX = MatrixShape.General ColumnMajor widthA widthB
+   let m = Shape.size heightA
+   let n = Shape.size widthA
+   let nrhs = Shape.size widthB
+   let aSize = m*n
+   let lda = m
+   evalContT $ do
+      aPtr <- ContT $ withForeignPtr a
+      atmpPtr <- Call.allocaArray aSize
+      liftIO $ copyToColumnMajor orderA m n aPtr atmpPtr
+      ldaPtr <- Call.cint lda
+      (x,(tmpPtr,ldtmp)) <- allocHigherArray shapeX m n nrhs
+      ldtmpPtr <- Call.cint ldtmp
+      bPtr <- ContT $ withForeignPtr b
+      liftIO $ copyToSubColumnMajor orderB m nrhs bPtr ldtmp tmpPtr
+      jpvtPtr <- Call.allocaArray n
+      rankPtr <- Call.alloca
+      gelsy m n nrhs atmpPtr ldaPtr tmpPtr ldtmpPtr jpvtPtr rcond rankPtr
+      rank <- liftIO $ fromIntegral <$> peek rankPtr
+      return (rank, x)
+
+
+type GELSY_ r ar a =
+   Int -> Int -> Int -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
+   Ptr CInt -> ar -> Ptr CInt -> ContT r IO ()
+
+newtype GELSY r a = GELSY {getGELSY :: GELSY_ r (RealOf a) a}
+
+gelsy :: (Class.Floating a) => GELSY_ r (RealOf a) a
+gelsy =
+   getGELSY $
+   Class.switchFloating
+      (GELSY gelsyReal)
+      (GELSY gelsyReal)
+      (GELSY gelsyComplex)
+      (GELSY gelsyComplex)
+
+gelsyReal :: (Class.Real a) => GELSY_ r a a
+gelsyReal m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do
+   mPtr <- Call.cint m
+   nPtr <- Call.cint n
+   nrhsPtr <- Call.cint nrhs
+   rcondPtr <- Call.real rcond
+   liftIO $ withAutoWorkspaceInfo "gelsy" $
+      LapackReal.gelsy mPtr nPtr nrhsPtr
+         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr
+
+gelsyComplex :: (Class.Real a) => GELSY_ r a (Complex a)
+gelsyComplex m n nrhs aPtr ldaPtr bPtr ldbPtr jpvtPtr rcond rankPtr = do
+   mPtr <- Call.cint m
+   nPtr <- Call.cint n
+   nrhsPtr <- Call.cint nrhs
+   rcondPtr <- Call.real rcond
+   rworkPtr <- Call.allocaArray (2*n)
+   liftIO $ withAutoWorkspaceInfo "gelsy" $ \workPtr lworkPtr infoPtr ->
+      LapackComplex.gelsy mPtr nPtr nrhsPtr
+         aPtr ldaPtr bPtr ldbPtr jpvtPtr rcondPtr rankPtr
+         workPtr lworkPtr rworkPtr infoPtr
+
+
+pseudoInverseRCond ::
+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   RealOf a -> General height width a -> (Int, General width height a)
+pseudoInverseRCond rcond a =
+   let (MatrixShape.General _ height width) = Array.shape a
+   in if Shape.size height < Shape.size width
+         then leastSquaresMinimumNorm rcond a $ identity height
+         else mapSnd transpose $
+              leastSquaresMinimumNorm rcond (transpose a) $
+              identity width
+
+
+type Householder height width = Array (MatrixShape.Householder height width)
+
+{-
+@(q,r) = householder a@
+means that @q@ is unitary and @r@ is upper triangular and @a = multiply q r@.
+-}
+householder ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a ->
+   (Square height a, General height width a)
+householder a =
+   let hh = householderDecompose a
+   in  (householderExtractQ hh, householderExtractR $ snd hh)
+
+householderDecompose ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> (Vector width a, Householder height width a)
+householderDecompose (Array shape@(MatrixShape.General order height width) a) =
+   unsafePerformIO $ do
+
+   let (m,n) = MatrixShape.dimensions shape
+   let lda = m
+   let mn = min m n
+   evalContT $ do
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      aPtr <- ContT $ withForeignPtr a
+      ldaPtr <- Call.cint lda
+      (qr,qrPtr) <- allocArray $ MatrixShape.Householder order height width
+      liftIO $ copyBlock (m*n) aPtr qrPtr
+      (tau,tauPtr) <- allocArray width
+      liftIO $ fill zero (n-mn) (advancePtr tauPtr mn)
+      liftIO $
+         case order of
+            RowMajor ->
+               withAutoWorkspaceInfo "gelqf" $
+                  LapackGen.gelqf mPtr nPtr qrPtr ldaPtr tauPtr
+            ColumnMajor ->
+               withAutoWorkspaceInfo "geqrf" $
+                  LapackGen.geqrf mPtr nPtr qrPtr ldaPtr tauPtr
+      return (tau, qr)
+
+householderDeterminant ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   Householder height width a -> a
+householderDeterminant
+      (Array (MatrixShape.Householder order height width) a) =
+   let m = Shape.size height
+       n = Shape.size width
+       k = case order of RowMajor -> n; ColumnMajor -> m
+   in unsafePerformIO $
+      withForeignPtr a $ \aPtr ->
+      Private.product (min m n) aPtr (k+1)
+
+
+{-|
+Generalized determinant - works also for non-square matrices.
+In contrast to the square root of the Gramian determinant
+it has the proper sign.
+-}
+determinant ::
+   (Shape.C height, Shape.C width, Eq a, Class.Floating a) =>
+   General height width a -> a
+determinant a =
+   let (tau,hh) = householderDecompose a
+   in  foldl (\x _ -> negate x)
+         (householderDeterminant hh)
+         (takeWhile (/=zero) $ Array.toList tau)
+
+
+householderExtractQ ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   (Vector width a, Householder height width a) -> Square height a
+householderExtractQ
+   (Array widthTau tau,
+    Array (MatrixShape.Householder order height width) qr) =
+
+   Array.unsafeCreate (MatrixShape.Square order height) $ \qPtr -> do
+
+   Call.assert "householderExtractQ: width shapes mismatch" (widthTau == width)
+
+   let m = Shape.size height
+   let k = min m $ Shape.size width
+   let lda = m
+   evalContT $ do
+      mPtr <- Call.cint m
+      kPtr <- Call.cint k
+      qrPtr <- ContT $ withForeignPtr qr
+      ldaPtr <- Call.cint lda
+      tauPtr <- ContT $ withForeignPtr tau
+      liftIO $
+         case order of
+            RowMajor -> do
+               copySubMatrix k m k qrPtr lda qPtr
+               withAutoWorkspaceInfo "unglq" $
+                  LapackGen.unglq mPtr mPtr kPtr qPtr ldaPtr tauPtr
+            ColumnMajor -> do
+               copyBlock (m*k) qrPtr qPtr
+               withAutoWorkspaceInfo "ungqr" $
+                  LapackGen.ungqr mPtr mPtr kPtr qPtr ldaPtr tauPtr
+
+householderExtractR ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   Householder height width a -> General height width a
+householderExtractR
+      (Array (MatrixShape.Householder order height width) qr) =
+
+   Array.unsafeCreate (MatrixShape.General order height width) $
+      \rPtr -> do
+
+   let (uplo, (m,n)) =
+         case order of
+            RowMajor -> ('L', (Shape.size width, Shape.size height))
+            ColumnMajor -> ('U', (Shape.size height, Shape.size width))
+   fill zero (m*n) rPtr
+   evalContT $ do
+      uploPtr <- Call.char uplo
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      qrPtr <- ContT $ withForeignPtr qr
+      ldqrPtr <- Call.cint m
+      ldrPtr <- Call.cint m
+      liftIO $ LapackGen.lacpy uploPtr mPtr nPtr qrPtr ldqrPtr rPtr ldrPtr
+
+{- |
+For an m-by-n-matrix @a@ with m>=n
+this function computes an m-by-(m-n)-matrix @b@
+such that @Matrix.multiply (transpose b) a@ is a zero matrix.
+The function does not try to compensate a rank deficiency of @a@.
+That is, @a|||b@ has full rank if and only if @a@ has full rank.
+
+For full-rank matrices you might also call this @kernel@ or @nullspace@.
+-}
+orthogonalComplement ::
+   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
+   General height width a -> General height ZeroInt a
+orthogonalComplement a =
+   dropColumns (Shape.size $ MatrixShape.generalWidth $ Array.shape a) $
+   Array.mapShape zeroIntWidth $
+   Square.toGeneral $ householderExtractQ $ householderDecompose a
+
+zeroIntWidth ::
+   (Shape.C width) =>
+   MatrixShape.General height width -> MatrixShape.General height ZeroInt
+zeroIntWidth (MatrixShape.General order height width) =
+   MatrixShape.General order height (zeroInt $ Shape.size width)
diff --git a/src/Numeric/LAPACK/Private.hs b/src/Numeric/LAPACK/Private.hs
--- a/src/Numeric/LAPACK/Private.hs
+++ b/src/Numeric/LAPACK/Private.hs
@@ -1,23 +1,38 @@
 {-# LANGUAGE TypeFamilies #-}
 module Numeric.LAPACK.Private where
 
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor, ColumnMajor), transposeFromOrder)
+
 import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
 import qualified Numeric.BLAS.FFI.Real as BlasReal
 import qualified Numeric.BLAS.FFI.Generic as BlasGen
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
 import Foreign.Marshal.Array (advancePtr)
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt)
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, mallocForeignPtrArray)
 import Foreign.Ptr (Ptr)
-import Foreign.Storable (Storable, peek)
+import Foreign.Storable (Storable, poke, peek)
 
-import Control.Monad.Trans.Cont (evalContT)
+import Text.Printf (printf)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
 import Control.Monad.IO.Class (liftIO)
+import Control.Monad (foldM)
+import Control.Applicative ((<$>))
 
 import Data.Functor.Identity (Identity(Identity, runIdentity))
 
-import Data.Complex (Complex)
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
 
+import qualified Data.Complex as Complex
+import Data.Complex (Complex((:+)))
+
 import Prelude hiding (sum)
 
 
@@ -25,10 +40,12 @@
 
 type instance RealOf Float = Float
 type instance RealOf Double = Double
-type instance RealOf (Complex Float) = Float
-type instance RealOf (Complex Double) = Double
+type instance RealOf (Complex a) = a
 
 
+type ComplexOf x = Complex (RealOf x)
+
+
 zero, one, minusOne :: Class.Floating a => a
 zero =
    runIdentity $
@@ -41,11 +58,8 @@
    Class.switchFloating
       (Identity (-1)) (Identity (-1)) (Identity (-1)) (Identity (-1))
 
-oneReal :: Class.Real a => a
-oneReal = runIdentity $ Class.switchReal (Identity 1) (Identity 1)
 
 
-
 fill :: (Class.Floating a) => a -> Int -> Ptr a -> IO ()
 fill a n dstPtr = evalContT $ do
    nPtr <- Call.cint n
@@ -62,12 +76,20 @@
    incyPtr <- Call.cint 1
    liftIO $ BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr
 
+copyToTemp :: (Class.Floating a) => Int -> ForeignPtr a -> ContT r IO (Ptr a)
+copyToTemp n fptr = do
+   ptr <- ContT $ withForeignPtr fptr
+   tmpPtr <- Call.allocaArray n
+   liftIO $ copyBlock n ptr tmpPtr
+   return tmpPtr
+
+
 {- |
 In ColumnMajor:
 Copy a m-by-n-matrix with lda>=m and ldb>=m.
 -}
 copySubMatrix ::
-   (Storable a, Class.Floating a) =>
+   (Class.Floating a) =>
    Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
 copySubMatrix m n lda aPtr ldb bPtr = evalContT $ do
    uploPtr <- Call.char 'A'
@@ -78,7 +100,7 @@
    liftIO $ LapackGen.lacpy uploPtr mPtr nPtr aPtr ldaPtr bPtr ldbPtr
 
 copyTransposed ::
-   (Storable a, Class.Floating a) =>
+   (Class.Floating a) =>
    Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
 copyTransposed n m aPtr ldb bPtr = evalContT $ do
    nPtr <- Call.cint n
@@ -90,10 +112,57 @@
          (pointerSeq 1 aPtr)
          (pointerSeq ldb bPtr)
 
+
+{- |
+Copy a m-by-n-matrix to ColumnMajor order.
+-}
+copyToColumnMajor ::
+   (Class.Floating a) =>
+   Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
+copyToColumnMajor order m n aPtr bPtr =
+   case order of
+      RowMajor -> copyTransposed m n aPtr m bPtr
+      ColumnMajor -> copyBlock (m*n) aPtr bPtr
+
+copyToSubColumnMajor ::
+   (Class.Floating a) =>
+   Order -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
+copyToSubColumnMajor order m n aPtr ldb bPtr =
+   case order of
+      RowMajor -> copyTransposed m n aPtr ldb bPtr
+      ColumnMajor ->
+         if m==ldb
+           then copyBlock (m*n) aPtr bPtr
+           else copySubMatrix m n m aPtr ldb bPtr
+
+
 pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a]
 pointerSeq k ptr = iterate (flip advancePtr k) ptr
 
 
+allocArray :: (Shape.C sh, Storable a) => sh -> ContT r IO (Array sh a, Ptr a)
+allocArray sh = do
+   fptr <- liftIO $ mallocForeignPtrArray $ Shape.size sh
+   ptr <- ContT $ withForeignPtr fptr
+   return (Array sh fptr, ptr)
+
+
+allocHigherArray ::
+   (Shape.C sh, Class.Floating a) =>
+   sh -> Int -> Int -> Int -> ContT r IO (Array sh a, (Ptr a, Int))
+allocHigherArray shapeX m n nrhs = do
+   (x,xPtr) <- allocArray shapeX
+   if m>n
+      then do
+         tmpPtr <- Call.allocaArray (m*nrhs)
+         ContT $ \act -> do
+            r <- act (x,(tmpPtr,m))
+            copySubMatrix n nrhs m tmpPtr n xPtr
+            return r
+      else return (x,(xPtr,n))
+
+
+
 newtype Sum a = Sum {runSum :: Int -> Ptr a -> Int -> IO a}
 
 sum :: Class.Floating a => Int -> Ptr a -> Int -> IO a
@@ -110,7 +179,7 @@
    evalContT $ do
       nPtr <- Call.cint n
       incxPtr <- Call.cint incx
-      yPtr <- Call.real oneReal
+      yPtr <- Call.real one
       incyPtr <- Call.cint 0
       liftIO $ BlasReal.dot nPtr xPtr incxPtr yPtr incyPtr
 
@@ -135,3 +204,110 @@
             transPtr mPtr nPtr alphaPtr aPtr ldaPtr
             xPtr incxPtr betaPtr yPtr incyPtr
       liftIO $ peek yPtr
+
+
+product :: Class.Floating a => Int -> Ptr a -> Int -> IO a
+product n xPtr incx =
+   foldM (\x ptr -> do y <- peek ptr; return $! x*y) one $
+   take n $ pointerSeq incx xPtr
+
+
+newtype LACGV a = LACGV {getLACGV :: Ptr CInt -> Ptr a -> Ptr CInt -> IO ()}
+
+lacgv :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
+lacgv =
+   getLACGV $
+   Class.switchFloating
+      (LACGV $ const $ const $ const $ return ())
+      (LACGV $ const $ const $ const $ return ())
+      (LACGV LapackComplex.lacgv)
+      (LACGV LapackComplex.lacgv)
+
+
+multiplyMatrix ::
+   (Class.Floating a) =>
+   Order -> Order -> Int -> Int -> Int ->
+   ForeignPtr a -> ForeignPtr a -> Ptr a -> IO ()
+multiplyMatrix orderA orderB m k n a b cPtr = do
+   let lda = case orderA of RowMajor -> k; ColumnMajor -> m
+   let ldb = case orderB of RowMajor -> n; ColumnMajor -> k
+   let ldc = m
+   evalContT $ do
+      transaPtr <- Call.char $ transposeFromOrder orderA
+      transbPtr <- Call.char $ transposeFromOrder orderB
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      kPtr <- Call.cint k
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr a
+      ldaPtr <- Call.cint lda
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint ldb
+      betaPtr <- Call.number zero
+      ldcPtr <- Call.cint ldc
+      liftIO $
+         BlasGen.gemm
+            transaPtr transbPtr mPtr nPtr kPtr alphaPtr aPtr ldaPtr
+            bPtr ldbPtr betaPtr cPtr ldcPtr
+
+
+
+withAutoWorkspaceInfo ::
+   (Class.Floating a) =>
+   String -> (Ptr a -> Ptr CInt -> Ptr CInt -> IO ()) -> IO ()
+withAutoWorkspaceInfo name computation = evalContT $ do
+   infoPtr <- Call.alloca
+   liftIO $ withAutoWorkspace $ \workPtr lworkPtr ->
+      computation workPtr lworkPtr infoPtr
+   info <- liftIO $ fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $ printf "%s: deficient rank %d" name info
+
+withAutoWorkspace ::
+   (Class.Floating a) =>
+   (Ptr a -> Ptr CInt -> IO ()) -> IO ()
+withAutoWorkspace computation = evalContT $ do
+   lworkPtr <- Call.cint (-1)
+   lwork <- liftIO $ alloca $ \workPtr -> do
+      computation workPtr lworkPtr
+      ceilingSize <$> peek workPtr
+   workPtr <- Call.allocaArray lwork
+   liftIO $ poke lworkPtr $ fromIntegral lwork
+   liftIO $ computation workPtr lworkPtr
+
+
+newtype FromReal a = FromReal {getFromReal :: RealOf a -> a}
+
+fromReal :: (Class.Floating a) => RealOf a -> a
+fromReal =
+   getFromReal $
+   Class.switchFloating
+      (FromReal id)
+      (FromReal id)
+      (FromReal (:+0))
+      (FromReal (:+0))
+
+newtype RealPart a = RealPart {getRealPart :: a -> RealOf a}
+
+realPart :: (Class.Floating a) => a -> RealOf a
+realPart =
+   getRealPart $
+   Class.switchFloating
+      (RealPart id)
+      (RealPart id)
+      (RealPart Complex.realPart)
+      (RealPart Complex.realPart)
+
+
+newtype FuncArg b a = FuncArg {runFuncArg :: a -> b}
+
+ceilingSize :: (Class.Floating a) => a -> Int
+ceilingSize =
+   runFuncArg $
+   Class.switchFloating
+      (FuncArg ceiling)
+      (FuncArg ceiling)
+      (FuncArg $ ceiling . Complex.realPart)
+      (FuncArg $ ceiling . Complex.realPart)
diff --git a/src/Numeric/LAPACK/Singular.hs b/src/Numeric/LAPACK/Singular.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Singular.hs
@@ -0,0 +1,406 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Singular (
+   values,
+   decompose,
+   decomposeNarrow,
+   decomposeSquat,
+   determinantAbsolute,
+   leastSquaresMinimumNormRCond,
+   pseudoInverseRCond,
+   RealOf,
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix.Square (Square)
+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor,ColumnMajor))
+import Numeric.LAPACK.Matrix.Private (General, ZeroInt, zeroInt)
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Private
+         (RealOf, withAutoWorkspace, fromReal, allocArray, allocHigherArray,
+          copyToTemp, copyToColumnMajor, copyToSubColumnMajor)
+
+import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
+import qualified Numeric.LAPACK.FFI.Real as LapackReal
+import qualified Numeric.Netlib.Utility as Call
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable.Internal as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable.Internal (Array(Array))
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Foreign.Marshal.Array (allocaArray)
+import Foreign.Marshal.Alloc (alloca)
+import Foreign.C.Types (CInt, CChar)
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr, nullPtr)
+import Foreign.Storable (Storable, peek)
+
+import Control.Monad.Trans.Cont (evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative ((<$>))
+
+import Text.Printf (printf)
+
+import Data.Complex (Complex)
+
+
+values ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> Vector ZeroInt (RealOf a)
+values =
+   getValues $
+   Class.switchFloating
+      (Values valuesAux) (Values valuesAux)
+      (Values valuesAux) (Values valuesAux)
+
+type Values_ height width a =
+   General height width a -> Vector ZeroInt (RealOf a)
+
+newtype Values height width a = Values {getValues :: Values_ height width a}
+
+valuesAux ::
+   (Shape.C height, Shape.C width,
+    Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   Values_ height width a
+valuesAux (Array shape@(MatrixShape.General _order height width) a) =
+   Array.unsafeCreateWithSize
+      (zeroInt $ min (Shape.size height) (Shape.size width)) $ \mn sPtr -> do
+   let (m,n) = MatrixShape.dimensions shape
+   let lda = m
+   evalContT $ do
+      jobuPtr <- Call.char 'N'
+      jobvtPtr <- Call.char 'N'
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      aPtr <- copyToTemp (m*n) a
+      ldaPtr <- Call.cint lda
+      let uPtr = nullPtr
+      let vtPtr = nullPtr
+      lduPtr <- Call.cint m
+      ldvtPtr <- Call.cint n
+      liftIO $
+         withInfo "gesvd" $ \infoPtr ->
+         gesvd jobuPtr jobvtPtr mPtr nPtr
+            aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr
+
+
+determinantAbsolute ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a -> RealOf a
+determinantAbsolute =
+   getDeterminantAbsolute $
+   Class.switchFloating
+      (DeterminantAbsolute determinantAbsoluteAux)
+      (DeterminantAbsolute determinantAbsoluteAux)
+      (DeterminantAbsolute determinantAbsoluteAux)
+      (DeterminantAbsolute determinantAbsoluteAux)
+
+newtype DeterminantAbsolute f a =
+   DeterminantAbsolute {
+      getDeterminantAbsolute :: f a -> RealOf a
+   }
+
+determinantAbsoluteAux ::
+   (Shape.C height, Shape.C width,
+    Class.Floating a, RealOf a ~ ar, Class.Floating ar) =>
+   General height width a -> ar
+determinantAbsoluteAux = Vector.product . values
+
+
+decompose ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a ->
+   (Square height a, Vector ZeroInt (RealOf a), Square width a)
+decompose =
+   getDecompose $
+   Class.switchFloating
+      (Decompose decomposeAux)
+      (Decompose decomposeAux)
+      (Decompose decomposeAux)
+      (Decompose decomposeAux)
+
+newtype Decompose m f v g a =
+   Decompose {
+      getDecompose :: m a -> (f a, v (RealOf a), g a)
+   }
+
+decomposeAux ::
+   (Shape.C height, Shape.C width,
+    Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   General height width a ->
+   (Square height a, Vector ZeroInt ar, Square width a)
+decomposeAux (Array (MatrixShape.General order height width) a) =
+   unsafePerformIO $ evalContT $ do
+      (u,uPtr0) <- allocArray (MatrixShape.Square order height)
+      (vt,vtPtr0) <- allocArray (MatrixShape.Square order width)
+      let ((m,n),(uPtr,vtPtr)) =
+            case order of
+               RowMajor ->
+                  ((Shape.size width, Shape.size height), (vtPtr0,uPtr0))
+               ColumnMajor ->
+                  ((Shape.size height, Shape.size width), (uPtr0,vtPtr0))
+      let mn = min m n
+      let lda = m
+      jobuPtr <- Call.char 'A'
+      jobvtPtr <- Call.char 'A'
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      aPtr <- copyToTemp (m*n) a
+      ldaPtr <- Call.cint lda
+      (s,sPtr) <- allocArray (zeroInt mn)
+      lduPtr <- Call.cint m
+      ldvtPtr <- Call.cint n
+      liftIO $
+         withInfo "gesvd" $ \infoPtr ->
+         gesvd jobuPtr jobvtPtr mPtr nPtr
+            aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr
+      return (u, s, vt)
+
+
+decomposeSquat ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a ->
+   (Square height a, Vector height (RealOf a), General height width a)
+decomposeSquat a =
+   let (u,s,vt) = decomposeNarrow $ Matrix.transpose a
+   in  (Square.transpose vt, s, Matrix.transpose u)
+
+decomposeNarrow ::
+   (Shape.C height, Shape.C width, Class.Floating a) =>
+   General height width a ->
+   (General height width a, Vector width (RealOf a), Square width a)
+decomposeNarrow =
+   getDecompose $
+   Class.switchFloating
+      (Decompose decomposeThin)
+      (Decompose decomposeThin)
+      (Decompose decomposeThin)
+      (Decompose decomposeThin)
+
+decomposeThin ::
+   (Shape.C height, Shape.C width,
+    Class.Floating a, RealOf a ~ ar, Storable ar) =>
+   General height width a ->
+   (General height width a, Vector width ar, Square width a)
+decomposeThin (Array (MatrixShape.General order height width) a) =
+   unsafePerformIO $ do
+     Call.assert "Singular.decomposeThin: matrix is wider than high"
+         (Shape.size height >= Shape.size width)
+     evalContT $ do
+      (u,uPtr0) <- allocArray (MatrixShape.General order height width)
+      (vt,vtPtr0) <- allocArray (MatrixShape.Square order width)
+      let ((m,n),(uPtr,vtPtr)) =
+            case order of
+               RowMajor ->
+                  ((Shape.size width, Shape.size height), (vtPtr0,uPtr0))
+               ColumnMajor ->
+                  ((Shape.size height, Shape.size width), (uPtr0,vtPtr0))
+      let mn = min m n
+      let lda = m
+      jobuPtr <- Call.char 'S'
+      jobvtPtr <- Call.char 'S'
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      aPtr <- copyToTemp (m*n) a
+      ldaPtr <- Call.cint lda
+      (s,sPtr) <- allocArray width
+      lduPtr <- Call.cint m
+      ldvtPtr <- Call.cint mn
+      liftIO $
+         withInfo "gesvd" $ \infoPtr ->
+         gesvd jobuPtr jobvtPtr mPtr nPtr
+            aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr
+      return (u, s, vt)
+
+
+type GESVD_ ar a =
+   Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr CInt ->
+   Ptr a -> Ptr CInt -> Ptr ar ->
+   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> Int -> Ptr CInt -> IO ()
+
+newtype GESVD a = GESVD {getGESVD :: GESVD_ (RealOf a) a}
+
+gesvd :: Class.Floating a => GESVD_ (RealOf a) a
+gesvd =
+   getGESVD $
+   Class.switchFloating
+      (GESVD gesvdReal)
+      (GESVD gesvdReal)
+      (GESVD gesvdComplex)
+      (GESVD gesvdComplex)
+
+gesvdReal :: (Class.Real a) => GESVD_ a a
+gesvdReal jobuPtr jobvtPtr mPtr nPtr
+      aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr _mn infoPtr =
+   withAutoWorkspace $ \workPtr lworkPtr ->
+   LapackReal.gesvd jobuPtr jobvtPtr
+      mPtr nPtr aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr
+      workPtr lworkPtr infoPtr
+
+gesvdComplex :: (Class.Real a) => GESVD_ a (Complex a)
+gesvdComplex jobuPtr jobvtPtr
+      mPtr nPtr aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr mn infoPtr =
+   allocaArray (5*mn) $ \rworkPtr ->
+   withAutoWorkspace $ \workPtr lworkPtr ->
+   LapackComplex.gesvd jobuPtr jobvtPtr
+      mPtr nPtr aPtr ldaPtr sPtr uPtr lduPtr vtPtr ldvtPtr
+      workPtr lworkPtr rworkPtr infoPtr
+
+
+leastSquaresMinimumNormRCond ::
+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs, Class.Floating a) =>
+   RealOf a -> General height width a -> General height nrhs a ->
+   (Int, General width nrhs a)
+leastSquaresMinimumNormRCond =
+   getLeastSquaresMinimumNormRCond $
+   Class.switchFloating
+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)
+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)
+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)
+      (LeastSquaresMinimumNormRCond leastSquaresMinimumNormRCondAux)
+
+newtype LeastSquaresMinimumNormRCond f g h a =
+   LeastSquaresMinimumNormRCond {
+      getLeastSquaresMinimumNormRCond ::
+         RealOf a -> f a -> g a -> (Int, h a)
+   }
+
+leastSquaresMinimumNormRCondAux ::
+   (Shape.C height, Eq height, Shape.C width, Shape.C nrhs,
+    Class.Floating a, RealOf a ~ ar, Class.Floating ar) =>
+   ar -> General height width a -> General height nrhs a ->
+   (Int, General width nrhs a)
+leastSquaresMinimumNormRCondAux rcond
+   (Array (MatrixShape.General orderA heightA widthA) a)
+   (Array (MatrixShape.General orderB heightB widthB) b) =
+      unsafePerformIO $ do
+   Call.assert "leastSquaresMinimumNorm: height shapes mismatch"
+      (heightA == heightB)
+   let shapeX = MatrixShape.General ColumnMajor widthA widthB
+   let m = Shape.size heightA
+   let n = Shape.size widthA
+   let nrhs = Shape.size widthB
+   let mn = min m n
+   let aSize = m*n
+   let lda = m
+   evalContT $ do
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      nrhsPtr <- Call.cint nrhs
+      aPtr <- Call.allocaArray aSize
+      liftIO $ withForeignPtr a $ \asrcPtr ->
+         copyToColumnMajor orderA m n asrcPtr aPtr
+      ldaPtr <- Call.cint lda
+      (x,(tmpPtr,ldtmp)) <- allocHigherArray shapeX m n nrhs
+      ldtmpPtr <- Call.cint ldtmp
+      liftIO $ withForeignPtr b $ \bPtr ->
+         copyToSubColumnMajor orderB m nrhs bPtr ldtmp tmpPtr
+
+      sPtr <- Call.allocaArray mn
+      rcondPtr <- Call.number rcond
+      rankPtr <- Call.alloca
+      liftIO $
+         withInfo "gelss" $ \infoPtr ->
+         gelss mPtr nPtr nrhsPtr aPtr ldaPtr tmpPtr ldtmpPtr sPtr rcondPtr
+            rankPtr mn infoPtr
+
+      rank <- liftIO $ fromIntegral <$> peek rankPtr
+      return (rank, x)
+
+
+type GELSS_ ar a =
+   Ptr CInt -> Ptr CInt -> Ptr CInt ->
+   Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt ->
+   Ptr ar -> Ptr ar -> Ptr CInt -> Int -> Ptr CInt -> IO ()
+
+newtype GELSS a = GELSS {getGELSS :: GELSS_ (RealOf a) a}
+
+gelss :: Class.Floating a => GELSS_ (RealOf a) a
+gelss =
+   getGELSS $
+   Class.switchFloating
+      (GELSS gelssReal)
+      (GELSS gelssReal)
+      (GELSS gelssComplex)
+      (GELSS gelssComplex)
+
+gelssReal :: (Class.Real a) => GELSS_ a a
+gelssReal mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr
+      rankPtr _mn infoPtr =
+   withAutoWorkspace $ \workPtr lworkPtr ->
+   LapackReal.gelss
+      mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr
+      rankPtr workPtr lworkPtr infoPtr
+
+gelssComplex :: (Class.Real a) => GELSS_ a (Complex a)
+gelssComplex mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr
+      rankPtr mn infoPtr =
+   allocaArray (5*mn) $ \rworkPtr ->
+   withAutoWorkspace $ \workPtr lworkPtr ->
+   LapackComplex.gelss
+      mPtr nPtr nrhsPtr aPtr ldaPtr bPtr ldbPtr sPtr rcondPtr
+      rankPtr workPtr lworkPtr rworkPtr infoPtr
+
+
+pseudoInverseRCond ::
+   (Shape.C height, Eq height, Shape.C width, Eq width, Class.Floating a) =>
+   RealOf a -> General height width a -> (Int, General width height a)
+pseudoInverseRCond =
+   getPseudoInverseRCond $
+   Class.switchFloating
+      (PseudoInverseRCond pseudoInverseRCondAux)
+      (PseudoInverseRCond pseudoInverseRCondAux)
+      (PseudoInverseRCond pseudoInverseRCondAux)
+      (PseudoInverseRCond pseudoInverseRCondAux)
+
+newtype PseudoInverseRCond f g a =
+   PseudoInverseRCond {
+      getPseudoInverseRCond :: RealOf a -> f a -> (Int, g a)
+   }
+
+pseudoInverseRCondAux ::
+   (Shape.C height, Eq height, Shape.C width, Eq width,
+    Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   RealOf a -> General height width a -> (Int, General width height a)
+pseudoInverseRCondAux rcond a =
+   let (MatrixShape.General _ height width) = Array.shape a
+   in if Shape.size height < Shape.size width
+        then
+          let (u,s,vt) = decomposeSquat a
+              (rank,recipS) = recipSigma rcond s
+          in  (rank,
+               Matrix.multiply (Matrix.adjoint vt) $
+               Matrix.scaleRows recipS $
+               Square.toGeneral $ Square.adjoint u)
+        else
+          let (u,s,vt) = decomposeNarrow a
+              (rank,recipS) = recipSigma rcond s
+          in  (rank,
+               Matrix.multiply (Square.toGeneral $ Square.adjoint vt) $
+               Matrix.scaleRows recipS $ Matrix.adjoint u)
+
+recipSigma ::
+   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
+   ar -> Array sh ar -> (Int, Array sh a)
+recipSigma rcond sigmas =
+   case Array.toList sigmas of
+      [] -> (0, Array.map fromReal sigmas)
+      xs@(x:_) ->
+         let smin = x * rcond
+         in (length (takeWhile (>=smin) xs),
+             Array.map (\s -> if s>=smin then fromReal (recip s) else 0) sigmas)
+
+
+withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
+withInfo name computation = alloca $ \infoPtr -> do
+   computation infoPtr
+   info <- fromIntegral <$> peek infoPtr
+   case compare info (0::Int) of
+      EQ -> return ()
+      LT -> error $ printf "%s: illegal value in %d-th argument" name (-info)
+      GT -> error $ printf "%s: %d superdiagonals did not converge" name info
diff --git a/src/Numeric/LAPACK/Vector.hs b/src/Numeric/LAPACK/Vector.hs
--- a/src/Numeric/LAPACK/Vector.hs
+++ b/src/Numeric/LAPACK/Vector.hs
@@ -1,6 +1,7 @@
 module Numeric.LAPACK.Vector (
    Vector,
    fromList,
+   autoFromList,
    constant,
    dot,
    sum,
@@ -8,6 +9,8 @@
    norm1,
    norm2,
    argAbsMaximum,
+   argAbs1Maximum,
+   product,
    scale,
    add, sub ,
    mac,
@@ -18,6 +21,7 @@
    ) where
 
 import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Private as Private
 import Numeric.LAPACK.Private (RealOf, zero, one, minusOne, fill)
 
@@ -48,7 +52,7 @@
 import Data.Word (Word64)
 import Data.Bits (shiftR, (.&.))
 
-import Prelude hiding (sum)
+import Prelude hiding (sum, product)
 
 
 type Vector = Array
@@ -57,11 +61,14 @@
 fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a
 fromList = Array.fromList
 
+autoFromList :: (Storable a) => [a] -> Vector (Shape.ZeroBased Int) a
+autoFromList = Array.vectorFromList
 
-constant :: (Shape.C sh, Storable a, Class.Floating a) => sh -> a -> Vector sh a
-constant sh a = Array.unsafeCreate sh $ fill a (Shape.size sh)
 
+constant :: (Shape.C sh, Class.Floating a) => sh -> a -> Vector sh a
+constant sh a = Array.unsafeCreateWithSize sh $ fill a
 
+
 newtype Dot sh a = Dot {runDot :: Vector sh a -> Vector sh a -> a}
 
 dot ::
@@ -119,14 +126,24 @@
 sum (Array sh x) = unsafePerformIO $
    withForeignPtr x $ \xPtr -> Private.sum (Shape.size sh) xPtr 1
 
-norm1 :: (Shape.C sh, Class.Real a) => Vector sh a -> a
+norm1 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
 norm1 (Array sh x) = unsafePerformIO $
    evalContT $ do
       nPtr <- Call.cint $ Shape.size sh
       sxPtr <- ContT $ withForeignPtr x
       incxPtr <- Call.cint 1
-      liftIO $ BlasReal.asum nPtr sxPtr incxPtr
+      liftIO $ csum1 nPtr sxPtr incxPtr
 
+csum1 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
+csum1 =
+   getNorm $
+   Class.switchFloating
+      (Norm BlasReal.asum)
+      (Norm BlasReal.asum)
+      (Norm LapackComplex.sum1)
+      (Norm LapackComplex.sum1)
+
+
 {- |
 Sum of the absolute values of real numbers or components of complex numbers.
 For real numbers it is equivalent to 'norm1'.
@@ -170,12 +187,10 @@
 
 {- |
 Returns the index and value of the element with the maximal absolute value.
-The function does not strictly compare the absolute value of a complex number
-but the sum of the absolute complex components.
 Caution: It actually returns the value of the element, not its absolute value!
 -}
 argAbsMaximum ::
-   (Shape.C sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Class.Floating a) =>
    Vector sh a -> (Shape.Index sh, a)
 argAbsMaximum (Array sh x) = unsafePerformIO $
    evalContT $ do
@@ -183,28 +198,64 @@
       sxPtr <- ContT $ withForeignPtr x
       incxPtr <- Call.cint 1
       liftIO $ do
+         k <- fromIntegral . subtract 1 <$> absMax nPtr sxPtr incxPtr
+         xmax <- peekElemOff sxPtr k
+         return (Shape.indices sh !! k, xmax)
+
+newtype ArgMaximum a =
+   ArgMaximum {runArgMaximum :: Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt}
+
+absMax :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO CInt
+absMax =
+   runArgMaximum $
+   Class.switchFloating
+      (ArgMaximum BlasGen.iamax)
+      (ArgMaximum BlasGen.iamax)
+      (ArgMaximum LapackComplex.imax1)
+      (ArgMaximum LapackComplex.imax1)
+
+
+{- |
+Returns the index and value of the element with the maximal absolute value.
+The function does not strictly compare the absolute value of a complex number
+but the sum of the absolute complex components.
+Caution: It actually returns the value of the element, not its absolute value!
+-}
+argAbs1Maximum ::
+   (Shape.C sh, Class.Floating a) =>
+   Vector sh a -> (Shape.Index sh, a)
+argAbs1Maximum (Array sh x) = unsafePerformIO $
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size sh
+      sxPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      liftIO $ do
          k <- fromIntegral . subtract 1 <$> BlasGen.iamax nPtr sxPtr incxPtr
          xmax <- peekElemOff sxPtr k
          return (Shape.indices sh !! k, xmax)
 
 
+product :: (Shape.C sh, Class.Floating a) => Vector sh a -> a
+product (Array sh x) = unsafePerformIO $
+   withForeignPtr x $ \xPtr -> Private.product (Shape.size sh) xPtr 1
+
+
 scale, _scale ::
-   (Shape.C sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Class.Floating a) =>
    a -> Vector sh a -> Vector sh a
-scale alpha (Array sh x) = Array.unsafeCreate sh $ \syPtr -> do
+scale alpha (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr -> do
    evalContT $ do
       alphaPtr <- Call.number alpha
-      nPtr <- Call.cint $ Shape.size sh
+      nPtr <- Call.cint n
       sxPtr <- ContT $ withForeignPtr x
       incxPtr <- Call.cint 1
       incyPtr <- Call.cint 1
       liftIO $ BlasGen.copy nPtr sxPtr incxPtr syPtr incyPtr
       liftIO $ BlasGen.scal nPtr alphaPtr syPtr incyPtr
 
-_scale a (Array sh b) = Array.unsafeCreate sh $ \cPtr -> do
+_scale a (Array sh b) = Array.unsafeCreateWithSize sh $ \n cPtr -> do
    let m = 1
    let k = 1
-   let n = Shape.size sh
    evalContT $ do
       transaPtr <- Call.char 'N'
       transbPtr <- Call.char 'N'
@@ -224,18 +275,19 @@
             aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
 
 add, sub ::
-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Eq sh, Class.Floating a) =>
    Vector sh a -> Vector sh a -> Vector sh a
 add = mac one
 sub x y = mac minusOne y x
 
 mac ::
-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Eq sh, Class.Floating a) =>
    a -> Vector sh a -> Vector sh a -> Vector sh a
-mac alpha (Array shX x) (Array shY y) = Array.unsafeCreate shX $ \szPtr -> do
+mac alpha (Array shX x) (Array shY y) =
+      Array.unsafeCreateWithSize shX $ \n szPtr -> do
    Call.assert "mac: shapes mismatch" (shX == shY)
    evalContT $ do
-      nPtr <- Call.cint $ Shape.size shX
+      nPtr <- Call.cint n
       saPtr <- Call.number alpha
       sxPtr <- ContT $ withForeignPtr x
       incxPtr <- Call.cint 1
@@ -246,11 +298,11 @@
       liftIO $ BlasGen.axpy nPtr saPtr sxPtr incxPtr szPtr inczPtr
 
 mul ::
-   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Eq sh, Class.Floating a) =>
    Vector sh a -> Vector sh a -> Vector sh a
-mul (Array shA a) (Array shX x) = Array.unsafeCreate shX $ \yPtr -> do
+mul (Array shA a) (Array shX x) =
+      Array.unsafeCreateWithSize shX $ \n yPtr -> do
    Call.assert "mul: shapes mismatch" (shA == shX)
-   let n = Shape.size shX
    evalContT $ do
       transPtr <- Call.char 'N'
       mPtr <- Call.cint n
@@ -271,9 +323,8 @@
 
 
 outer ::
-   (Shape.C shx, Eq shx, Shape.C shy, Eq shy,
-    Storable a, Class.Floating a) =>
-   Vector shx a -> Vector shy a -> Array (MatrixShape.General shx shy) a
+   (Shape.C shx, Eq shx, Shape.C shy, Eq shy, Class.Floating a) =>
+   Vector shx a -> Vector shy a -> Matrix.General shx shy a
 outer (Array shX x) (Array shY y) =
    Array.unsafeCreate (MatrixShape.General MatrixShape.ColumnMajor shX shY) $
       \cPtr -> do
@@ -301,7 +352,7 @@
 newtype Conjugate sh a = Conjugate {getConjugate :: Vector sh a -> Vector sh a}
 
 conjugate ::
-   (Shape.C sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Class.Floating a) =>
    Vector sh a -> Vector sh a
 conjugate =
    getConjugate $
@@ -312,11 +363,11 @@
       (Conjugate complexConjugate)
 
 complexConjugate ::
-   (Shape.C sh, Storable a, Class.Real a) =>
+   (Shape.C sh, Class.Real a) =>
    Vector sh (Complex a) -> Vector sh (Complex a)
-complexConjugate (Array sh x) = Array.unsafeCreate sh $ \syPtr ->
+complexConjugate (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr ->
    evalContT $ do
-      nPtr <- Call.cint $ Shape.size sh
+      nPtr <- Call.cint n
       sxPtr <- ContT $ withForeignPtr x
       incxPtr <- Call.cint 1
       incyPtr <- Call.cint 1
@@ -339,11 +390,11 @@
 Only the least significant 47 bits of @seed@ are used.
 -}
 random ::
-   (Shape.C sh, Storable a, Class.Floating a) =>
+   (Shape.C sh, Class.Floating a) =>
    RandomDistribution -> sh -> Word64 -> Vector sh a
-random dist sh seed = Array.unsafeCreate sh $ \xPtr ->
+random dist sh seed = Array.unsafeCreateWithSize sh $ \n xPtr ->
    evalContT $ do
-      nPtr <- Call.cint $ Shape.size sh
+      nPtr <- Call.cint n
       distPtr <-
          Call.cint $
          case (getConst $ isComplexInFunctor xPtr, dist) of