diff --git a/lapack.cabal b/lapack.cabal
--- a/lapack.cabal
+++ b/lapack.cabal
@@ -1,6 +1,6 @@
 Cabal-Version:    2.2
 Name:             lapack
-Version:          0.5.0.3
+Version:          0.5.1
 License:          BSD-3-Clause
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -53,7 +53,7 @@
   test-module.list
 
 Source-Repository this
-  Tag:         0.5.0.3
+  Tag:         0.5.1
   Type:        darcs
   Location:    https://hub.darcs.net/thielema/lapack/
 
@@ -67,11 +67,12 @@
 
 Library
   Build-Depends:
+    comfort-blas >=0.0 && <0.1,
     lapack-ffi >=0.0.3 && <0.1,
     blas-ffi >=0.0 && <0.2,
     netlib-ffi >=0.1.1 && <0.2,
     comfort-array-shape >=0.0 && <0.1,
-    comfort-array >=0.5 && <0.6,
+    comfort-array >=0.5.2 && <0.6,
     guarded-allocation >=0.0 && <0.1,
     hyper >=0.1 && <0.3,
     blaze-html >=0.7 && <0.10,
@@ -186,7 +187,6 @@
     Numeric.LAPACK.Matrix.Array.Private
     Numeric.LAPACK.Matrix.Lazy.UpperTriangular
     Numeric.LAPACK.Matrix.ModifierTyped
-    Numeric.LAPACK.Matrix.Modifier
     Numeric.LAPACK.Matrix.RowMajor
     Numeric.LAPACK.Matrix.Private
     Numeric.LAPACK.Vector.Private
diff --git a/src/Numeric/LAPACK/Linear/LowerUpper.hs b/src/Numeric/LAPACK/Linear/LowerUpper.hs
--- a/src/Numeric/LAPACK/Linear/LowerUpper.hs
+++ b/src/Numeric/LAPACK/Linear/LowerUpper.hs
@@ -43,7 +43,7 @@
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Shape as ExtShape
 import Numeric.LAPACK.Matrix.Array.Private (Full)
-import Numeric.LAPACK.Matrix.Modifier (Transposition, Conjugation, Inversion)
+import Numeric.BLAS.Matrix.Modifier (Transposition, Conjugation, Inversion)
 
 import qualified Numeric.Netlib.Class as Class
 
diff --git a/src/Numeric/LAPACK/Linear/Plain.hs b/src/Numeric/LAPACK/Linear/Plain.hs
--- a/src/Numeric/LAPACK/Linear/Plain.hs
+++ b/src/Numeric/LAPACK/Linear/Plain.hs
@@ -57,7 +57,7 @@
 import Numeric.LAPACK.Matrix.Triangular.Basic (Lower, Upper)
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), Triangle(Triangle))
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition(NonTransposed, Transposed),
           Conjugation(NonConjugated, Conjugated),
           Inversion(NonInverted, Inverted))
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
@@ -136,7 +136,7 @@
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Plain as Plain
-import qualified Numeric.LAPACK.Matrix.Modifier as Mod
+import qualified Numeric.BLAS.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Matrix.Divide as Divide
 import qualified Numeric.LAPACK.Matrix.Multiply as Multiply
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
diff --git a/src/Numeric/LAPACK/Matrix/Array/Multiply.hs b/src/Numeric/LAPACK/Matrix/Array/Multiply.hs
--- a/src/Numeric/LAPACK/Matrix/Array/Multiply.hs
+++ b/src/Numeric/LAPACK/Matrix/Array/Multiply.hs
@@ -32,7 +32,7 @@
 import Numeric.LAPACK.Matrix.Array.Banded (Banded)
 import Numeric.LAPACK.Matrix.Array.Private (ArrayMatrix, diagTag)
 import Numeric.LAPACK.Matrix.Type.Private (Matrix)
-import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed,Transposed))
+import Numeric.BLAS.Matrix.Modifier (Transposition(NonTransposed,Transposed))
 import Numeric.LAPACK.Matrix.Shape.Omni (StripSingleton(StripBands, StripFilled))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked, deconsUnchecked))
diff --git a/src/Numeric/LAPACK/Matrix/Banded/Basic.hs b/src/Numeric/LAPACK/Matrix/Banded/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Banded/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Banded/Basic.hs
@@ -46,7 +46,7 @@
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor,ColumnMajor), transposeFromOrder, swapOnRowMajor,
           UnaryProxy, addOffDiagonals)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (zero, one)
 import Numeric.LAPACK.Matrix.Extent.Private (Extent)
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitian.hs b/src/Numeric/LAPACK/Matrix/BandedHermitian.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitian.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitian.hs
@@ -47,7 +47,7 @@
 import Numeric.LAPACK.Matrix.Array.Mosaic (FlexHermitian)
 import Numeric.LAPACK.Matrix.Array.Private (Full)
 import Numeric.LAPACK.Matrix.Layout.Private (Order, UnaryProxy, natFromProxy)
-import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed, Transposed))
+import Numeric.BLAS.Matrix.Modifier (Transposition(NonTransposed, Transposed))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf)
 
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs b/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitian/Basic.hs
@@ -37,7 +37,7 @@
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor,ColumnMajor), uploFromOrder,
           UnaryProxy, natFromProxy)
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition(NonTransposed, Transposed), transposeOrder,
           Conjugation(NonConjugated), conjugatedOnRowMajor)
 import Numeric.LAPACK.Vector (Vector)
@@ -434,10 +434,13 @@
          case transposed of
             Transposed -> (Banded.transpose upper, Banded.transpose lower)
             NonTransposed -> (lower,upper)
-   in VectorPriv.mac one
-         (Banded.multiplyFull
+   in VectorPriv.recheck $
+      Vector.mac one
+         (VectorPriv.uncheck $
+          Banded.multiplyFull
             (Banded.mapExtent ExtentPriv.fromSquare lowerT) b)
-         (Banded.multiplyFull
+         (VectorPriv.uncheck $
+          Banded.multiplyFull
             (Banded.mapExtent ExtentPriv.fromSquare upperT) b)
 
 takeUpper ::
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs b/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitian/Eigen.hs
@@ -12,7 +12,7 @@
 import qualified Numeric.LAPACK.Shape as ExtShape
 import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
 import Numeric.LAPACK.Matrix.Layout.Private (uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier (conjugatedOnRowMajor)
+import Numeric.BLAS.Matrix.Modifier (conjugatedOnRowMajor)
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf)
 import Numeric.LAPACK.Private
diff --git a/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs b/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/BandedHermitianPositiveDefinite/Linear.hs
@@ -16,7 +16,7 @@
 import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))
 import Numeric.LAPACK.Matrix.Mosaic.Private (copyTriangleToTemp)
 import Numeric.LAPACK.Matrix.Layout.Private (uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(Conjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(Conjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
 import Numeric.LAPACK.Scalar (RealOf, realPart)
 import Numeric.LAPACK.Private (copyBlock, withInfo, rankMsg, definiteMsg)
diff --git a/src/Numeric/LAPACK/Matrix/Basic.hs b/src/Numeric/LAPACK/Matrix/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Basic.hs
@@ -9,7 +9,7 @@
 import qualified Numeric.LAPACK.Private as Private
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), transposeFromOrder, flipOrder)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Matrix.Private
          (Full, Tall, Wide, Square, General, fromFull, ShapeInt, revealOrder)
 import Numeric.LAPACK.Vector (Vector)
@@ -417,7 +417,7 @@
       (RowMajor.recomplex .
        RowMajor.scaleColumns
          (RowMajor.tensorProduct (Left NonConjugated) x
-            (Vector.one Shape.Enumeration)) .
+            (Vector.one Shape.static)) .
        RowMajor.decomplex)
 
 scaleColumnsComplex ::
