diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) Henning Thielemann 2018
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/lapack.cabal b/lapack.cabal
new file mode 100644
--- /dev/null
+++ b/lapack.cabal
@@ -0,0 +1,66 @@
+Name:             lapack
+Version:          0.0
+License:          BSD3
+License-File:     LICENSE
+Author:           Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         http://hub.darcs.net/thielema/lapack/
+Category:         Data Structures
+Synopsis:         Numerical Linear Algebra using LAPACK
+Description:
+  This is a high-level interface to LAPACK.
+  .
+  Features:
+  .
+  * Based on @comfort-array@:
+    Allows to precisely express one-column or one-row matrices,
+    as well as dense, square, triangular, banded and symmetric matrices.
+  .
+  * Support all data types that are supported by LAPACK,
+    i.e. Float, Double, Complex Float, Complex Double
+  .
+  * No need for c2hs, hsc, Template Haskell or C helper functions
+  .
+  * Dependency only on BLAS and LAPACK, no GSL
+  .
+  * Separate formatting operator @(##)@:
+    Works better for tuples of matrices and vectors than 'show'.
+    'Show' is used for code one-liners
+    that can be copied back into Haskell modules.
+  .
+  See also: @hmatrix@.
+Tested-With:      GHC==7.4.2, GHC==7.8.4, GHC==8.2.2
+Cabal-Version:    >=1.6
+Build-Type:       Simple
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    http://hub.darcs.net/thielema/lapack/
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://hub.darcs.net/thielema/lapack/
+
+Library
+  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,
+    transformers >=0.3 && <0.6,
+    non-empty >=0.3 && <0.4,
+    utility-ht >=0.0.10 && <0.1,
+    base >=4.5 && <5
+
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   src
+  Exposed-Modules:
+    Numeric.LAPACK.Matrix
+    Numeric.LAPACK.Matrix.Shape
+    Numeric.LAPACK.Vector
+    Numeric.LAPACK.LinearSystem
+  Other-Modules:
+    Numeric.LAPACK.Matrix.Shape.Private
+    Numeric.LAPACK.Private
+    Numeric.LAPACK.Format
diff --git a/src/Numeric/LAPACK/Format.hs b/src/Numeric/LAPACK/Format.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Format.hs
@@ -0,0 +1,125 @@
+module Numeric.LAPACK.Format where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import Numeric.LAPACK.Matrix.Shape.Private (Order(RowMajor, ColumnMajor))
+
+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 (Array)
+
+import Foreign.Storable (Storable)
+
+import Text.Printf (PrintfArg, printf)
+
+import qualified Data.List.HT as ListHT
+import Data.Complex (Complex((:+)))
+
+
+infix 0 ##
+
+(##) :: (Format a) => a -> String -> IO ()
+a ## fmt = putStr $ unlines $ format fmt a
+
+
+class Format a where
+   format :: String -> a -> [String]
+
+instance Format Int where
+   format _fmt a = [show a]
+
+instance Format Float where
+   format fmt a = [printf fmt a]
+
+instance Format Double where
+   format fmt a = [printf fmt a]
+
+instance (PrintfArg a) => Format (Complex a) where
+   format fmt a = [printfComplex fmt a]
+
+instance (Format a, Format b) => Format (a,b) where
+   format fmt (a,b) = format fmt a ++ [""] ++ format fmt b
+
+instance (Format a, Format b, Format c) => Format (a,b,c) where
+   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
+   format = formatArray
+
+
+class (Shape.C sh) => FormatArray sh where
+   formatArray ::
+      (Storable a, 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]
+
+instance (Integral i) => FormatArray (Shape.OneBased i) where
+   formatArray fmt m = [unwords $ map (printfFloating fmt) $ Array.toList m]
+
+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]
+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
+
+instance
+   (Shape.C height, Shape.C width) =>
+      FormatArray (MatrixShape.Householder height width) where
+   formatArray = formatHouseholder
+
+formatHouseholder ::
+   (Shape.C height, Shape.C width, Storable a, 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  zipWith
+         (\row xs ->
+            concat $
+            zipWith (\col cell -> (if row==col then '|' else ' '):cell) [0..] $
+            zipWith (ListHT.padLeft ' ') strWidths xs)
+         [(0::Int)..] xss
+
+formatRows ::
+   (Class.Floating a, Shape.C height, Shape.C width) =>
+   String -> Order -> (height, width) -> [a] -> [[String]]
+formatRows fmt order (height,width) =
+   (case order of
+      RowMajor -> ListHT.sliceVertical (Shape.size width)
+      ColumnMajor -> ListHT.sliceHorizontal (Shape.size height)) .
+   map (printfFloating fmt)
+
+columnWidths :: [[[a]]] -> [Int]
+columnWidths xss =
+   case map (map length) xss of
+      [] -> []
+      w:ws -> foldl (zipWith max) w ws
+
+
+newtype Printf a = Printf {runPrintf :: String -> a -> String}
+
+printfFloating :: (Class.Floating a) => String -> a -> String
+printfFloating =
+   runPrintf $
+   Class.switchFloating
+      (Printf printf)
+      (Printf printf)
+      (Printf printfComplex)
+      (Printf printfComplex)
+
+printfComplex :: (PrintfArg a) => String -> Complex a -> String
+printfComplex fmt (r:+i) = printf (fmt ++ "+i" ++ fmt) r i
diff --git a/src/Numeric/LAPACK/LinearSystem.hs b/src/Numeric/LAPACK/LinearSystem.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/LinearSystem.hs
@@ -0,0 +1,478 @@
+{-# 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
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix.hs
@@ -0,0 +1,596 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+module Numeric.LAPACK.Matrix (
+   General,
+   (Format.##),
+   Format.Format,
+   Format.FormatArray,
+   ZeroInt, zeroInt,
+   transpose,
+   fromList,
+   identity,
+   diagonal, getDiagonal,
+   fromRows,    fromRowsWithSize,
+   fromColumns, fromColumnsWithSize,
+   singleRow,   singleColumn,
+   flattenRow,  flattenColumn,
+   pickRow, pickColumn,
+   takeRows, takeColumns,
+   dropRows, dropColumns,
+   reverseRows, reverseColumns,
+   fromRowMajor, toRowMajor, flatten,
+   (|||),
+   (===),
+
+   rowSums, columnSums,
+   scaleRows, scaleColumns,
+   multiply,
+   multiplyVector,
+
+   trace,
+   ) 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.Vector as Vector
+import Numeric.LAPACK.Matrix.Shape.Private
+         (Order(RowMajor, ColumnMajor), charFromOrder)
+import Numeric.LAPACK.Private
+         (zero, one, pointerSeq, copyTransposed, copySubMatrix, copyBlock)
+
+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 Data.Array.Comfort.Shape ((:+:)((:+:)))
+
+import Foreign.Marshal.Array (copyArray, advancePtr, pokeArray)
+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (Storable, poke, peek)
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+import qualified Data.NonEmpty as NonEmpty
+import Data.Foldable (forM_)
+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
+
+
+type ZeroInt = Shape.ZeroBased Int
+
+zeroInt :: Int -> ZeroInt
+zeroInt = Shape.ZeroBased
+
+
+fromList ::
+   (Shape.C height, Shape.C width, Storable a) =>
+   height -> width -> [a] -> General height width a
+fromList height width =
+   Array.fromList (MatrixShape.General RowMajor height width)
+
+
+type Vector = Array
+
+
+identity, _identity ::
+   (Shape.C sh, Storable a, 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
+
+diagonal ::
+   (Shape.C sh, Storable a, 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
+
+getDiagonal ::
+   (Shape.C sh, Eq sh, Storable a, 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
+
+
+singleRow :: Vector width a -> General () width a
+singleRow (Array sh fptr) =
+   Array (MatrixShape.General RowMajor () sh) fptr
+
+singleColumn :: Vector width a -> General width () a
+singleColumn (Array sh fptr) =
+   Array (MatrixShape.General ColumnMajor sh ()) fptr
+
+flattenRow :: General () width a -> Vector width a
+flattenRow (Array (MatrixShape.General _ () sh) fptr) = Array sh fptr
+
+flattenColumn :: General width () a -> Vector width a
+flattenColumn (Array (MatrixShape.General _ sh ()) fptr) = Array sh fptr
+
+
+fromRows ::
+   (Shape.C width, Eq width, Storable a) =>
+   NonEmpty.T [] (Vector width a) -> General ZeroInt width a
+fromRows (NonEmpty.Cons row rows) =
+   fromRowsWithSize (Array.shape row) (row:rows)
+
+fromRowsWithSize ::
+   (Shape.C width, Eq width, Storable a) =>
+   width -> [Vector width a] -> General ZeroInt width a
+fromRowsWithSize width rows =
+   Array.unsafeCreate
+      (MatrixShape.General RowMajor (zeroInt $ length rows) width)
+      (gather width rows)
+
+fromColumns ::
+   (Shape.C height, Eq height, Storable a) =>
+   NonEmpty.T [] (Vector height a) -> General height ZeroInt a
+fromColumns (NonEmpty.Cons column columns) =
+   fromColumnsWithSize (Array.shape column) (column:columns)
+
+fromColumnsWithSize ::
+   (Shape.C height, Eq height, Storable a) =>
+   height -> [Vector height a] -> General height ZeroInt a
+fromColumnsWithSize height columns =
+   Array.unsafeCreate
+      (MatrixShape.General ColumnMajor height (zeroInt $ length columns))
+      (gather height columns)
+
+gather ::
+   (Shape.C width, Eq width, Storable a) =>
+   width -> [Array width a] -> Ptr a -> IO ()
+gather width rows dstPtr =
+   let widthSize = Shape.size width
+   in forM_ (zip (pointerSeq widthSize dstPtr) rows) $
+         \(dstRowPtr, Array.Array rowWidth srcFPtr) ->
+         withForeignPtr srcFPtr $ \srcPtr -> do
+            Call.assert
+               "Matrix.fromRows/fromColumns: non-matching vector size"
+               (width == rowWidth)
+            copyArray dstRowPtr srcPtr widthSize
+
+
+pickRow ::
+   (Shape.C height, Shape.C width, Shape.Index height ~ ix,
+    Storable a, Class.Floating a) =>
+   General height width a -> ix -> Vector width a
+pickRow (Array (MatrixShape.General order height width) x) ix =
+   case order of
+      RowMajor -> pickConsecutive height width x ix
+      ColumnMajor -> pickScattered width height x ix
+
+pickColumn ::
+   (Shape.C height, Shape.C width, Shape.Index width ~ ix,
+    Storable a, Class.Floating a) =>
+   General height width a -> ix -> Vector height a
+pickColumn (Array (MatrixShape.General order height width) x) ix =
+   case order of
+      RowMajor -> pickScattered height width x ix
+      ColumnMajor -> pickConsecutive width height x ix
+
+pickConsecutive ::
+   (Shape.C height, Shape.C width, Shape.Index height ~ ix,
+    Storable 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
+      let offset = Shape.offset height ix
+      nPtr <- Call.cint n
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      incyPtr <- Call.cint 1
+      liftIO $
+         BlasGen.copy nPtr (advancePtr xPtr (n*offset)) incxPtr yPtr incyPtr
+
+pickScattered ::
+   (Shape.C height, Shape.C width, Shape.Index width ~ ix,
+    Storable 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
+      let offset = Shape.offset width ix
+      nPtr <- Call.cint n
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint $ Shape.size width
+      incyPtr <- Call.cint 1
+      liftIO $
+         BlasGen.copy nPtr (advancePtr xPtr offset) incxPtr yPtr incyPtr
+
+
+takeRows, dropRows ::
+   (Shape.C width, Storable a, 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 ->
+      withForeignPtr a $ \aPtr ->
+      case order of
+         RowMajor -> copyBlock (heightB*n) aPtr bPtr
+         ColumnMajor -> copySubMatrix heightB n heightA aPtr heightB bPtr
+
+dropRows k0
+      (Array (MatrixShape.General order (Shape.ZeroBased heightA) width) a) =
+   let k = min k0 heightA
+       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 ->
+      withForeignPtr a $ \aPtr ->
+      case order of
+         RowMajor -> copyBlock (heightB*n) (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) =>
+   Int -> General height ZeroInt a -> General height ZeroInt a
+takeColumns k = transpose . takeRows k . transpose
+dropColumns k = transpose . dropRows k . transpose
+
+
+-- alternative: laswp
+reverseRows ::
+   (Shape.C width, Storable a, 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
+      let n = Shape.size height
+      let m = Shape.size width
+      fwdPtr <- Call.bool True
+      nPtr <- Call.cint n
+      mPtr <- Call.cint m
+      kPtr <- Call.allocaArray n
+      aPtr <- ContT $ withForeignPtr a
+      liftIO $ do
+         copyBlock (n*m) 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) =>
+   General height ZeroInt a -> General height ZeroInt a
+reverseColumns = transpose . reverseRows . transpose
+
+
+fromRowMajor ::
+   (Shape.C height, Shape.C width, Storable a, 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) =>
+   General height width a -> Array (height,width) a
+toRowMajor (Array (MatrixShape.General order height width) x) =
+   let shape = (height, width)
+   in case order of
+         RowMajor -> Array shape x
+         ColumnMajor -> Array.unsafeCreate shape $ \yPtr -> evalContT $ do
+            let n = Shape.size width
+            let m = Shape.size height
+            nPtr <- Call.cint n
+            xPtr <- ContT $ withForeignPtr x
+            incxPtr <- Call.cint m
+            incyPtr <- Call.cint 1
+            liftIO $ sequence_ $ take m $
+               zipWith
+                  (\xkPtr ykPtr ->
+                     BlasGen.copy nPtr xkPtr incxPtr ykPtr incyPtr)
+                  (pointerSeq 1 xPtr)
+                  (pointerSeq n yPtr)
+
+flatten ::
+   (Shape.C height, Shape.C width, Storable a, 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
+
+
+infixl 3 |||
+infixl 2 ===
+
+(|||) ::
+   (Shape.C height, Eq height, Shape.C widtha, Shape.C widthb,
+    Storable a, Class.Floating a) =>
+   General height widtha a ->
+   General height widthb a ->
+   General height (widtha:+:widthb) a
+(|||)
+      (Array (MatrixShape.General orderA heightA widthA) a)
+      (Array (MatrixShape.General orderB heightB widthB) b) =
+   if heightA /= heightB
+      then error "(|||): mismatching heights"
+      else
+         case (orderA,orderB) of
+            (RowMajor,RowMajor) ->
+               Array.unsafeCreate
+                  (MatrixShape.General RowMajor heightA (widthA:+:widthB)) $
+               \cPtr -> evalContT $ do
+                  let n = Shape.size heightA
+                  let ma = Shape.size widthA
+                  let mb = Shape.size widthB
+                  let m = ma+mb
+                  maPtr <- Call.cint ma
+                  mbPtr <- Call.cint mb
+                  aPtr <- ContT $ withForeignPtr a
+                  bPtr <- ContT $ withForeignPtr b
+                  incxPtr <- Call.cint 1
+                  incyPtr <- Call.cint 1
+                  liftIO $
+                     sequence_ $ take n $
+                     zipWith3
+                        (\akPtr bkPtr ckPtr -> do
+                           BlasGen.copy maPtr akPtr incxPtr ckPtr incyPtr
+                           BlasGen.copy mbPtr bkPtr incxPtr
+                              (ckPtr `advancePtr` ma) incyPtr)
+                        (pointerSeq ma aPtr)
+                        (pointerSeq mb bPtr)
+                        (pointerSeq m cPtr)
+            (RowMajor,ColumnMajor) ->
+               Array.unsafeCreate
+                  (MatrixShape.General ColumnMajor heightA (widthA:+:widthB)) $
+               \cPtr -> evalContT $ do
+                  let n = Shape.size heightA
+                  let ma = Shape.size widthA
+                  let mb = Shape.size widthB
+                  aPtr <- ContT $ withForeignPtr a
+                  bPtr <- ContT $ withForeignPtr b
+                  liftIO $ do
+                     copyTransposed n ma aPtr n cPtr
+                     copyBlock (n*mb) bPtr (advancePtr cPtr (n*ma))
+            (ColumnMajor,RowMajor) ->
+               Array.unsafeCreate
+                  (MatrixShape.General ColumnMajor heightA (widthA:+:widthB)) $
+               \cPtr -> evalContT $ do
+                  let n = Shape.size heightA
+                  let ma = Shape.size widthA
+                  let mb = Shape.size widthB
+                  let volA = n*ma
+                  aPtr <- ContT $ withForeignPtr a
+                  bPtr <- ContT $ withForeignPtr b
+                  liftIO $ do
+                     copyBlock volA aPtr cPtr
+                     copyTransposed n mb bPtr n (advancePtr cPtr volA)
+            (ColumnMajor,ColumnMajor) ->
+               Array.unsafeCreate
+                  (MatrixShape.General ColumnMajor heightA (widthA:+:widthB)) $
+               \cPtr -> evalContT $ do
+                  let n = Shape.size heightA
+                  let na = n * Shape.size widthA
+                  let nb = n * Shape.size widthB
+                  naPtr <- Call.cint na
+                  nbPtr <- Call.cint nb
+                  aPtr <- ContT $ withForeignPtr a
+                  bPtr <- ContT $ withForeignPtr b
+                  incxPtr <- Call.cint 1
+                  incyPtr <- Call.cint 1
+                  liftIO $ do
+                     BlasGen.copy naPtr aPtr incxPtr cPtr incyPtr
+                     BlasGen.copy nbPtr bPtr incxPtr
+                        (cPtr `advancePtr` na) incyPtr
+
+(===) ::
+   (Shape.C width, Eq width, Shape.C heighta, Shape.C heightb,
+    Storable a, Class.Floating a) =>
+   General heighta width a ->
+   General heightb width a ->
+   General (heighta:+:heightb) width a
+(===) a b = transpose (transpose a ||| transpose b)
+
+
+rowSums ::
+   (Shape.C height, Shape.C width, Storable a, 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) =>
+   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) =>
+   Vector height a -> General height width a -> General height width a
+scaleRows
+   (Array heightX x) (Array shape@(MatrixShape.General order height width) a) =
+      Array.unsafeCreate shape $ \bPtr -> do
+   Call.assert "scaleRows: sizes mismatch" (heightX == height)
+   case order of
+      RowMajor -> evalContT $ do
+         let m = Shape.size height
+         let n = Shape.size width
+         alphaPtr <- Call.alloca
+         nPtr <- Call.cint n
+         xPtr <- ContT $ withForeignPtr x
+         aPtr <- ContT $ withForeignPtr a
+         incaPtr <- Call.cint 1
+         incbPtr <- Call.cint 1
+         liftIO $ sequence_ $ take m $
+            zipWith3
+               (\xkPtr akPtr bkPtr -> do
+                  poke alphaPtr =<< peek xkPtr
+                  BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr
+                  BlasGen.scal nPtr alphaPtr bkPtr incbPtr)
+               (pointerSeq 1 xPtr)
+               (pointerSeq n aPtr)
+               (pointerSeq n bPtr)
+      ColumnMajor -> evalContT $ do
+         let m = Shape.size width
+         let n = Shape.size height
+         transPtr <- Call.char 'N'
+         nPtr <- Call.cint n
+         klPtr <- Call.cint 0
+         kuPtr <- Call.cint 0
+         alphaPtr <- Call.number one
+         xPtr <- ContT $ withForeignPtr x
+         ldxPtr <- Call.cint 1
+         aPtr <- ContT $ withForeignPtr a
+         incaPtr <- Call.cint 1
+         betaPtr <- Call.number zero
+         incbPtr <- Call.cint 1
+         liftIO $ sequence_ $ take m $
+            zipWith
+               (\akPtr bkPtr ->
+                  BlasGen.gbmv transPtr
+                     nPtr nPtr klPtr kuPtr alphaPtr xPtr ldxPtr
+                     akPtr incaPtr betaPtr bkPtr incbPtr)
+               (pointerSeq n aPtr)
+               (pointerSeq n bPtr)
+
+scaleColumns ::
+   (Shape.C height, Shape.C width, Eq width, Storable a, 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/Shape.hs b/src/Numeric/LAPACK/Matrix/Shape.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Shape.hs
@@ -0,0 +1,5 @@
+module Numeric.LAPACK.Matrix.Shape (
+   General,
+   ) where
+
+import Numeric.LAPACK.Matrix.Shape.Private
diff --git a/src/Numeric/LAPACK/Matrix/Shape/Private.hs b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Shape/Private.hs
@@ -0,0 +1,127 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Matrix.Shape.Private where
+
+import qualified Data.Array.Comfort.Shape as Shape
+
+import Data.Tuple.HT (swap)
+
+
+data Order = RowMajor | ColumnMajor
+   deriving (Eq, Show)
+
+flipOrder :: Order -> Order
+flipOrder RowMajor = ColumnMajor
+flipOrder ColumnMajor = RowMajor
+
+charFromOrder :: Order -> Char
+charFromOrder RowMajor = 'T'
+charFromOrder ColumnMajor = 'N'
+
+
+data General height width =
+   General {
+      generalOrder :: Order,
+      generalHeight :: height,
+      generalWidth :: width
+   } deriving (Eq, Show)
+
+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)
+
+   offset (General RowMajor height width) =
+      Shape.offset (height,width)
+   offset (General ColumnMajor height width) =
+      Shape.offset (width,height) . swap
+   uncheckedOffset (General RowMajor height width) =
+      Shape.uncheckedOffset (height,width)
+   uncheckedOffset (General ColumnMajor height width) =
+      Shape.uncheckedOffset (width,height) . swap
+
+   sizeOffset (General RowMajor height width) =
+      Shape.sizeOffset (height,width)
+   sizeOffset (General ColumnMajor height width) =
+      Shape.sizeOffset (width,height) . swap
+   uncheckedSizeOffset (General RowMajor height width) =
+      Shape.uncheckedSizeOffset (height,width)
+   uncheckedSizeOffset (General ColumnMajor height width) =
+      Shape.uncheckedSizeOffset (width,height) . swap
+
+   inBounds (General _ height width) = Shape.inBounds (height,width)
+   size (General _ height width) = Shape.size (height,width)
+   uncheckedSize (General _ height width) = Shape.uncheckedSize (height,width)
+
+
+transpose :: General height width -> General width height
+transpose (General order height width) = General (flipOrder order) width height
+
+dimensions ::
+   (Shape.C height, Shape.C width) => General height width -> (Int, Int)
+dimensions (General order height width) =
+   case order of
+      RowMajor -> (Shape.size width, Shape.size height)
+      ColumnMajor -> (Shape.size height, Shape.size width)
+
+
+data Householder height width =
+   Householder {
+      householderOrder :: Order,
+      householderHeight :: height,
+      householderWidth :: width
+   } deriving (Eq, Show)
+
+data Reflector = Reflector deriving (Eq)
+data Triangular = Triangular deriving (Eq)
+
+householderPart ::
+   (Shape.C height, Shape.C width) =>
+   Householder height width ->
+   (Shape.Index height, Shape.Index width) -> Either Reflector Triangular
+householderPart (Householder _ height width) (r,c) =
+   if Shape.offset height r > Shape.offset width c
+     then Left Reflector
+     else Right Triangular
+
+instance
+   (Shape.C height, Shape.C width) =>
+      Shape.C (Householder height width) where
+
+   type Index (Householder height width) =
+            (Either Reflector Triangular,
+             (Shape.Index height, Shape.Index width))
+
+   indices sh@(Householder _ height width) =
+      map (\ix -> (householderPart sh ix, ix)) $
+      Shape.indices (height,width)
+
+   offset sh@(Householder order height width) (part,ix) =
+      if part == householderPart sh ix
+        then
+            case order of
+               RowMajor -> Shape.offset (height,width) ix
+               ColumnMajor -> Shape.offset (width,height) (swap ix)
+        else error "Shape.Householder.offset: wrong matrix part"
+   uncheckedOffset (Householder RowMajor height width) =
+      Shape.uncheckedOffset (height,width) . snd
+   uncheckedOffset (Householder ColumnMajor height width) =
+      Shape.uncheckedOffset (width,height) . swap . snd
+
+   sizeOffset sh@(Householder order height width) (part,ix) =
+      if part == householderPart sh ix
+        then
+            case order of
+               RowMajor -> Shape.sizeOffset (height,width) ix
+               ColumnMajor -> Shape.sizeOffset (width,height) (swap ix)
+        else error "Shape.Householder.sizeOffset: wrong matrix part"
+   uncheckedSizeOffset (Householder RowMajor height width) =
+      Shape.uncheckedSizeOffset (height,width) . snd
+   uncheckedSizeOffset (Householder ColumnMajor height width) =
+      Shape.uncheckedSizeOffset (width,height) . swap . snd
+
+   size (Householder _ height width) = Shape.size (height,width)
+   uncheckedSize (Householder _ height width) =
+      Shape.uncheckedSize (height,width)
+   inBounds sh@(Householder _ height width) (part,ix) =
+      Shape.inBounds (height,width) ix
+      &&
+      part == householderPart sh ix
diff --git a/src/Numeric/LAPACK/Private.hs b/src/Numeric/LAPACK/Private.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Private.hs
@@ -0,0 +1,137 @@
+{-# LANGUAGE TypeFamilies #-}
+module Numeric.LAPACK.Private where
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+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.Ptr (Ptr)
+import Foreign.Storable (Storable, peek)
+
+import Control.Monad.Trans.Cont (evalContT)
+import Control.Monad.IO.Class (liftIO)
+
+import Data.Functor.Identity (Identity(Identity, runIdentity))
+
+import Data.Complex (Complex)
+
+import Prelude hiding (sum)
+
+
+type family RealOf x
+
+type instance RealOf Float = Float
+type instance RealOf Double = Double
+type instance RealOf (Complex Float) = Float
+type instance RealOf (Complex Double) = Double
+
+
+zero, one, minusOne :: Class.Floating a => a
+zero =
+   runIdentity $
+   Class.switchFloating (Identity 0) (Identity 0) (Identity 0) (Identity 0)
+one =
+   runIdentity $
+   Class.switchFloating (Identity 1) (Identity 1) (Identity 1) (Identity 1)
+minusOne =
+   runIdentity $
+   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
+   srcPtr <- Call.number a
+   incxPtr <- Call.cint 0
+   incyPtr <- Call.cint 1
+   liftIO $ BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr
+
+
+copyBlock :: (Class.Floating a) => Int -> Ptr a -> Ptr a -> IO ()
+copyBlock n srcPtr dstPtr = evalContT $ do
+   nPtr <- Call.cint n
+   incxPtr <- Call.cint 1
+   incyPtr <- Call.cint 1
+   liftIO $ BlasGen.copy nPtr srcPtr incxPtr dstPtr incyPtr
+
+{- |
+In ColumnMajor:
+Copy a m-by-n-matrix with lda>=m and ldb>=m.
+-}
+copySubMatrix ::
+   (Storable 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'
+   mPtr <- Call.cint m
+   nPtr <- Call.cint n
+   ldaPtr <- Call.cint lda
+   ldbPtr <- Call.cint ldb
+   liftIO $ LapackGen.lacpy uploPtr mPtr nPtr aPtr ldaPtr bPtr ldbPtr
+
+copyTransposed ::
+   (Storable a, Class.Floating a) =>
+   Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
+copyTransposed n m aPtr ldb bPtr = evalContT $ do
+   nPtr <- Call.cint n
+   incaPtr <- Call.cint m
+   incbPtr <- Call.cint 1
+   liftIO $ sequence_ $ take m $
+      zipWith
+         (\akPtr bkPtr -> BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr)
+         (pointerSeq 1 aPtr)
+         (pointerSeq ldb bPtr)
+
+pointerSeq :: (Storable a) => Int -> Ptr a -> [Ptr a]
+pointerSeq k ptr = iterate (flip advancePtr k) ptr
+
+
+newtype Sum a = Sum {runSum :: Int -> Ptr a -> Int -> IO a}
+
+sum :: Class.Floating a => Int -> Ptr a -> Int -> IO a
+sum =
+   runSum $
+   Class.switchFloating
+      (Sum sumReal)
+      (Sum sumReal)
+      (Sum sumComplex)
+      (Sum sumComplex)
+
+sumReal :: Class.Real a => Int -> Ptr a -> Int -> IO a
+sumReal n xPtr incx =
+   evalContT $ do
+      nPtr <- Call.cint n
+      incxPtr <- Call.cint incx
+      yPtr <- Call.real oneReal
+      incyPtr <- Call.cint 0
+      liftIO $ BlasReal.dot nPtr xPtr incxPtr yPtr incyPtr
+
+sumComplex :: Class.Real a => Int -> Ptr (Complex a) -> Int -> IO (Complex a)
+sumComplex n xPtr incx =
+   evalContT $ do
+      transPtr <- Call.char 'N'
+      mPtr <- Call.cint 1
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number one
+      onePtr <- Call.number one
+      zeroincPtr <- Call.cint 0
+      aPtr <- Call.allocaArray n
+      ldaPtr <- Call.cint 1
+      incxPtr <- Call.cint incx
+      betaPtr <- Call.number zero
+      yPtr <- Call.alloca
+      incyPtr <- Call.cint 1
+      liftIO $ BlasGen.copy nPtr onePtr zeroincPtr aPtr incyPtr
+      liftIO $
+         BlasGen.gemv
+            transPtr mPtr nPtr alphaPtr aPtr ldaPtr
+            xPtr incxPtr betaPtr yPtr incyPtr
+      liftIO $ peek yPtr
diff --git a/src/Numeric/LAPACK/Vector.hs b/src/Numeric/LAPACK/Vector.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Vector.hs
@@ -0,0 +1,369 @@
+module Numeric.LAPACK.Vector (
+   Vector,
+   fromList,
+   constant,
+   dot,
+   sum,
+   absSum,
+   norm1,
+   norm2,
+   argAbsMaximum,
+   scale,
+   add, sub ,
+   mac,
+   mul,
+   outer,
+   conjugate,
+   random, RandomDistribution(..),
+   ) where
+
+import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
+import qualified Numeric.LAPACK.Private as Private
+import Numeric.LAPACK.Private (RealOf, zero, one, minusOne, fill)
+
+import qualified Numeric.LAPACK.FFI.Generic as LapackGen
+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 Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr)
+import Foreign.Storable (Storable, peek, peekElemOff, pokeElemOff)
+import Foreign.C.Types (CInt)
+
+import System.IO.Unsafe (unsafePerformIO)
+
+import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
+import Control.Monad.IO.Class (liftIO)
+import Control.Applicative (Const(Const,getConst), (<$>))
+
+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 Data.Complex (Complex)
+import Data.Word (Word64)
+import Data.Bits (shiftR, (.&.))
+
+import Prelude hiding (sum)
+
+
+type Vector = Array
+
+
+fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a
+fromList = Array.fromList
+
+
+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)
+
+
+newtype Dot sh a = Dot {runDot :: Vector sh a -> Vector sh a -> a}
+
+dot ::
+   (Shape.C sh, Eq sh, Class.Floating a) =>
+   Vector sh a -> Vector sh a -> a
+dot =
+   runDot $
+   Class.switchFloating
+      (Dot dotReal)
+      (Dot dotReal)
+      (Dot dotComplex)
+      (Dot dotComplex)
+
+dotReal ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   Vector sh a -> Vector sh a -> a
+dotReal (Array shX x) (Array shY y) = unsafePerformIO $ do
+   Call.assert "dot: shapes mismatch" (shX == shY)
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size shX
+      sxPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      syPtr <- ContT $ withForeignPtr y
+      incyPtr <- Call.cint 1
+      liftIO $ BlasReal.dot nPtr sxPtr incxPtr syPtr incyPtr
+
+{-
+We cannot use 'cdot' because Haskell's FFI
+does not support Complex numbers as return values.
+-}
+dotComplex ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
+dotComplex (Array shX x) (Array shY y) = unsafePerformIO $ do
+   Call.assert "dot: shapes mismatch" (shX == shY)
+   evalContT $ do
+      transPtr <- Call.char 'N'
+      mPtr <- Call.cint 1
+      nPtr <- Call.cint $ Shape.size shX
+      alphaPtr <- Call.number one
+      xPtr <- ContT $ withForeignPtr x
+      ldxPtr <- Call.cint 1
+      yPtr <- ContT $ withForeignPtr y
+      incyPtr <- Call.cint 1
+      betaPtr <- Call.number zero
+      zPtr <- Call.alloca
+      inczPtr <- Call.cint 1
+      liftIO $
+         BlasGen.gemv
+            transPtr mPtr nPtr alphaPtr xPtr ldxPtr
+            yPtr incyPtr betaPtr zPtr inczPtr
+      liftIO $ peek zPtr
+
+sum :: (Shape.C sh, Class.Floating a) => Vector sh a -> a
+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 (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
+
+{- |
+Sum of the absolute values of real numbers or components of complex numbers.
+For real numbers it is equivalent to 'norm1'.
+-}
+absSum :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
+absSum (Array sh x) = unsafePerformIO $
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size sh
+      sxPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      liftIO $ asum nPtr sxPtr incxPtr
+
+asum :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
+asum =
+   getNorm $
+   Class.switchFloating
+      (Norm BlasReal.asum) (Norm BlasReal.asum)
+      (Norm BlasComplex.casum) (Norm BlasComplex.casum)
+
+
+{- |
+Euclidean norm of a vector or Frobenius norm of a matrix.
+-}
+norm2 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
+norm2 (Array sh x) = unsafePerformIO $
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size sh
+      sxPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      liftIO $ nrm2 nPtr sxPtr incxPtr
+
+nrm2 :: Class.Floating a => Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)
+nrm2 =
+   getNorm $
+   Class.switchFloating
+      (Norm BlasReal.nrm2) (Norm BlasReal.nrm2)
+      (Norm BlasComplex.cnrm2) (Norm BlasComplex.cnrm2)
+
+newtype Norm a =
+   Norm {getNorm :: Ptr CInt -> Ptr a -> Ptr CInt -> IO (RealOf a)}
+
+{- |
+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) =>
+   Vector sh a -> (Shape.Index sh, a)
+argAbsMaximum (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)
+
+
+scale, _scale ::
+   (Shape.C sh, Storable a, Class.Floating a) =>
+   a -> Vector sh a -> Vector sh a
+scale alpha (Array sh x) = Array.unsafeCreate sh $ \syPtr -> do
+   evalContT $ do
+      alphaPtr <- Call.number alpha
+      nPtr <- Call.cint $ Shape.size sh
+      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
+   let m = 1
+   let k = 1
+   let n = Shape.size sh
+   evalContT $ do
+      transaPtr <- Call.char 'N'
+      transbPtr <- Call.char 'N'
+      mPtr <- Call.cint m
+      kPtr <- Call.cint k
+      nPtr <- Call.cint n
+      alphaPtr <- Call.number one
+      aPtr <- Call.number a
+      ldaPtr <- Call.cint m
+      bPtr <- ContT $ withForeignPtr b
+      ldbPtr <- Call.cint k
+      betaPtr <- Call.number zero
+      ldcPtr <- Call.cint m
+      liftIO $
+         BlasGen.gemm
+            transaPtr transbPtr mPtr nPtr kPtr alphaPtr
+            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
+
+add, sub ::
+   (Shape.C sh, Eq sh, Storable a, 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) =>
+   a -> Vector sh a -> Vector sh a -> Vector sh a
+mac alpha (Array shX x) (Array shY y) = Array.unsafeCreate shX $ \szPtr -> do
+   Call.assert "mac: shapes mismatch" (shX == shY)
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size shX
+      saPtr <- Call.number alpha
+      sxPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      syPtr <- ContT $ withForeignPtr y
+      incyPtr <- Call.cint 1
+      inczPtr <- Call.cint 1
+      liftIO $ BlasGen.copy nPtr syPtr incyPtr szPtr inczPtr
+      liftIO $ BlasGen.axpy nPtr saPtr sxPtr incxPtr szPtr inczPtr
+
+mul ::
+   (Shape.C sh, Eq sh, Storable a, Class.Floating a) =>
+   Vector sh a -> Vector sh a -> Vector sh a
+mul (Array shA a) (Array shX x) = Array.unsafeCreate shX $ \yPtr -> do
+   Call.assert "mul: shapes mismatch" (shA == shX)
+   let n = Shape.size shX
+   evalContT $ do
+      transPtr <- Call.char 'N'
+      mPtr <- Call.cint n
+      nPtr <- Call.cint n
+      klPtr <- Call.cint 0
+      kuPtr <- Call.cint 0
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr a
+      ldaPtr <- Call.cint 1
+      xPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      betaPtr <- Call.number zero
+      incyPtr <- Call.cint 1
+      liftIO $
+         BlasGen.gbmv transPtr
+            mPtr nPtr klPtr kuPtr alphaPtr aPtr ldaPtr
+            xPtr incxPtr betaPtr yPtr incyPtr
+
+
+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
+outer (Array shX x) (Array shY y) =
+   Array.unsafeCreate (MatrixShape.General MatrixShape.ColumnMajor shX shY) $
+      \cPtr -> do
+   let m = Shape.size shX
+   let n = Shape.size shY
+   evalContT $ do
+      transaPtr <- Call.char 'N'
+      transbPtr <- Call.char 'N'
+      mPtr <- Call.cint m
+      nPtr <- Call.cint n
+      kPtr <- Call.cint 1
+      alphaPtr <- Call.number one
+      aPtr <- ContT $ withForeignPtr x
+      ldaPtr <- Call.cint m
+      bPtr <- ContT $ withForeignPtr y
+      ldbPtr <- Call.cint 1
+      betaPtr <- Call.number zero
+      ldcPtr <- Call.cint m
+      liftIO $
+         BlasGen.gemm
+            transaPtr transbPtr mPtr nPtr kPtr alphaPtr
+            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
+
+
+newtype Conjugate sh a = Conjugate {getConjugate :: Vector sh a -> Vector sh a}
+
+conjugate ::
+   (Shape.C sh, Storable a, Class.Floating a) =>
+   Vector sh a -> Vector sh a
+conjugate =
+   getConjugate $
+   Class.switchFloating
+      (Conjugate id)
+      (Conjugate id)
+      (Conjugate complexConjugate)
+      (Conjugate complexConjugate)
+
+complexConjugate ::
+   (Shape.C sh, Storable a, Class.Real a) =>
+   Vector sh (Complex a) -> Vector sh (Complex a)
+complexConjugate (Array sh x) = Array.unsafeCreate sh $ \syPtr ->
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size sh
+      sxPtr <- ContT $ withForeignPtr x
+      incxPtr <- Call.cint 1
+      incyPtr <- Call.cint 1
+      liftIO $ do
+         BlasGen.copy nPtr sxPtr incxPtr syPtr incyPtr
+         LapackComplex.lacgv nPtr syPtr incyPtr
+
+
+data RandomDistribution =
+     UniformBox01
+   | UniformBoxPM1
+   | Normal
+   | UniformDisc
+   | UniformCircle
+   deriving (Eq, Ord, Show, Enum)
+
+{-
+@random distribution shape seed@
+
+Only the least significant 47 bits of @seed@ are used.
+-}
+random ::
+   (Shape.C sh, Storable a, Class.Floating a) =>
+   RandomDistribution -> sh -> Word64 -> Vector sh a
+random dist sh seed = Array.unsafeCreate sh $ \xPtr ->
+   evalContT $ do
+      nPtr <- Call.cint $ Shape.size sh
+      distPtr <-
+         Call.cint $
+         case (getConst $ isComplexInFunctor xPtr, dist) of
+            (_, UniformBox01) -> 1
+            (_, UniformBoxPM1) -> 2
+            (_, Normal) -> 3
+            (True, UniformDisc) -> 4
+            (True, UniformCircle) -> 5
+            (False, UniformDisc) -> 2
+            (False, UniformCircle) ->
+               error
+                  "Vector.random: UniformCircle not supported for real numbers"
+      iseedPtr <- Call.allocaArray 4
+      liftIO $ do
+         pokeElemOff iseedPtr 0 $ fromIntegral ((seed `shiftR` 35) .&. 0xFFF)
+         pokeElemOff iseedPtr 1 $ fromIntegral ((seed `shiftR` 23) .&. 0xFFF)
+         pokeElemOff iseedPtr 2 $ fromIntegral ((seed `shiftR` 11) .&. 0xFFF)
+         pokeElemOff iseedPtr 3 $ fromIntegral ((seed.&.0x7FF)*2+1)
+         LapackGen.larnv distPtr iseedPtr nPtr xPtr
+
+isComplexInFunctor :: (Class.Floating a) => f a -> Const Bool a
+isComplexInFunctor _ =
+   Class.switchFloating (Const False) (Const False) (Const True) (Const True)