diff --git a/src/Numeric/LAPACK/Matrix/Divide.hs b/src/Numeric/LAPACK/Matrix/Divide.hs
--- a/src/Numeric/LAPACK/Matrix/Divide.hs
+++ b/src/Numeric/LAPACK/Matrix/Divide.hs
@@ -20,7 +20,7 @@
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Array.Private (Full)
 import Numeric.LAPACK.Matrix.Type.Private (scaleWithCheck)
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition(NonTransposed,Transposed),
           Inversion(Inverted))
 import Numeric.LAPACK.Vector (Vector)
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian.hs b/src/Numeric/LAPACK/Matrix/Hermitian.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian.hs
@@ -87,7 +87,7 @@
           SymmetricP)
 import Numeric.LAPACK.Matrix.Array.Private (Full, General, Square, packTag)
 import Numeric.LAPACK.Matrix.Layout.Private (Order)
-import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed, Transposed))
+import Numeric.BLAS.Matrix.Modifier (Transposition(NonTransposed, Transposed))
 import Numeric.LAPACK.Matrix.Private (ShapeInt)
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf, one)
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Basic.hs
@@ -22,7 +22,7 @@
           withPacking, noLabel, applyFuncPair, triArg)
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor,ColumnMajor), uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition(NonTransposed, Transposed),
           Conjugation(Conjugated), conjugatedOnRowMajor)
 import Numeric.LAPACK.Vector (Vector)
diff --git a/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs b/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs
--- a/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs
+++ b/src/Numeric/LAPACK/Matrix/Hermitian/Eigen.hs
@@ -12,7 +12,7 @@
 import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian, HermitianP)
 import Numeric.LAPACK.Matrix.Square.Basic (Square)
 import Numeric.LAPACK.Matrix.Layout.Private (Order(ColumnMajor), uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier (conjugatedOnRowMajor)
+import Numeric.BLAS.Matrix.Modifier (conjugatedOnRowMajor)
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf)
 import Numeric.LAPACK.Private
diff --git a/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs b/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
--- a/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
+++ b/src/Numeric/LAPACK/Matrix/HermitianPositiveDefinite/Linear.hs
@@ -18,7 +18,7 @@
          (withPackingLinear, label, applyFuncPair, triArg, copyTriangleToTemp)
 import Numeric.LAPACK.Matrix.Mosaic.Basic (takeDiagonal)
 import Numeric.LAPACK.Matrix.Layout.Private (uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(Conjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(Conjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
 import Numeric.LAPACK.Linear.Private (solver)
 import Numeric.LAPACK.Scalar (RealOf, realPart, zero)
diff --git a/src/Numeric/LAPACK/Matrix/Layout/Private.hs b/src/Numeric/LAPACK/Matrix/Layout/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Layout/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Layout/Private.hs
@@ -3,10 +3,14 @@
 {-# LANGUAGE EmptyDataDecls #-}
 {-# LANGUAGE StandaloneDeriving #-}
 {-# LANGUAGE GADTs #-}
-module Numeric.LAPACK.Matrix.Layout.Private where
+module Numeric.LAPACK.Matrix.Layout.Private (
+   module Numeric.LAPACK.Matrix.Layout.Private,
+   module Numeric.BLAS.Matrix.Layout,
+   ) where
 
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import Numeric.LAPACK.Matrix.Extent.Private (Extent)
+import Numeric.BLAS.Matrix.Layout (Order(..), flipOrder, transposeFromOrder)
 
 import qualified Type.Data.Num.Unary.Literal as TypeNum
 import qualified Type.Data.Num.Unary.Proof as Proof
@@ -25,21 +29,6 @@
 import Data.Tuple.HT (mapSnd, swap)
 import Data.Bool.HT (if')
 
-
-data Order = RowMajor | ColumnMajor
-   deriving (Eq, Show)
-
-instance NFData Order where
-   rnf RowMajor = ()
-   rnf ColumnMajor = ()
-
-flipOrder :: Order -> Order
-flipOrder RowMajor = ColumnMajor
-flipOrder ColumnMajor = RowMajor
-
-transposeFromOrder :: Order -> Char
-transposeFromOrder RowMajor = 'T'
-transposeFromOrder ColumnMajor = 'N'
 
 swapOnRowMajor :: Order -> (a,a) -> (a,a)
 swapOnRowMajor order =
diff --git a/src/Numeric/LAPACK/Matrix/Modifier.hs b/src/Numeric/LAPACK/Matrix/Modifier.hs
deleted file mode 100644
--- a/src/Numeric/LAPACK/Matrix/Modifier.hs
+++ /dev/null
@@ -1,50 +0,0 @@
-module Numeric.LAPACK.Matrix.Modifier where
-
-import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor,ColumnMajor), flipOrder)
-
-import Data.Monoid (Monoid, mempty, mappend)
-import Data.Semigroup (Semigroup, (<>))
-
-
-
-data Transposition = NonTransposed | Transposed
-   deriving (Eq, Show, Enum, Bounded)
-
-instance Semigroup Transposition where
-   x<>y = if x==y then NonTransposed else Transposed
-
-instance Monoid Transposition where
-   mempty = NonTransposed
-   mappend = (<>)
-
-transposeOrder :: Transposition -> Order -> Order
-transposeOrder NonTransposed = id
-transposeOrder Transposed = flipOrder
-
-
-
-data Conjugation = NonConjugated | Conjugated
-   deriving (Eq, Show, Enum, Bounded)
-
-instance Semigroup Conjugation where
-   x<>y = if x==y then NonConjugated else Conjugated
-
-instance Monoid Conjugation where
-   mempty = NonConjugated
-   mappend = (<>)
-
-conjugatedOnRowMajor :: Order -> Conjugation
-conjugatedOnRowMajor RowMajor = Conjugated
-conjugatedOnRowMajor ColumnMajor = NonConjugated
-
-
-data Inversion = NonInverted | Inverted
-   deriving (Eq, Show, Enum, Bounded)
-
-instance Semigroup Inversion where
-   x<>y = if x==y then NonInverted else Inverted
-
-instance Monoid Inversion where
-   mempty = NonInverted
-   mappend = (<>)
diff --git a/src/Numeric/LAPACK/Matrix/ModifierTyped.hs b/src/Numeric/LAPACK/Matrix/ModifierTyped.hs
--- a/src/Numeric/LAPACK/Matrix/ModifierTyped.hs
+++ b/src/Numeric/LAPACK/Matrix/ModifierTyped.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Matrix.ModifierTyped where
 
-import qualified Numeric.LAPACK.Matrix.Modifier as Mod
+import qualified Numeric.BLAS.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import Numeric.LAPACK.Matrix.Extent.Private (Extent, Small)
 
diff --git a/src/Numeric/LAPACK/Matrix/Mosaic/Private.hs b/src/Numeric/LAPACK/Matrix/Mosaic/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Mosaic/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Mosaic/Private.hs
@@ -9,7 +9,7 @@
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor,ColumnMajor), flipOrder, uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
 import Numeric.LAPACK.Scalar (zero)
 import Numeric.LAPACK.Shape.Private (Unchecked(Unchecked))
diff --git a/src/Numeric/LAPACK/Matrix/Multiply.hs b/src/Numeric/LAPACK/Matrix/Multiply.hs
--- a/src/Numeric/LAPACK/Matrix/Multiply.hs
+++ b/src/Numeric/LAPACK/Matrix/Multiply.hs
@@ -14,7 +14,7 @@
 import qualified Numeric.LAPACK.Matrix.Class as MatrixClass
 import qualified Numeric.LAPACK.Matrix.Type.Private as Matrix
 import qualified Numeric.LAPACK.Matrix.Basic as FullBasic
-import qualified Numeric.LAPACK.Matrix.Modifier as Mod
+import qualified Numeric.BLAS.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
@@ -24,7 +24,7 @@
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Array.Private (Full)
 import Numeric.LAPACK.Matrix.Type.Private (Matrix, scaleWithCheck)
-import Numeric.LAPACK.Matrix.Modifier (Transposition(NonTransposed,Transposed))
+import Numeric.BLAS.Matrix.Modifier (Transposition(NonTransposed,Transposed))
 import Numeric.LAPACK.Vector (Vector)
 
 import qualified Numeric.Netlib.Class as Class
diff --git a/src/Numeric/LAPACK/Matrix/Permutation.hs b/src/Numeric/LAPACK/Matrix/Permutation.hs
--- a/src/Numeric/LAPACK/Matrix/Permutation.hs
+++ b/src/Numeric/LAPACK/Matrix/Permutation.hs
@@ -20,7 +20,7 @@
 import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
-import qualified Numeric.LAPACK.Matrix.Modifier as Mod
+import qualified Numeric.BLAS.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Permutation.Private as Plain
 import qualified Numeric.LAPACK.Permutation as Perm
 import Numeric.LAPACK.Matrix.Type.Private (Matrix(Permutation))
diff --git a/src/Numeric/LAPACK/Matrix/Plain.hs b/src/Numeric/LAPACK/Matrix/Plain.hs
--- a/src/Numeric/LAPACK/Matrix/Plain.hs
+++ b/src/Numeric/LAPACK/Matrix/Plain.hs
@@ -66,7 +66,7 @@
 import Numeric.LAPACK.Matrix.Private
          (Full, Tall, Wide, General, ShapeInt, shapeInt,
           mapExtent, fromFull, generalizeTall, generalizeWide)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated,Conjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated,Conjugated))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (zero, one)
 import Numeric.LAPACK.Private (pointerSeq, fill, copySubMatrix, copyBlock)
diff --git a/src/Numeric/LAPACK/Matrix/RowMajor.hs b/src/Numeric/LAPACK/Matrix/RowMajor.hs
--- a/src/Numeric/LAPACK/Matrix/RowMajor.hs
+++ b/src/Numeric/LAPACK/Matrix/RowMajor.hs
@@ -1,23 +1,33 @@
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
-module Numeric.LAPACK.Matrix.RowMajor where
+module Numeric.LAPACK.Matrix.RowMajor (
+   Matrix,
+   Matrix.takeRow,
+   Matrix.takeColumn,
+   Matrix.fromRows,
+   Matrix.tensorProduct,
+   Matrix.decomplex,
+   Matrix.recomplex,
+   Matrix.scaleRows,
+   Matrix.scaleColumns,
+   kronecker,
+   ) where
 
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
-import qualified Numeric.LAPACK.Private as Private
-import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor, ColumnMajor))
-import Numeric.LAPACK.Matrix.Private (Full, ShapeInt, shapeInt)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated,Conjugated))
-import Numeric.LAPACK.Scalar (zero, one)
-import Numeric.LAPACK.Private (ComplexPart, pointerSeq)
+import Numeric.LAPACK.Matrix.Private (Full)
 
+import qualified Numeric.BLAS.Matrix.RowMajor as Matrix
+import Numeric.BLAS.Matrix.RowMajor (Matrix)
+import Numeric.BLAS.Matrix.Layout (Order(RowMajor, ColumnMajor))
+import Numeric.BLAS.Scalar (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 Foreign.Marshal.Array (copyArray, advancePtr)
-import Foreign.ForeignPtr (withForeignPtr, castForeignPtr)
-import Foreign.Storable (Storable)
+import Foreign.Marshal.Array (advancePtr)
+import Foreign.ForeignPtr (withForeignPtr)
 
 import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
 import Control.Monad.IO.Class (liftIO)
@@ -27,150 +37,7 @@
 import qualified Data.Array.Comfort.Shape as Shape
 import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
-import Data.Complex (Complex)
 import Data.Foldable (forM_)
-
-
-type Matrix height width = Array (height,width)
-type Vector = Array
-
-takeRow ::
-   (Shape.Indexed height, Shape.C width, Shape.Index height ~ ix,
-    Storable a) =>
-   ix -> Matrix height width a -> Vector width a
-takeRow ix (Array (height,width) x) =
-   Array.unsafeCreateWithSize width $ \n yPtr ->
-   withForeignPtr x $ \xPtr ->
-      copyArray yPtr (advancePtr xPtr (n * Shape.offset height ix)) n
-
-takeColumn ::
-   (Shape.C height, Shape.Indexed width, Shape.Index width ~ ix,
-    Class.Floating a) =>
-   ix -> Matrix height width a -> Vector height a
-takeColumn ix (Array (height,width) x) =
-   Array.unsafeCreateWithSize height $ \n yPtr -> evalContT $ do
-      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
-
-
-fromRows ::
-   (Shape.C width, Eq width, Storable a) =>
-   width -> [Vector width a] -> Matrix ShapeInt width a
-fromRows width rows =
-   Array.unsafeCreate (shapeInt $ length rows, width) $ \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: non-matching vector size"
-               (width == rowWidth)
-            copyArray dstRowPtr srcPtr widthSize
-
-
-tensorProduct ::
-   (Shape.C height, Shape.C width, Class.Floating a) =>
-   Either Conjugation Conjugation ->
-   Vector height a -> Vector width a -> Matrix height width a
-tensorProduct side (Array height x) (Array width y) =
-   Array.unsafeCreate (height,width) $ \cPtr -> do
-   let m = Shape.size width
-   let n = Shape.size height
-   let trans conjugated =
-         case conjugated of NonConjugated -> 'T'; Conjugated -> 'C'
-   let ((transa,transb),(lda,ldb)) =
-         case side of
-            Left c -> ((trans c, 'N'),(1,1))
-            Right c -> (('N', trans c),(m,n))
-   evalContT $ do
-      transaPtr <- Call.char transa
-      transbPtr <- Call.char transb
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      kPtr <- Call.cint 1
-      alphaPtr <- Call.number one
-      aPtr <- ContT $ withForeignPtr y
-      ldaPtr <- Call.leadingDim lda
-      bPtr <- ContT $ withForeignPtr x
-      ldbPtr <- Call.leadingDim ldb
-      betaPtr <- Call.number zero
-      ldcPtr <- Call.leadingDim m
-      liftIO $
-         BlasGen.gemm
-            transaPtr transbPtr mPtr nPtr kPtr alphaPtr
-            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
-
-
-decomplex ::
-   (Class.Real a) =>
-   Matrix height width (Complex a) ->
-   Matrix height (width, Shape.Enumeration ComplexPart) a
-decomplex (Array (height,width) a) =
-   Array (height, (width, Shape.Enumeration)) (castForeignPtr a)
-
-recomplex ::
-   (Class.Real a) =>
-   Matrix height (width, Shape.Enumeration ComplexPart) a ->
-   Matrix height width (Complex a)
-recomplex (Array (height, (width, Shape.Enumeration)) a) =
-   Array (height,width) (castForeignPtr a)
-
-
-scaleRows ::
-   (Shape.C height, Eq height, Shape.C width, Class.Floating a) =>
-   Vector height a -> Matrix height width a -> Matrix height width a
-scaleRows (Array heightX x) (Array shape@(height,width) a) =
-      Array.unsafeCreate shape $ \bPtr -> do
-   Call.assert "scaleRows: sizes mismatch" (heightX == height)
-   evalContT $ do
-      let m = Shape.size height
-      let n = Shape.size width
-      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
-               BlasGen.copy nPtr akPtr incaPtr bkPtr incbPtr
-               BlasGen.scal nPtr xkPtr bkPtr incbPtr)
-            (pointerSeq 1 xPtr)
-            (pointerSeq n aPtr)
-            (pointerSeq n bPtr)
-
-scaleColumns ::
-   (Shape.C height, Shape.C width, Eq width, Class.Floating a) =>
-   Vector width a -> Matrix height width a -> Matrix height width a
-scaleColumns (Array widthX x) (Array shape@(height,width) a) =
-      Array.unsafeCreate shape $ \bPtr -> do
-   Call.assert "scaleColumns: sizes mismatch" (widthX == width)
-   evalContT $ do
-      let m = Shape.size height
-      let n = Shape.size width
-      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.leadingDim 1
-      aPtr <- ContT $ withForeignPtr a
-      incaPtr <- Call.cint 1
-      betaPtr <- Call.number zero
-      incbPtr <- Call.cint 1
-      liftIO $ sequence_ $ take m $
-         zipWith
-            (\akPtr bkPtr ->
-               Private.gbmv transPtr
-                  nPtr nPtr klPtr kuPtr alphaPtr xPtr ldxPtr
-                  akPtr incaPtr betaPtr bkPtr incbPtr)
-            (pointerSeq n aPtr)
-            (pointerSeq n bPtr)
 
 
 kronecker ::
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs b/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Symmetric/Basic.hs
@@ -18,7 +18,7 @@
          (withPacking, noLabel, applyFuncPair, triArg)
 import Numeric.LAPACK.Matrix.Layout.Private
          (MirrorSingleton(SimpleMirror), Order, uploFromOrder)
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (zero)
 import Numeric.LAPACK.Private (fill)
diff --git a/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs b/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs
--- a/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs
+++ b/src/Numeric/LAPACK/Matrix/Symmetric/Unified.hs
@@ -21,7 +21,7 @@
           Labelled(Labelled))
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor,ColumnMajor), uploFromOrder, uploOrder)
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Conjugation(NonConjugated, Conjugated), conjugatedOnRowMajor,
           Transposition(NonTransposed, Transposed), transposeOrder)
 import Numeric.LAPACK.Matrix.Private (Full, General, Square)
diff --git a/src/Numeric/LAPACK/Orthogonal/Householder.hs b/src/Numeric/LAPACK/Orthogonal/Householder.hs
--- a/src/Numeric/LAPACK/Orthogonal/Householder.hs
+++ b/src/Numeric/LAPACK/Orthogonal/Householder.hs
@@ -31,12 +31,12 @@
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
 import qualified Numeric.LAPACK.Matrix.Extent.Strict as ExtentStrict
 import qualified Numeric.LAPACK.Matrix.Extent as Extent
-import qualified Numeric.LAPACK.Matrix.Modifier as Mod
+import qualified Numeric.BLAS.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Shape as ExtShape
 import Numeric.LAPACK.Orthogonal.Plain (Householder)
 import Numeric.LAPACK.Matrix.Array.Mosaic (Upper)
 import Numeric.LAPACK.Matrix.Array.Private (Full, Square)
-import Numeric.LAPACK.Matrix.Modifier (Transposition, Conjugation)
+import Numeric.BLAS.Matrix.Modifier (Transposition, Conjugation)
 
 import qualified Numeric.Netlib.Class as Class
 
diff --git a/src/Numeric/LAPACK/Orthogonal/Plain.hs b/src/Numeric/LAPACK/Orthogonal/Plain.hs
--- a/src/Numeric/LAPACK/Orthogonal/Plain.hs
+++ b/src/Numeric/LAPACK/Orthogonal/Plain.hs
@@ -29,7 +29,7 @@
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), sideSwapFromOrder)
 import Numeric.LAPACK.Matrix.Extent.Private (Extent)
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition(NonTransposed, Transposed),
           Conjugation(NonConjugated, Conjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
diff --git a/src/Numeric/LAPACK/Permutation.hs b/src/Numeric/LAPACK/Permutation.hs
--- a/src/Numeric/LAPACK/Permutation.hs
+++ b/src/Numeric/LAPACK/Permutation.hs
@@ -18,7 +18,7 @@
 
 import qualified Numeric.LAPACK.Permutation.Private as Plain
 import qualified Numeric.LAPACK.Matrix.Array.Private as ArrMatrix
-import qualified Numeric.LAPACK.Matrix.Modifier as Mod
+import qualified Numeric.BLAS.Matrix.Modifier as Mod
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import Numeric.LAPACK.Permutation.Private (Permutation)
 import Numeric.LAPACK.Matrix.Array.Private (Full, Square)
diff --git a/src/Numeric/LAPACK/Permutation/Private.hs b/src/Numeric/LAPACK/Permutation/Private.hs
--- a/src/Numeric/LAPACK/Permutation/Private.hs
+++ b/src/Numeric/LAPACK/Permutation/Private.hs
@@ -7,7 +7,7 @@
 import qualified Numeric.LAPACK.Output as Output
 import Numeric.LAPACK.Output (Output, formatAligned)
 import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor, ColumnMajor))
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition(NonTransposed,Transposed),
           Inversion(NonInverted,Inverted))
 import Numeric.LAPACK.Matrix.Private (Full, Square)
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
@@ -11,7 +11,7 @@
 import qualified Numeric.BLAS.FFI.Generic as BlasGen
 import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))
+import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))
 import Numeric.LAPACK.Scalar (RealOf, zero, one, isZero)
 
 import qualified Foreign.Marshal.Array.Guarded as ForeignArray
@@ -36,7 +36,7 @@
 import Data.Array.Comfort.Storable (Array)
 
 import qualified Data.Complex as Complex
-import Data.Complex (Complex)
+import Data.Complex (Complex((:+)))
 import Data.Tuple.HT (swap)
 
 import Prelude hiding (sum)
@@ -514,5 +514,17 @@
 getConstFunc _ = getConst
 
 
-data ComplexPart = RealPart | ImaginaryPart
-   deriving (Eq, Ord, Show, Enum, Bounded)
+{-
+ToDo:
+
+type ComplexShape =
+         Shape.NestedTuple Shape.TupleAccessor (Complex Shape.Element)
+
+This would allow the use of Complex.realPart as accessor,
+but it requires GHC>7.6.3 or so, where realPart has no RealFloat constraint.
+-}
+type ComplexShape = Shape.NestedTuple Shape.TupleIndex (Complex Shape.Element)
+
+ixReal, ixImaginary :: Shape.ElementIndex (Complex Shape.Element)
+ixReal :+ ixImaginary =
+   Shape.indexTupleFromShape (Shape.static :: ComplexShape)
diff --git a/src/Numeric/LAPACK/Scalar.hs b/src/Numeric/LAPACK/Scalar.hs
--- a/src/Numeric/LAPACK/Scalar.hs
+++ b/src/Numeric/LAPACK/Scalar.hs
@@ -1,6 +1,3 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE TypeOperators #-}
-{-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Scalar (
    RealOf,
    ComplexOf,
@@ -31,128 +28,4 @@
    conjugate,
    ) where
 
-import qualified Numeric.Netlib.Class as Class
-
-import Data.Functor.Identity (Identity(Identity, runIdentity))
-
-import qualified Data.Complex as Complex
-import Data.Complex (Complex((:+)))
-
-
-type family RealOf x
-
-type instance RealOf Float = Float
-type instance RealOf Double = Double
-type instance RealOf (Complex.Complex a) = a
-
-
-type ComplexOf x = Complex.Complex (RealOf x)
-
-
-data ComplexSingleton a where
-   Real :: (Class.Real a, RealOf a ~ a) => ComplexSingleton a
-   Complex :: (Class.Real a) => ComplexSingleton (Complex.Complex a)
-
-complexSingleton :: (Class.Floating a) => ComplexSingleton a
-complexSingleton = Class.switchFloating Real Real Complex Complex
-
-complexSingletonOf :: (Class.Floating a) => a -> ComplexSingleton a
-complexSingletonOf = const complexSingleton
-
-complexSingletonOfFunctor :: (Class.Floating a) => f a -> ComplexSingleton a
-complexSingletonOfFunctor = const complexSingleton
-
-withComplexSingleton :: (Class.Floating a) => (ComplexSingleton a -> a) -> a
-withComplexSingleton f = f complexSingleton
-
-
-data PrecisionSingleton a where
-   Float :: PrecisionSingleton Float
-   Double :: PrecisionSingleton Double
-
-precisionSingleton :: (Class.Real a) => PrecisionSingleton a
-precisionSingleton = Class.switchReal Float Double
-
-precisionOf :: (Class.Real a) => a -> PrecisionSingleton a
-precisionOf _ = precisionSingleton
-
-precisionOfFunctor :: (Class.Real a) => f a -> PrecisionSingleton a
-precisionOfFunctor _ = precisionSingleton
-
-
--- move to netlib-carray:Utility or netlib-ffi:Class
-zero, one, minusOne :: Class.Floating a => a
-zero = selectFloating 0 0 0 0
-one = selectFloating 1 1 1 1
-minusOne = selectFloating (-1) (-1) (-1) (-1)
-
-selectReal :: (Class.Real a) => Float -> Double -> a
-selectReal rf rd = runIdentity $ Class.switchReal (Identity rf) (Identity rd)
-
-selectFloating ::
-   (Class.Floating a) =>
-   Float -> Double -> Complex.Complex Float -> Complex.Complex Double -> a
-selectFloating rf rd cf cd =
-   withComplexSingleton $ \sw ->
-      case sw of
-         Real -> selectReal rf rd
-         Complex -> Class.switchReal cf cd
-
-
-
-equal :: (Class.Floating a) => a -> a -> Bool
-equal a b =
-   case complexSingletonOf a of
-      Real -> a==b
-      Complex -> a==b
-
-
-isZero :: (Class.Floating a) => a -> Bool
-isZero = equal zero
-
-
-fromReal :: (Class.Floating a) => RealOf a -> a
-fromReal a =
-   withComplexSingleton $ \sw ->
-      case sw of
-         Real -> a
-         Complex -> a:+0
-
-toComplex :: (Class.Floating a) => a -> ComplexOf a
-toComplex a =
-   case complexSingletonOf a of
-      Real -> a:+0
-      Complex -> a
-
-realPart :: (Class.Floating a) => a -> RealOf a
-realPart a =
-   case complexSingletonOf a of
-      Real -> a
-      Complex -> Complex.realPart a
-
-absolute :: (Class.Floating a) => a -> RealOf a
-absolute a =
-   case complexSingletonOf a of
-      Real -> abs a
-      Complex -> Complex.magnitude a
-
-
-norm1 :: (Class.Floating a) => a -> RealOf a
-norm1 a =
-   case complexSingletonOf a of
-      Real -> abs a
-      Complex -> case a of r:+i -> abs r + abs i
-
-
-absoluteSquared :: (Class.Floating a) => a -> RealOf a
-absoluteSquared a =
-   case complexSingletonOf a of
-      Real -> a*a
-      Complex -> case a of r:+i -> r*r+i*i
-
-
-conjugate :: (Class.Floating a) => a -> a
-conjugate a =
-   case complexSingletonOf a of
-      Real -> a
-      Complex -> Complex.conjugate a
+import Numeric.BLAS.Scalar
diff --git a/src/Numeric/LAPACK/Split.hs b/src/Numeric/LAPACK/Split.hs
--- a/src/Numeric/LAPACK/Split.hs
+++ b/src/Numeric/LAPACK/Split.hs
@@ -14,7 +14,7 @@
           swapOnRowMajor, sideSwapFromOrder,
           Triangle, uploFromOrder, flipOrder)
 import Numeric.LAPACK.Matrix.Extent.Private (Extent)
-import Numeric.LAPACK.Matrix.Modifier
+import Numeric.BLAS.Matrix.Modifier
          (Transposition, transposeOrder,
           Conjugation(NonConjugated, Conjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
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
@@ -4,245 +4,101 @@
    Vector,
    RealOf,
    ComplexOf,
-   toList,
-   fromList,
-   autoFromList,
-   CheckedArray.append, (+++),
+   Vector.toList,
+   Vector.fromList,
+   Vector.autoFromList,
+   CheckedArray.append, (Vector.+++),
    CheckedArray.take, CheckedArray.drop,
    CheckedArray.takeLeft, CheckedArray.takeRight,
-   swap,
+   Vector.swap,
    CheckedArray.singleton,
-   constant,
-   zero,
-   one,
-   unit,
-   dot, inner, (-*|),
-   sum,
-   absSum,
+   Vector.constant,
+   Vector.zero,
+   Vector.one,
+   Vector.unit,
+   Vector.dot, Vector.inner, (Vector.-*|),
+   Vector.sum,
+   Vector.absSum,
    norm1,
-   norm2,
-   norm2Squared,
+   Vector.norm2,
+   Vector.norm2Squared,
    normInf,
-   normInf1,
+   Vector.normInf1,
    argAbsMaximum,
-   argAbs1Maximum,
-   product,
-   scale, scaleReal, (.*|),
-   add, sub, (|+|), (|-|),
-   negate, raise,
-   mac,
-   mul, mulConj,
+   Vector.argAbs1Maximum,
+   Vector.product,
+   Vector.scale, Vector.scaleReal, (Vector..*|),
+   Vector.add, Vector.sub, (Vector.|+|), (Vector.|-|),
+   Vector.negate, Vector.raise,
+   Vector.mac,
+   Vector.mul, Vector.mulConj,
    divide, recip,
-   minimum, argMinimum,
-   maximum, argMaximum,
-   limits, argLimits,
+   Vector.minimum, Vector.argMinimum,
+   Vector.maximum, Vector.argMaximum,
+   Vector.limits, Vector.argLimits,
    CheckedArray.foldl,
    CheckedArray.foldl1,
    CheckedArray.foldMap,
 
    conjugate,
-   fromReal,
-   toComplex,
-   realPart,
-   imaginaryPart,
-   zipComplex,
-   unzipComplex,
+   Vector.fromReal,
+   Vector.toComplex,
+   Vector.realPart,
+   Vector.imaginaryPart,
+   Vector.zipComplex,
+   Vector.unzipComplex,
 
    random, RandomDistribution(..),
    ) where
 
-import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor
-import qualified Numeric.LAPACK.Vector.Private as Vector
+import qualified Numeric.LAPACK.Vector.Private as VectorPriv
 import qualified Numeric.LAPACK.Scalar as Scalar
 import qualified Numeric.LAPACK.Private as Private
-import Numeric.LAPACK.Matrix.Hermitian.Private
-         (Determinant(Determinant, getDeterminant))
-import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))
+import qualified Numeric.BLAS.Vector as Vector
 import Numeric.LAPACK.Linear.Private (withInfo)
-import Numeric.LAPACK.Scalar (ComplexOf, RealOf, minusOne, absolute)
-import Numeric.LAPACK.Private
-         (ComplexPart(RealPart, ImaginaryPart), fill, copyConjugate, realPtr)
+import Numeric.LAPACK.Scalar (ComplexOf, RealOf, absolute)
+import Numeric.LAPACK.Private (copyConjugate)
 
 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.Marshal.Array (advancePtr)
-import Foreign.ForeignPtr (withForeignPtr, castForeignPtr)
+import Foreign.ForeignPtr (withForeignPtr)
 import Foreign.Ptr (Ptr)
-import Foreign.Storable (Storable, peek, peekElemOff, pokeElemOff)
+import Foreign.Storable (Storable, 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.Monad.ST (runST)
 import Control.Monad (fmap, return, (=<<))
 import Control.Applicative (liftA3, (<$>))
 
-import qualified Data.Array.Comfort.Storable.Mutable.Unchecked as UMutArray
-import qualified Data.Array.Comfort.Storable.Mutable as MutArray
 import qualified Data.Array.Comfort.Storable.Unchecked as Array
 import qualified Data.Array.Comfort.Storable as CheckedArray
 import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Unchecked (Array(Array), append, (!))
-import Data.Array.Comfort.Shape ((::+))
+import Data.Array.Comfort.Storable.Unchecked (Array(Array))
 
-import Data.Function (id, flip, ($), (.))
+import Data.Function (id, ($), (.))
 import Data.Complex (Complex)
 import Data.Maybe (Maybe(Nothing,Just), maybe)
 import Data.Tuple.HT (mapFst, uncurry3)
-import Data.Tuple (fst, snd)
+import Data.Tuple (snd)
 import Data.Word (Word64)
 import Data.Bits (shiftR, (.&.))
-import Data.Ord (Ord, (>=))
+import Data.Ord (Ord)
 import Data.Eq (Eq, (==))
 import Data.Bool (Bool(False,True))
 
-import Prelude (Int, fromIntegral, (+), (-), (*), Char, Show, Enum, error, IO)
+import Prelude (Int, fromIntegral, (+), (-), (*), Show, Enum, error, IO)
 
 
 type Vector = Array
 
 
-toList :: (Shape.C sh, Storable a) => Vector sh a -> [a]
-toList = Array.toList
-
-fromList :: (Shape.C sh, Storable a) => sh -> [a] -> Vector sh a
-fromList = CheckedArray.fromList
-
-autoFromList :: (Storable a) => [a] -> Vector (Shape.ZeroBased Int) a
-autoFromList = Array.vectorFromList
-
-
-{- |
-> constant () = singleton
-
-However, singleton does not need 'Class.Floating' constraint.
--}
-constant :: (Shape.C sh, Class.Floating a) => sh -> a -> Vector sh a
-constant sh a = Array.unsafeCreateWithSize sh $ fill a
-
-zero :: (Shape.C sh, Class.Floating a) => sh -> Vector sh a
-zero = flip constant Scalar.zero
-
-one :: (Shape.C sh, Class.Floating a) => sh -> Vector sh a
-one = flip constant Scalar.one
-
-unit ::
-   (Shape.Indexed sh, Class.Floating a) =>
-   sh -> Shape.Index sh -> Vector sh a
-unit sh ix = Array.unsafeCreateWithSize sh $ \n xPtr -> do
-   fill Scalar.zero n xPtr
-   pokeElemOff xPtr (Shape.offset sh ix) Scalar.one
-
-
-{- |
-Precedence and associativity (right) of (List.++).
-This also matches '(::+)'.
--}
-infixr 5 +++
-
-(+++) ::
-   (Shape.C shx, Shape.C shy, Storable a) =>
-   Vector shx a -> Vector shy a -> Vector (shx::+shy) a
-(+++) = append
-
-
-swap ::
-   (Shape.Indexed sh, Storable a) =>
-   Shape.Index sh -> Shape.Index sh -> Vector sh a -> Vector sh a
-swap i j x =
-   runST (do
-      y <- MutArray.thaw x
-      xi <- MutArray.read y i
-      xj <- MutArray.read y j
-      MutArray.write y i xj
-      MutArray.write y j xi
-      UMutArray.unsafeFreeze y)
-
-
-infixl 7 -*|, .*|
-
-newtype Dot f a = Dot {runDot :: f a -> f a -> a}
-
-{- |
-> dot x y = Matrix.toScalar (singleRow x <#> singleColumn y)
--}
-dot, (-*|) ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   Vector sh a -> Vector sh a -> a
-(-*|) = dot
-dot =
-   runDot $
-   Class.switchFloating
-      (Dot dotReal)
-      (Dot dotReal)
-      (Dot $ dotComplex 'T')
-      (Dot $ dotComplex 'T')
-
-{- |
-> inner x y = dot (conjugate x) y
--}
-inner ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   Vector sh a -> Vector sh a -> a
-inner =
-   runDot $
-   Class.switchFloating
-      (Dot dotReal)
-      (Dot dotReal)
-      (Dot $ dotComplex 'C')
-      (Dot $ dotComplex 'C')
-
-dotReal ::
-   (Shape.C sh, Eq sh, Class.Real a) =>
-   Vector sh a -> Vector sh a -> a
-dotReal arrX@(Array shX _x) (Array shY y) = unsafePerformIO $ do
-   Call.assert "dot: shapes mismatch" (shX == shY)
-   evalContT $ do
-      (nPtr, sxPtr, incxPtr) <- vectorArgs arrX
-      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) =>
-   Char -> Vector sh (Complex a) -> Vector sh (Complex a) -> Complex a
-dotComplex trans (Array shX x) (Array shY y) = unsafePerformIO $ do
-   Call.assert "dot: shapes mismatch" (shX == shY)
-   evalContT $ do
-      let m = Shape.size shX
-      transPtr <- Call.char trans
-      mPtr <- Call.cint m
-      nPtr <- Call.cint 1
-      alphaPtr <- Call.number Scalar.one
-      xPtr <- ContT $ withForeignPtr x
-      ldxPtr <- Call.leadingDim m
-      yPtr <- ContT $ withForeignPtr y
-      incyPtr <- Call.cint 1
-      betaPtr <- Call.number Scalar.zero
-      zPtr <- Call.alloca
-      inczPtr <- Call.cint 1
-      liftIO $
-         Private.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.Floating a) => Vector sh a -> RealOf a
 norm1 arr = unsafePerformIO $
    evalContT $ liftIO . uncurry3 csum1 =<< vectorArgs arr
@@ -256,78 +112,17 @@
       (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'.
--}
-absSum :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
-absSum arr = unsafePerformIO $
-   evalContT $ liftIO . uncurry3 asum =<< vectorArgs arr
-
-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 arr = unsafePerformIO $
-   evalContT $ liftIO . uncurry3 nrm2 =<< vectorArgs arr
-
-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)}
 
 
-norm2Squared :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
-norm2Squared =
-   getDeterminant $
-   Class.switchFloating
-      (Determinant norm2SquaredReal)
-      (Determinant norm2SquaredReal)
-      (Determinant $ norm2SquaredReal . decomplex)
-      (Determinant $ norm2SquaredReal . decomplex)
-
-norm2SquaredReal :: (Shape.C sh, Class.Real a) => Vector sh a -> a
-norm2SquaredReal arr =
-   unsafePerformIO $ evalContT $ do
-      (nPtr, sxPtr, incxPtr) <- vectorArgs arr
-      liftIO $ BlasReal.dot nPtr sxPtr incxPtr sxPtr incxPtr
-
-
 normInf :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
 normInf arr = unsafePerformIO $
    evalContT $ do
       (nPtr, sxPtr, incxPtr) <- vectorArgs arr
       liftIO $
          fmap (absolute . maybe Scalar.zero snd) $
-         peekElemOff1 sxPtr =<< Vector.absMax nPtr sxPtr incxPtr
-
-{- |
-Computes (almost) the infinity norm of the vector.
-For complex numbers every element is replaced
-by the sum of the absolute component values first.
--}
-normInf1 :: (Shape.C sh, Class.Floating a) => Vector sh a -> RealOf a
-normInf1 arr = unsafePerformIO $
-   evalContT $ do
-      (nPtr, sxPtr, incxPtr) <- vectorArgs arr
-      liftIO $
-         fmap (Scalar.norm1 . maybe Scalar.zero snd) $
-         peekElemOff1 sxPtr =<< BlasGen.iamax nPtr sxPtr incxPtr
-
+         peekElemOff1 sxPtr =<< VectorPriv.absMax nPtr sxPtr incxPtr
 
 {- |
 Returns the index and value of the element with the maximal absolute value.
@@ -344,28 +139,9 @@
             (maybe
                (error "Vector.argAbsMaximum: empty vector")
                (mapFst (Shape.uncheckedIndexFromOffset $ Array.shape arr))) $
-         peekElemOff1 sxPtr =<< Vector.absMax nPtr sxPtr incxPtr
+         peekElemOff1 sxPtr =<< VectorPriv.absMax nPtr sxPtr incxPtr
 
 
-{- |
-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.InvIndexed sh, Class.Floating a) =>
-   Vector sh a -> (Shape.Index sh, a)
-argAbs1Maximum arr = unsafePerformIO $
-   evalContT $ do
-      (nPtr, sxPtr, incxPtr) <- vectorArgs arr
-      liftIO $
-         fmap
-            (maybe
-               (error "Vector.argAbs1Maximum: empty vector")
-               (mapFst (Shape.uncheckedIndexFromOffset $ Array.shape arr))) $
-         peekElemOff1 sxPtr =<< BlasGen.iamax nPtr sxPtr incxPtr
-
 vectorArgs ::
    (Shape.C sh) => Array sh a -> ContT r IO (Ptr CInt, Ptr a, Ptr CInt)
 vectorArgs (Array sh x) =
@@ -383,184 +159,6 @@
          else Just . (,) ki <$> peekElemOff ptr ki
 
 
-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
-
-
-{- |
-For restrictions see 'limits'.
--}
-minimum, maximum :: (Shape.C sh, Class.Real a) => Vector sh a -> a
-minimum = fst . limits
-maximum = snd . limits
-
-{- |
-For restrictions see 'limits'.
--}
-argMinimum, argMaximum ::
-   (Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>
-   Vector sh a -> (ix,a)
-argMinimum = fst . argLimits
-argMaximum = snd . argLimits
-
-{- |
-It should hold @limits x = Array.limits x@.
-The function is based on fast BLAS functions.
-It should be faster than @Array.minimum@ and @Array.maximum@
-although it is certainly not as fast as possible.
-It is less precise if minimum and maximum differ considerably in magnitude
-and there are several minimum or maximum candidates of similar value.
-E.g. you cannot rely on the property
-that @raise (- minimum x) x@ has only non-negative elements.
--}
-limits :: (Shape.C sh, Class.Real a) => Vector sh a -> (a,a)
-limits xs0 =
-   let xs = Array.mapShape Shape.Deferred xs0
-       x0 = snd $ argAbsMaximum xs
-       x1 = xs ! fst (argAbsMaximum (raise (-x0) xs))
-   in if x0>=0 then (x1,x0) else (x0,x1)
-
-argLimits ::
-   (Shape.InvIndexed sh, Shape.Index sh ~ ix, Class.Real a) =>
-   Vector sh a -> ((ix,a),(ix,a))
-argLimits xs =
-   let p0@(_i0,x0) = argAbsMaximum xs
-       p1 = (i1,xs!i1); i1 = fst $ argAbsMaximum $ raise (-x0) xs
-   in if x0>=0 then (p1,p0) else (p0,p1)
-
-
-scale, _scale, (.*|) ::
-   (Shape.C sh, Class.Floating a) =>
-   a -> Vector sh a -> Vector sh a
-(.*|) = scale
-
-scale alpha (Array sh x) = Array.unsafeCreateWithSize sh $ \n syPtr -> do
-   evalContT $ do
-      alphaPtr <- Call.number alpha
-      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.unsafeCreateWithSize sh $ \n cPtr -> do
-   let m = 1
-   let k = 1
-   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 Scalar.one
-      aPtr <- Call.number a
-      ldaPtr <- Call.leadingDim m
-      bPtr <- ContT $ withForeignPtr b
-      ldbPtr <- Call.leadingDim k
-      betaPtr <- Call.number Scalar.zero
-      ldcPtr <- Call.leadingDim m
-      liftIO $
-         BlasGen.gemm
-            transaPtr transbPtr mPtr nPtr kPtr alphaPtr
-            aPtr ldaPtr bPtr ldbPtr betaPtr cPtr ldcPtr
-
-
-scaleReal ::
-   (Shape.C sh, Class.Floating a) =>
-   RealOf a -> Vector sh a -> Vector sh a
-scaleReal =
-   getScaleReal $
-   Class.switchFloating
-      (ScaleReal scale)
-      (ScaleReal scale)
-      (ScaleReal $ \x -> recomplex . scale x . decomplex)
-      (ScaleReal $ \x -> recomplex . scale x . decomplex)
-
-newtype ScaleReal f a = ScaleReal {getScaleReal :: RealOf a -> f a -> f a}
-
-
-decomplex ::
-   (Class.Real a) =>
-   Vector sh (Complex a) -> Vector (sh, Shape.Enumeration ComplexPart) a
-decomplex (Array sh a) = Array (sh, Shape.Enumeration) (castForeignPtr a)
-
-recomplex ::
-   (Class.Real a) =>
-   Vector (sh, Shape.Enumeration ComplexPart) a -> Vector sh (Complex a)
-recomplex (Array (sh, Shape.Enumeration) a) = Array sh (castForeignPtr a)
-
-
-
-infixl 6 |+|, |-|, `add`, `sub`
-
-
-add, sub, (|+|), (|-|) ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   Vector sh a -> Vector sh a -> Vector sh a
-add = mac Scalar.one
-sub x y = mac minusOne y x
-
-(|+|) = add
-(|-|) = sub
-
-mac ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   a -> Vector sh a -> Vector sh a -> Vector sh a
-mac alpha x y =
-   if Array.shape x == Array.shape y
-      then Vector.mac alpha x y
-      else error "mac: shapes mismatch"
-
-negate :: (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh a
-negate =
-   getConjugate $
-   Class.switchFloating
-      (Conjugate $ scaleReal Scalar.minusOne)
-      (Conjugate $ scaleReal Scalar.minusOne)
-      (Conjugate $ scaleReal Scalar.minusOne)
-      (Conjugate $ scaleReal Scalar.minusOne)
-
-
-raise :: (Shape.C sh, Class.Floating a) => a -> Array sh a -> Array sh a
-raise c (Array shX x) =
-   Array.unsafeCreateWithSize shX $ \n yPtr -> evalContT $ do
-      nPtr <- Call.cint n
-      cPtr <- Call.number c
-      onePtr <- Call.number Scalar.one
-      inccPtr <- Call.cint 0
-      xPtr <- ContT $ withForeignPtr x
-      inc1Ptr <- Call.cint 1
-      liftIO $ do
-         BlasGen.copy nPtr xPtr inc1Ptr yPtr inc1Ptr
-         BlasGen.axpy nPtr onePtr cPtr inccPtr yPtr inc1Ptr
-
-
-mul ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   Vector sh a -> Vector sh a -> Vector sh a
-mul = mulConjugation NonConjugated
-
-{- |
-prop> \xs ys -> mulConj xs ys == mul (conjugate xs) ys
--}
-mulConj ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   Vector sh a -> Vector sh a -> Vector sh a
-mulConj = mulConjugation Conjugated
-
-mulConjugation ::
-   (Shape.C sh, Eq sh, Class.Floating a) =>
-   Conjugation -> Vector sh a -> Vector sh a -> Vector sh a
-mulConjugation conj (Array shA a) (Array shX x) =
-      Array.unsafeCreateWithSize shX $ \n yPtr -> do
-   Call.assert "mul: shapes mismatch" (shA == shX)
-   evalContT $ do
-      aPtr <- ContT $ withForeignPtr a
-      xPtr <- ContT $ withForeignPtr x
-      liftIO $ Private.mul conj n aPtr 1 xPtr 1 yPtr 1
-
 divide ::
    (Shape.C sh, Eq sh, Class.Floating a) =>
    Vector sh a -> Vector sh a -> Vector sh a
@@ -585,8 +183,10 @@
 
 recip :: (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh a
 recip x =
-   Vector.recheck $
-   divide (Vector.uncheck $ one $ Array.shape x) (Vector.uncheck x)
+   VectorPriv.recheck $
+   divide
+      (VectorPriv.uncheck $ Vector.one $ Array.shape x)
+      (VectorPriv.uncheck x)
 
 
 newtype Conjugate f a = Conjugate {getConjugate :: f a -> f a}
@@ -612,86 +212,6 @@
       incxPtr <- Call.cint 1
       incyPtr <- Call.cint 1
       liftIO $ copyConjugate nPtr sxPtr incxPtr syPtr incyPtr
-
-
-fromReal ::
-   (Shape.C sh, Class.Floating a) => Vector sh (RealOf a) -> Vector sh a
-fromReal =
-   getFromReal $
-   Class.switchFloating
-      (FromReal id)
-      (FromReal id)
-      (FromReal complexFromReal)
-      (FromReal complexFromReal)
-
-newtype FromReal f a = FromReal {getFromReal :: f (RealOf a) -> f a}
-
-toComplex ::
-   (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh (ComplexOf a)
-toComplex =
-   getToComplex $
-   Class.switchFloating
-      (ToComplex complexFromReal)
-      (ToComplex complexFromReal)
-      (ToComplex id)
-      (ToComplex id)
-
-newtype ToComplex f a = ToComplex {getToComplex :: f a -> f (ComplexOf a)}
-
-complexFromReal ::
-   (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh (Complex a)
-complexFromReal (Array sh x) =
-   Array.unsafeCreateWithSize sh $ \n yPtr ->
-   case realPtr yPtr of
-      yrPtr -> evalContT $ do
-         nPtr <- Call.cint n
-         xPtr <- ContT $ withForeignPtr x
-         incxPtr <- Call.cint 1
-         incyPtr <- Call.cint 2
-         inczPtr <- Call.cint 0
-         zPtr <- Call.number Scalar.zero
-         liftIO $ do
-            BlasGen.copy nPtr xPtr incxPtr yrPtr incyPtr
-            BlasGen.copy nPtr zPtr inczPtr (advancePtr yrPtr 1) incyPtr
-
-
-realPart ::
-   (Shape.C sh, Class.Floating a) => Vector sh a -> Vector sh (RealOf a)
-realPart =
-   getToReal $
-   Class.switchFloating
-      (ToReal id)
-      (ToReal id)
-      (ToReal $ RowMajor.takeColumn RealPart . decomplex)
-      (ToReal $ RowMajor.takeColumn RealPart . decomplex)
-
-newtype ToReal f a = ToReal {getToReal :: f a -> f (RealOf a)}
-
-imaginaryPart ::
-   (Shape.C sh, Class.Real a) => Vector sh (Complex a) -> Vector sh a
-imaginaryPart = RowMajor.takeColumn ImaginaryPart . decomplex
-
-
-zipComplex ::
-   (Shape.C sh, Eq sh, Class.Real a) =>
-   Vector sh a -> Vector sh a -> Vector sh (Complex a)
-zipComplex (Array shr xr) (Array shi xi) =
-   Array.unsafeCreateWithSize shr $ \n yPtr -> evalContT $ do
-      liftIO $ Call.assert "zipComplex: shapes mismatch" (shr==shi)
-      nPtr <- Call.cint n
-      xrPtr <- ContT $ withForeignPtr xr
-      xiPtr <- ContT $ withForeignPtr xi
-      let yrPtr = realPtr yPtr
-      incxPtr <- Call.cint 1
-      incyPtr <- Call.cint 2
-      liftIO $ do
-         BlasGen.copy nPtr xrPtr incxPtr yrPtr incyPtr
-         BlasGen.copy nPtr xiPtr incxPtr (advancePtr yrPtr 1) incyPtr
-
-unzipComplex ::
-   (Shape.C sh, Class.Real a) =>
-   Vector sh (Complex a) -> (Vector sh a, Vector sh a)
-unzipComplex x = (realPart x, imaginaryPart x)
 
 
 data RandomDistribution =
diff --git a/src/Numeric/LAPACK/Vector/Private.hs b/src/Numeric/LAPACK/Vector/Private.hs
--- a/src/Numeric/LAPACK/Vector/Private.hs
+++ b/src/Numeric/LAPACK/Vector/Private.hs
@@ -4,35 +4,13 @@
 
 import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
 import qualified Numeric.BLAS.FFI.Generic as BlasGen
-import qualified Numeric.Netlib.Utility as Call
 import qualified Numeric.Netlib.Class as Class
 
-import Foreign.ForeignPtr (withForeignPtr)
 import Foreign.Ptr (Ptr)
 import Foreign.C.Types (CInt)
 
-import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
-import Control.Monad.IO.Class (liftIO)
-
 import qualified Data.Array.Comfort.Storable.Unchecked as Array
-import qualified Data.Array.Comfort.Shape as Shape
-import Data.Array.Comfort.Storable.Unchecked (Array(Array))
-
-
-mac ::
-   (Shape.C sh, Class.Floating a) =>
-   a -> Array sh a -> Array sh a -> Array sh a
-mac alpha (Array shX x) (Array _shY y) =
-   Array.unsafeCreateWithSize shX $ \n szPtr -> evalContT $ do
-      nPtr <- Call.cint n
-      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
+import Data.Array.Comfort.Storable.Unchecked (Array)
 
 
 newtype ArgMaximum a =