diff --git a/Makefile b/Makefile
--- a/Makefile
+++ b/Makefile
@@ -9,8 +9,11 @@
 run-test-cabal:	update-test
 	cabal test --show-details=streaming
 
-update-test:
-	doctest-extract-0.1 -i src/ -o test/ --module-prefix DocTest --library-main=Main $$(cat test-module.list)
+update-test:	test-module.list
+	doctest-extract-0.1 -i src/ -o test/ --module-prefix DocTest --library-main=Main $$(cat $<)
+
+test-module.list:	lapack.cabal
+	grep '^ \+DocTest\.Numeric\.LAPACK\.' $< | cut -d. -f2- >$@
 
 %.html: %.md
 	pandoc --standalone --metadata pagetitle="$*" --output=$@ $<
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.1.1
+Version:          0.5.2
 License:          BSD-3-Clause
 License-File:     LICENSE
 Author:           Henning Thielemann <haskell@henning-thielemann.de>
@@ -57,7 +57,7 @@
   test-module.list
 
 Source-Repository this
-  Tag:         0.5.1.1
+  Tag:         0.5.2
   Type:        darcs
   Location:    https://hub.darcs.net/thielema/lapack/
 
@@ -71,10 +71,10 @@
 
 Library
   Build-Depends:
-    comfort-blas >=0.0 && <0.1,
+    comfort-blas >=0.0.3 && <0.1,
     lapack-ffi >=0.0.3 && <0.1,
     blas-ffi >=0.0 && <0.2,
-    netlib-ffi >=0.1.1 && <0.2,
+    netlib-ffi >=0.1.2 && <0.2,
     comfort-array-shape >=0.0 && <0.1,
     comfort-array >=0.5.2 && <0.6,
     guarded-allocation >=0.0 && <0.1,
@@ -86,6 +86,7 @@
     lazyio >=0.1 && <0.2,
     transformers >=0.4 && <0.7,
     tfp >=1.0.2 && <1.1,
+    containers >=0.4 && <0.8,
     semigroups >=0.18.3 && <1.0,
     non-empty >=0.3 && <0.4,
     Stream >=0.4.7 && <0.5,
@@ -166,6 +167,7 @@
     Numeric.LAPACK.Matrix.BandedHermitian.Basic
     Numeric.LAPACK.Matrix.BandedHermitian.Eigen
     Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite.Linear
+    Numeric.LAPACK.Matrix.Layout.Basic
     Numeric.LAPACK.Matrix.Layout.Private
     Numeric.LAPACK.Matrix.Extent.Private
     Numeric.LAPACK.Matrix.Extent.Strict
@@ -191,7 +193,6 @@
     Numeric.LAPACK.Matrix.Array.Private
     Numeric.LAPACK.Matrix.Lazy.UpperTriangular
     Numeric.LAPACK.Matrix.ModifierTyped
-    Numeric.LAPACK.Matrix.RowMajor
     Numeric.LAPACK.Matrix.Private
     Numeric.LAPACK.Vector.Private
     Numeric.LAPACK.Private
diff --git a/src/Numeric/LAPACK/Format.hs b/src/Numeric/LAPACK/Format.hs
--- a/src/Numeric/LAPACK/Format.hs
+++ b/src/Numeric/LAPACK/Format.hs
@@ -79,7 +79,7 @@
    format cfg (a,b,c) = format cfg a /+/ format cfg b /+/ format cfg c
 
 instance (Shape.C sh) => Format (Perm.Permutation sh) where
-   format _cfg = Perm.format
+   format = Perm.format
 
 instance (FormatArray sh, Class.Floating a) => Format (Array sh a) where
    format = formatArray
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.BLAS.Matrix.Modifier (Transposition, Conjugation, Inversion)
+import Numeric.Netlib.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.BLAS.Matrix.Modifier
+import Numeric.Netlib.Modifier
          (Transposition(NonTransposed, Transposed),
           Conjugation(NonConjugated, Conjugated),
           Inversion(NonInverted, Inverted))
@@ -127,7 +127,7 @@
 instance Matrix.Format LU where
    type FormatExtra LU extra = extra ~ ()
    format fmt lu@(LowerUpper _ipiv m) =
-      Perm.format (extractP NonInverted lu)
+      Perm.format fmt (extractP NonInverted lu)
       /+/
       formatArray fmt m
 
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
@@ -59,6 +59,8 @@
    takeTop, takeBottom,
    takeLeft, takeRight,
    takeRowArray, takeColumnArray,
+   takeRowSet, takeColumnSet,
+   takeRowIntSet, takeColumnIntSet,
    swapRows, swapColumns,
    reverseRows, reverseColumns,
    fromRowMajor, toRowMajor,
@@ -136,7 +138,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.BLAS.Matrix.Modifier as Mod
+import qualified Numeric.Netlib.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
@@ -165,6 +167,8 @@
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Either.HT as EitherHT
 import Data.Function.HT (Id)
+import Data.IntSet (IntSet)
+import Data.Set (Set)
 
 import Prelude hiding (map)
 
@@ -497,11 +501,48 @@
    General height width a -> General sh width a
 takeRowArray = ArrMatrix.lift1 . Plain.takeRowArray
 
+takeRowSet ::
+   (Extent.C vert,
+    Shape.Indexed height, Shape.Index height ~ ix, Ord ix,
+    Shape.C width, Class.Floating a) =>
+   Set ix ->
+   Full Extent.Size vert Extent.Big height width a ->
+   Full Extent.Size vert Extent.Big (Set ix) width a
+takeRowSet = ArrMatrix.lift1 . Plain.takeRowSet
+
+takeRowIntSet ::
+   (Extent.C vert,
+    Shape.Indexed height, Shape.Index height ~ Int,
+    Shape.C width, Class.Floating a) =>
+   IntSet ->
+   Full Extent.Size vert Extent.Big height width a ->
+   Full Extent.Size vert Extent.Big IntSet width a
+takeRowIntSet = ArrMatrix.lift1 . Plain.takeRowIntSet
+
+
 takeColumnArray ::
    (Shape.C height, Shape.Indexed width, Shape.C sh, Class.Floating a) =>
    BoxedArray.Array sh (Shape.Index width) ->
    General height width a -> General height sh a
 takeColumnArray = ArrMatrix.lift1 . Plain.takeColumnArray
+
+takeColumnSet ::
+   (Extent.C horiz,
+    Shape.Indexed width, Shape.Index width ~ ix, Ord ix,
+    Shape.C height, Class.Floating a) =>
+   Set ix ->
+   Full Extent.Size Extent.Big horiz height width a ->
+   Full Extent.Size Extent.Big horiz height (Set ix) a
+takeColumnSet = ArrMatrix.lift1 . Plain.takeColumnSet
+
+takeColumnIntSet ::
+   (Extent.C horiz,
+    Shape.Indexed width, Shape.Index width ~ Int,
+    Shape.C height, Class.Floating a) =>
+   IntSet ->
+   Full Extent.Size Extent.Big horiz height width a ->
+   Full Extent.Size Extent.Big horiz height IntSet a
+takeColumnIntSet = ArrMatrix.lift1 . Plain.takeColumnIntSet
 
 
 
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.BLAS.Matrix.Modifier (Transposition(NonTransposed,Transposed))
+import Numeric.Netlib.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
@@ -38,15 +38,15 @@
 import qualified Numeric.LAPACK.Matrix.Triangular.Basic as Triangular
 import qualified Numeric.LAPACK.Matrix.Mosaic.Generic as Mosaic
 import qualified Numeric.LAPACK.Matrix.Mosaic.Private as Mos
-import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor
 import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Private as Private
+import qualified Numeric.BLAS.Matrix.RowMajor as RowMajor
 import qualified Data.Array.Comfort.Shape.Static as ShapeStatic
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor,ColumnMajor), transposeFromOrder, swapOnRowMajor,
           UnaryProxy, addOffDiagonals)
-import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.Netlib.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.BLAS.Matrix.Modifier (Transposition(NonTransposed, Transposed))
+import Numeric.Netlib.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
@@ -27,19 +27,19 @@
 import qualified Numeric.LAPACK.Matrix.Extent as Extent
 import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
 import qualified Numeric.LAPACK.Matrix.Mosaic.Private as Mos
-import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor
 import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Vector.Private as VectorPriv
 import qualified Numeric.LAPACK.Vector as Vector
+import qualified Numeric.BLAS.Matrix.RowMajor as RowMajor
 import qualified Data.Array.Comfort.Shape.Static as ShapeStatic
 import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
 import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)
 import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor,ColumnMajor), uploFromOrder,
+         (Order(RowMajor,ColumnMajor), uploFromOrder, conjugatedOnRowMajor,
           UnaryProxy, natFromProxy)
-import Numeric.BLAS.Matrix.Modifier
+import Numeric.Netlib.Modifier
          (Transposition(NonTransposed, Transposed), transposeOrder,
-          Conjugation(NonConjugated), conjugatedOnRowMajor)
+          Conjugation(NonConjugated))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf, zero, one)
 import Numeric.LAPACK.Private
@@ -74,6 +74,7 @@
 
 import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
 import Control.Monad.IO.Class (liftIO)
+import Control.Monad (when)
 
 import Data.Foldable (for_)
 import Data.Tuple.HT (mapPair)
@@ -192,10 +193,10 @@
                   (n-1) k
                   lda (advancePtr aPtr (k+1))
                   ldb (advancePtr bPtr (k+1))
-               fill zero k (advancePtr bPtr (ldb*n-k))
+               when (n>0) $ fill zero k (advancePtr bPtr (ldb*n-k))
             RowMajor -> do
                copySubMatrix lda n lda aPtr ldb (advancePtr bPtr k)
-               fill zero k bPtr
+               when (n>0) $ fill zero k bPtr
                rowToColumnMajor copyConjugate
                   (n-1) k
                   lda (advancePtr aPtr 1)
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
@@ -11,8 +11,8 @@
 import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Shape as ExtShape
 import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
-import Numeric.LAPACK.Matrix.Layout.Private (uploFromOrder)
-import Numeric.BLAS.Matrix.Modifier (conjugatedOnRowMajor)
+import Numeric.LAPACK.Matrix.Layout.Private
+         (uploFromOrder, 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.BLAS.Matrix.Modifier (Conjugation(Conjugated))
+import Numeric.Netlib.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
@@ -4,12 +4,12 @@
 
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
-import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Private as Private
+import qualified Numeric.BLAS.Matrix.RowMajor as RowMajor
 import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor, ColumnMajor), transposeFromOrder, flipOrder)
-import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
+         (Order(RowMajor, ColumnMajor), transposeFromOrder)
+import Numeric.Netlib.Modifier (Conjugation(NonConjugated))
 import Numeric.LAPACK.Matrix.Private
          (Full, Tall, Wide, Square, General, fromFull, ShapeInt, revealOrder)
 import Numeric.LAPACK.Vector (Vector)
@@ -37,6 +37,11 @@
 import Data.Complex (Complex)
 
 
+
+fullExtent ::
+   Full meas vert horiz height width a -> Extent meas vert horiz height width
+fullExtent = Layout.fullExtent . Array.shape
+
 caseTallWide ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.C height, Shape.C width) =>
@@ -366,11 +371,7 @@
    General heightB widthA a -> Full meas vert horiz heightB widthB a ->
    Full meas vert horiz (heightA::+heightB) (widthA::+widthB) a
 stackBiased vertBias horizBias a b c d =
-   mapExtent
-      (\ _ ->
-         Extent.stack
-            (Layout.fullExtent $ Array.shape a)
-            (Layout.fullExtent $ Array.shape d)) $
+   mapExtent (\ _ -> Extent.stack (fullExtent a) (fullExtent d)) $
    above vertBias Extent.appendAny
       (beside horizBias Extent.appendAny (fromFull a) b)
       (beside horizBias Extent.appendAny c (fromFull d))
@@ -511,37 +512,47 @@
    Full meas vert horiz fuse width a ->
    Full meas vert horiz height width a
 -- preserve order of the right factor
-multiply
-   (Array (Layout.Full orderA extentA) a)
-   (Array (Layout.Full orderB extentB) b) =
-   case Extent.fuse extentA extentB of
-      Nothing -> error "multiply: fuse shapes mismatch"
-      Just extent ->
-         Array.unsafeCreate (Layout.Full orderB extent) $ \cPtr -> do
-
-      let (height,fuse) = Extent.dimensions extentA
-      let width = Extent.width extentB
-      let m = Shape.size height
-      let n = Shape.size width
-      let k = Shape.size fuse
-      case orderB of
-         RowMajor ->
-            Private.multiplyMatrix (flipOrder orderB) (flipOrder orderA)
-               n k m b a cPtr
-         ColumnMajor -> Private.multiplyMatrix orderA orderB m k n a b cPtr
+multiply a b =
+   fullFromTransposable (Extent.fuseUnchecked (fullExtent a) (fullExtent b)) $
+   let at = transposableFromFull a in
+   let bt = transposableFromFull b in
+   case bt of
+      RowMajor.NonTransposed _ ->
+         RowMajor.NonTransposed $
+         RowMajor.multiplyTransposable at bt
+      RowMajor.Transposed _ ->
+         RowMajor.Transposed $
+         RowMajor.multiplyTransposable
+            (RowMajor.transposeTransposable bt)
+            (RowMajor.transposeTransposable at)
 
 -- always return ColumnMajor
-multiplyColumnMajor
-   (Array (Layout.Full orderA extentA) a)
-   (Array (Layout.Full orderB extentB) b) =
-   case Extent.fuse extentA extentB of
-      Nothing -> error "multiply: fuse shapes mismatch"
-      Just extent ->
-         Array.unsafeCreate (Layout.Full ColumnMajor extent) $ \cPtr -> do
+multiplyColumnMajor a b =
+   fullFromTransposable (Extent.fuseUnchecked (fullExtent a) (fullExtent b)) $
+   RowMajor.Transposed $
+   RowMajor.multiplyTransposable
+      (RowMajor.transposeTransposable $ transposableFromFull b)
+      (RowMajor.transposeTransposable $ transposableFromFull a)
 
-      let (height,fuse) = Extent.dimensions extentA
-      let width = Extent.width extentB
-      let m = Shape.size height
-      let n = Shape.size width
-      let k = Shape.size fuse
-      Private.multiplyMatrix orderA orderB m k n a b cPtr
+
+transposableFromFull ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   Full meas vert horiz height width a ->
+   RowMajor.Transposable height width a
+transposableFromFull (Array (Layout.Full order extent) a) =
+   let (height,width) = Extent.dimensions extent in
+   case order of
+      RowMajor -> RowMajor.NonTransposed $ Array (height,width) a
+      ColumnMajor -> RowMajor.Transposed $ Array (width,height) a
+
+fullFromTransposable ::
+   (Extent.Measure meas, Extent.C vert, Extent.C horiz) =>
+   Extent meas vert horiz height width ->
+   RowMajor.Transposable height width a ->
+   Full meas vert horiz height width a
+fullFromTransposable extent at =
+   case at of
+      RowMajor.NonTransposed a ->
+         Array.reshape (Layout.Full RowMajor extent) a
+      RowMajor.Transposed a ->
+         Array.reshape (Layout.Full ColumnMajor extent) a
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.BLAS.Matrix.Modifier
+import Numeric.Netlib.Modifier
          (Transposition(NonTransposed,Transposed),
           Inversion(Inverted))
 import Numeric.LAPACK.Vector (Vector)
diff --git a/src/Numeric/LAPACK/Matrix/Extent/Private.hs b/src/Numeric/LAPACK/Matrix/Extent/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Extent/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Extent/Private.hs
@@ -482,6 +482,14 @@
 fuse (Square s0) (Square s1) = toMaybe (s0==s1) $ Square s0
 fuse (Separate h f0) (Separate f1 w) = toMaybe (f0==f1) $ Separate h w
 
+fuseUnchecked ::
+   (Measure meas, C vert, C horiz) =>
+   Extent meas vert horiz height fuse ->
+   Extent meas vert horiz fuse width ->
+   Extent meas vert horiz height width
+fuseUnchecked (Square s0) (Square _s1) = Square s0
+fuseUnchecked (Separate h _f0) (Separate _f1 w) = Separate h w
+
 relaxMeasureWith ::
    (Measure measA, Measure measB,
     MultiplyMeasure measA measB ~ measC,
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.BLAS.Matrix.Modifier (Transposition(NonTransposed, Transposed))
+import Numeric.Netlib.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
@@ -21,10 +21,9 @@
           rowMajorPointers, columnMajorPointers,
           withPacking, noLabel, applyFuncPair, triArg)
 import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor,ColumnMajor), uploFromOrder)
-import Numeric.BLAS.Matrix.Modifier
-         (Transposition(NonTransposed, Transposed),
-          Conjugation(Conjugated), conjugatedOnRowMajor)
+         (Order(RowMajor,ColumnMajor), uploFromOrder, conjugatedOnRowMajor)
+import Numeric.Netlib.Modifier
+         (Transposition(NonTransposed, Transposed), Conjugation(Conjugated))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (RealOf, zero)
 import Numeric.LAPACK.Private (fill, realPtr, condConjugate)
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
@@ -11,8 +11,8 @@
 import qualified Numeric.LAPACK.Shape as ExtShape
 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.BLAS.Matrix.Modifier (conjugatedOnRowMajor)
+import Numeric.LAPACK.Matrix.Layout.Private
+         (Order(ColumnMajor), uploFromOrder, 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.BLAS.Matrix.Modifier (Conjugation(Conjugated))
+import Numeric.Netlib.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.hs b/src/Numeric/LAPACK/Matrix/Layout.hs
--- a/src/Numeric/LAPACK/Matrix/Layout.hs
+++ b/src/Numeric/LAPACK/Matrix/Layout.hs
@@ -2,7 +2,7 @@
    General,
    Tall,
    Wide,
-   Square,
+   Square, squareSize,
    Full(..), fullHeight, fullWidth,
    Order(..), flipOrder,
    general,
diff --git a/src/Numeric/LAPACK/Matrix/Layout/Basic.hs b/src/Numeric/LAPACK/Matrix/Layout/Basic.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/LAPACK/Matrix/Layout/Basic.hs
@@ -0,0 +1,32 @@
+module Numeric.LAPACK.Matrix.Layout.Basic where
+
+import qualified Numeric.Netlib.Modifier as Modi
+import Numeric.Netlib.Layout (Order(RowMajor,ColumnMajor))
+
+import Data.Tuple.HT (swap)
+
+
+transposeFromOrder :: Order -> Char
+transposeFromOrder RowMajor = 'T'
+transposeFromOrder ColumnMajor = 'N'
+
+swapOnRowMajor :: Order -> (a,a) -> (a,a)
+swapOnRowMajor order =
+   case order of
+      RowMajor -> swap
+      ColumnMajor -> id
+
+sideSwapFromOrder :: Order -> (a,a) -> (Char, (a,a))
+sideSwapFromOrder order (m0,n0) =
+   let ((side,m), (_,n)) = swapOnRowMajor order (('L', m0), ('R', n0))
+   in (side,(m,n))
+
+
+conjugatedOnRowMajor :: Order -> Modi.Conjugation
+conjugatedOnRowMajor RowMajor = Modi.Conjugated
+conjugatedOnRowMajor ColumnMajor = Modi.NonConjugated
+
+
+uploFromOrder :: Order -> Char
+uploFromOrder RowMajor = 'L'
+uploFromOrder ColumnMajor = 'U'
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
@@ -5,12 +5,16 @@
 {-# LANGUAGE GADTs #-}
 module Numeric.LAPACK.Matrix.Layout.Private (
    module Numeric.LAPACK.Matrix.Layout.Private,
-   module Numeric.BLAS.Matrix.Layout,
+   module Numeric.LAPACK.Matrix.Layout.Basic,
+   Order(..), flipOrder,
    ) 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 Numeric.LAPACK.Matrix.Layout.Basic
+         (transposeFromOrder, swapOnRowMajor, uploFromOrder,
+          sideSwapFromOrder, conjugatedOnRowMajor)
+import Numeric.Netlib.Layout (Order(..), flipOrder)
 
 import qualified Type.Data.Num.Unary.Literal as TypeNum
 import qualified Type.Data.Num.Unary.Proof as Proof
@@ -30,18 +34,6 @@
 import Data.Bool.HT (if')
 
 
-swapOnRowMajor :: Order -> (a,a) -> (a,a)
-swapOnRowMajor order =
-   case order of
-      RowMajor -> swap
-      ColumnMajor -> id
-
-sideSwapFromOrder :: Order -> (a,a) -> (Char, (a,a))
-sideSwapFromOrder order (m0,n0) =
-   let ((side,m), (_,n)) = swapOnRowMajor order (('L', m0), ('R', n0))
-   in (side,(m,n))
-
-
 mapChecked ::
    (Shape.C sha, Shape.C shb) =>
    String -> (sha -> shb) -> sha -> shb
@@ -131,7 +123,10 @@
    Full meas vert horiz height width -> width
 fullWidth = Extent.width . fullExtent
 
+squareSize :: Square sh -> sh
+squareSize = Extent.squareSize . fullExtent
 
+
 fullIndices ::
    (Extent.Measure meas, Extent.C vert, Extent.C horiz,
     Shape.Indexed a, Shape.Indexed b) =>
@@ -455,10 +450,6 @@
 hermitianFromSymmetric :: SymmetricP pack size -> HermitianP pack size
 hermitianFromSymmetric (Mosaic pack SimpleMirror upper order size) =
    Mosaic pack ConjugateMirror upper order size
-
-uploFromOrder :: Order -> Char
-uploFromOrder RowMajor = 'L'
-uploFromOrder ColumnMajor = 'U'
 
 
 
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.BLAS.Matrix.Modifier as Mod
+import qualified Numeric.Netlib.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.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.Netlib.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.BLAS.Matrix.Modifier as Mod
+import qualified Numeric.Netlib.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.BLAS.Matrix.Modifier (Transposition(NonTransposed,Transposed))
+import Numeric.Netlib.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.BLAS.Matrix.Modifier as Mod
+import qualified Numeric.Netlib.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
@@ -31,6 +31,8 @@
    Basic.takeTop, Basic.takeBottom,
    Basic.takeLeft, Basic.takeRight,
    takeRowArray, takeColumnArray,
+   takeRowSet, takeColumnSet,
+   takeRowIntSet, takeColumnIntSet,
    swapRows, swapColumns,
    reverseRows, reverseColumns,
    fromRowMajor, toRowMajor,
@@ -54,11 +56,11 @@
 import qualified Numeric.LAPACK.Matrix.Square.Basic as Square
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Matrix.Basic as Basic
-import qualified Numeric.LAPACK.Matrix.RowMajor as RowMajor
 import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Vector.Private as VectorPriv
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Shape as ExtShape
+import qualified Numeric.BLAS.Matrix.RowMajor as RowMajor
 import Numeric.LAPACK.Matrix.Layout.Private
          (Order(RowMajor, ColumnMajor), transposeFromOrder)
 import Numeric.LAPACK.Matrix.Basic
@@ -66,7 +68,7 @@
 import Numeric.LAPACK.Matrix.Private
          (Full, Tall, Wide, General, ShapeInt, shapeInt,
           mapExtent, fromFull, generalizeTall, generalizeWide)
-import Numeric.BLAS.Matrix.Modifier (Conjugation(NonConjugated,Conjugated))
+import Numeric.Netlib.Modifier (Conjugation(NonConjugated,Conjugated))
 import Numeric.LAPACK.Vector (Vector)
 import Numeric.LAPACK.Scalar (zero, one)
 import Numeric.LAPACK.Private (pointerSeq, fill, copySubMatrix, copyBlock)
@@ -96,6 +98,10 @@
 import qualified Data.NonEmpty.Mixed as NonEmptyM
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Foldable as Fold
+import qualified Data.IntSet as IntSet
+import qualified Data.Set as Set
+import Data.IntSet (IntSet)
+import Data.Set (Set)
 import Data.Foldable (forM_)
 import Data.Maybe (listToMaybe)
 import Data.Bool.HT (if')
@@ -393,19 +399,25 @@
 For scattered rows in column major order
 the function has quite ugly memory access patterns.
 -}
-takeRowArray ::
-   (Shape.Indexed height, Shape.C width, Shape.C sh, Class.Floating a) =>
-   BoxedArray.Array sh (Shape.Index height) ->
-   General height width a -> General sh width a
-takeRowArray ixs (Array (Layout.Full order extent) a) =
+takeSelectedRows ::
+   (Shape.Indexed heightA, Shape.Index heightA ~ ix,
+    Shape.C heightB, Shape.C width, Class.Floating a,
+    Extent.C vert, Extent.C horiz) =>
+   (Extent.Extent Extent.Size vert horiz heightA width ->
+    Extent.Extent Extent.Size vert horiz heightB width) ->
+   [ix] ->
+   Full Extent.Size vert horiz heightA width a ->
+   Full Extent.Size vert horiz heightB width a
+takeSelectedRows makeExtent ixs (Array (Layout.Full order extent) a) =
    let (heightA, width) = Extent.dimensions extent
-       heightB = BoxedArray.shape ixs
-       offsets = map (Shape.offset heightA) $ BoxedArray.toList ixs
+       newExtent = makeExtent extent
+       heightB = Extent.height newExtent
+       offsets = map (Shape.offset heightA) ixs
        startBlocks blocks = zip (scanl (+) 0 $ map fst blocks) blocks
        ma = Shape.size heightA
        mb = Shape.size heightB
        n = Shape.size width
-   in Array.unsafeCreate (Layout.general order heightB width) $ \bPtr ->
+   in Array.unsafeCreate (Layout.Full order newExtent) $ \bPtr ->
       withForeignPtr a $ \aPtr ->
       case order of
          RowMajor -> do
@@ -439,14 +451,66 @@
    NonEmptyM.groupBy (\i j -> i+1 == j)
 
 
+takeRowArray ::
+   (Shape.Indexed heightA, Shape.C width, Shape.C heightB, Class.Floating a) =>
+   BoxedArray.Array heightB (Shape.Index heightA) ->
+   General heightA width a -> General heightB width a
+takeRowArray ixs =
+   takeSelectedRows
+      (\extent -> Extent.general (BoxedArray.shape ixs) (Extent.width extent))
+      (BoxedArray.toList ixs)
+
+takeRowSet ::
+   (Extent.C vert,
+    Shape.Indexed height, Shape.Index height ~ ix, Ord ix,
+    Shape.C width, Class.Floating a) =>
+   Set ix ->
+   Full Extent.Size vert Extent.Big height width a ->
+   Full Extent.Size vert Extent.Big (Set ix) width a
+takeRowSet ixs =
+   takeSelectedRows
+      (Extent.reduceWideHeight ixs)
+      (Set.toAscList ixs)
+
+takeRowIntSet ::
+   (Extent.C vert,
+    Shape.Indexed height, Shape.Index height ~ Int,
+    Shape.C width, Class.Floating a) =>
+   IntSet ->
+   Full Extent.Size vert Extent.Big height width a ->
+   Full Extent.Size vert Extent.Big IntSet width a
+takeRowIntSet ixs =
+   takeSelectedRows
+      (Extent.reduceWideHeight ixs)
+      (IntSet.toAscList ixs)
+
+
 takeColumnArray ::
    (Shape.C height, Shape.Indexed width, Shape.C sh, Class.Floating a) =>
    BoxedArray.Array sh (Shape.Index width) ->
    General height width a -> General height sh a
 takeColumnArray ixs = transpose . takeRowArray ixs . transpose
 
+takeColumnSet ::
+   (Extent.C horiz,
+    Shape.Indexed width, Shape.Index width ~ ix, Ord ix,
+    Shape.C height, Class.Floating a) =>
+   Set ix ->
+   Full Extent.Size Extent.Big horiz height width a ->
+   Full Extent.Size Extent.Big horiz height (Set ix) a
+takeColumnSet ixs = transpose . takeRowSet ixs . transpose
 
+takeColumnIntSet ::
+   (Extent.C horiz,
+    Shape.Indexed width, Shape.Index width ~ Int,
+    Shape.C height, Class.Floating a) =>
+   IntSet ->
+   Full Extent.Size Extent.Big horiz height width a ->
+   Full Extent.Size Extent.Big horiz height IntSet a
+takeColumnIntSet ixs = transpose . takeRowIntSet ixs . transpose
 
+
+
 fromRowMajor ::
    (Shape.C height, Shape.C width, Storable a) =>
    Array (height,width) a -> General height width a
@@ -600,9 +664,12 @@
        shc =
          Layout.Full orderB $
          Extent.kronecker (Layout.fullExtent $ Array.shape a) extentB
+       at = Basic.transposableFromFull a
    in either
-         (Array.reshape shc . RowMajor.kronecker a)
-         (Array.reshape shc . RowMajor.kronecker (transpose a)) $
+         (Array.reshape shc . RowMajor.kroneckerLeftTransposable at)
+         (Array.reshape shc .
+            RowMajor.kroneckerLeftTransposable
+               (RowMajor.transposeTransposable at)) $
       Matrix.revealOrder b
 
 
diff --git a/src/Numeric/LAPACK/Matrix/RowMajor.hs b/src/Numeric/LAPACK/Matrix/RowMajor.hs
deleted file mode 100644
--- a/src/Numeric/LAPACK/Matrix/RowMajor.hs
+++ /dev/null
@@ -1,80 +0,0 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE TypeOperators #-}
-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 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 (advancePtr)
-import Foreign.ForeignPtr (withForeignPtr)
-
-import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
-import Control.Monad.IO.Class (liftIO)
-import Control.Applicative (liftA2)
-
-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))
-
-import Data.Foldable (forM_)
-
-
-kronecker ::
-   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
-    Shape.C heightA, Shape.C widthA, Shape.C heightB, Shape.C widthB,
-    Class.Floating a) =>
-   Full meas vert horiz heightA widthA a ->
-   Matrix heightB widthB a ->
-   Matrix (heightA,heightB) (widthA,widthB) a
-kronecker
-      (Array (Layout.Full orderA extentA) a) (Array (heightB,widthB) b) =
-   let (heightA,widthA) = Extent.dimensions extentA
-   in Array.unsafeCreate ((heightA,heightB), (widthA,widthB)) $ \cPtr ->
-      evalContT $ do
-   let (ma,na) = (Shape.size heightA, Shape.size widthA)
-   let (mb,nb) = (Shape.size heightB, Shape.size widthB)
-   let (lda,istep) =
-         case orderA of
-            RowMajor -> (1,na)
-            ColumnMajor -> (ma,1)
-   transaPtr <- Call.char 'N'
-   transbPtr <- Call.char 'T'
-   mPtr <- Call.cint na
-   nPtr <- Call.cint nb
-   kPtr <- Call.cint 1
-   alphaPtr <- Call.number one
-   aPtr <- ContT $ withForeignPtr a
-   ldaPtr <- Call.leadingDim lda
-   bPtr <- ContT $ withForeignPtr b
-   ldbPtr <- Call.leadingDim 1
-   betaPtr <- Call.number zero
-   ldcPtr <- Call.leadingDim nb
-   liftIO $
-      forM_ (liftA2 (,) (take ma [0..]) (take mb [0..])) $ \(i,j) -> do
-         let aiPtr = advancePtr aPtr (istep*i)
-         let bjPtr = advancePtr bPtr (nb*j)
-         let cijPtr = advancePtr cPtr (na*nb*(j+mb*i))
-         BlasGen.gemm
-            transbPtr transaPtr nPtr mPtr kPtr alphaPtr
-            bjPtr ldbPtr aiPtr ldaPtr betaPtr cijPtr ldcPtr
diff --git a/src/Numeric/LAPACK/Matrix/Square/Basic.hs b/src/Numeric/LAPACK/Matrix/Square/Basic.hs
--- a/src/Numeric/LAPACK/Matrix/Square/Basic.hs
+++ b/src/Numeric/LAPACK/Matrix/Square/Basic.hs
@@ -41,8 +41,8 @@
 import qualified Numeric.LAPACK.Matrix.Private as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Private as Private
-import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor, ColumnMajor), swapOnRowMajor)
+import qualified Numeric.BLAS.Matrix.RowMajor as RowMajor
+import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor, ColumnMajor))
 import Numeric.LAPACK.Matrix.Private
          (General, argGeneral, LiberalSquare, Square, argSquare,
           Full, mapExtent, ShapeInt, shapeInt)
@@ -71,6 +71,7 @@
 import Control.Monad.IO.Class (liftIO)
 
 import Data.Function.HT (powerAssociative)
+import Data.Tuple.HT (double)
 
 
 mapSize :: (sh0 -> sh1) -> Square sh0 a -> Square sh1 a
@@ -145,7 +146,7 @@
 identityFromHeight =
    argGeneral $ \order height _ _ -> identityOrder order height
 
-identityOrder, _identityOrder ::
+identityOrder ::
    (Shape.C sh, Class.Floating a) => Order -> sh -> Square sh a
 identityOrder order sh =
    Array.unsafeCreate (Layout.square order sh) $ \aPtr ->
@@ -156,22 +157,6 @@
       betaPtr <- Call.number one
       liftIO $ LapackGen.laset uploPtr nPtr nPtr alphaPtr betaPtr aPtr nPtr
 
-_identityOrder order sh =
-   Array.unsafeCreateWithSize (Layout.square order sh) $ \blockSize yPtr ->
-   evalContT $ do
-      nPtr <- Call.alloca
-      xPtr <- Call.number zero
-      incxPtr <- Call.cint 0
-      incyPtr <- Call.cint 1
-      liftIO $ do
-         pokeCInt nPtr blockSize
-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
-         let n = fromIntegral $ Shape.size sh
-         poke nPtr n
-         poke xPtr one
-         poke incyPtr (n+1)
-         BlasGen.copy nPtr xPtr incxPtr yPtr incyPtr
-
 diagonal :: (Shape.C sh, Class.Floating a) => Vector sh a -> Square sh a
 diagonal (Array sh x) =
    Array.unsafeCreateWithSize (Layout.square ColumnMajor sh) $
@@ -252,10 +237,11 @@
 multiplyCommutativeUnchecked ::
    (Shape.C sh, Class.Floating a) =>
    Square sh a -> Square sh a -> Square sh a
-multiplyCommutativeUnchecked
-   (Array shape@(Layout.Full order extent) a)
-   (Array _ b) =
-      Array.unsafeCreate shape $ \cPtr ->
-   let n = Shape.size $ Extent.height extent
-       (at,bt) = swapOnRowMajor order (a,b)
-   in  Private.multiplyMatrix ColumnMajor ColumnMajor n n n at bt cPtr
+multiplyCommutativeUnchecked a b =
+   Array.reshape (Array.shape b) $
+   RowMajor.multiply
+      (array2FromSquare $ Basic.uncheck a)
+      (array2FromSquare $ Basic.uncheck b)
+
+array2FromSquare :: Square sh a -> RowMajor.Square sh a
+array2FromSquare = Array.mapShape (double . Layout.squareSize)
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.BLAS.Matrix.Modifier (Conjugation(NonConjugated))
+import Numeric.Netlib.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
@@ -20,9 +20,10 @@
           noLabel, label, applyFuncPair, triArg,
           Labelled(Labelled))
 import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor,ColumnMajor), uploFromOrder, uploOrder)
-import Numeric.BLAS.Matrix.Modifier
-         (Conjugation(NonConjugated, Conjugated), conjugatedOnRowMajor,
+         (Order(RowMajor,ColumnMajor), uploFromOrder, uploOrder,
+          conjugatedOnRowMajor)
+import Numeric.Netlib.Modifier
+         (Conjugation(NonConjugated, Conjugated),
           Transposition(NonTransposed, Transposed), transposeOrder)
 import Numeric.LAPACK.Matrix.Private (Full, General, Square)
 import Numeric.LAPACK.Linear.Private (solver, withDeterminantInfo, diagonalMsg)
diff --git a/src/Numeric/LAPACK/Matrix/Type/Private.hs b/src/Numeric/LAPACK/Matrix/Type/Private.hs
--- a/src/Numeric/LAPACK/Matrix/Type/Private.hs
+++ b/src/Numeric/LAPACK/Matrix/Type/Private.hs
@@ -191,7 +191,7 @@
 
 instance Format Permutation where
    type FormatExtra Permutation extra = ()
-   format _fmt (Permutation perm) = Perm.format perm
+   format cfg (Permutation perm) = Perm.format cfg perm
 
 
 {- |
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.BLAS.Matrix.Modifier as Mod
+import qualified Numeric.Netlib.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.BLAS.Matrix.Modifier (Transposition, Conjugation)
+import Numeric.Netlib.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.BLAS.Matrix.Modifier
+import Numeric.Netlib.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.BLAS.Matrix.Modifier as Mod
+import qualified Numeric.Netlib.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
@@ -5,9 +5,10 @@
 import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
 import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
 import qualified Numeric.LAPACK.Output as Output
+import Numeric.LAPACK.Matrix.Plain.Format (Config(configEmpty))
 import Numeric.LAPACK.Output (Output, formatAligned)
 import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor, ColumnMajor))
-import Numeric.BLAS.Matrix.Modifier
+import Numeric.Netlib.Modifier
          (Transposition(NonTransposed,Transposed),
           Inversion(NonInverted,Inverted))
 import Numeric.LAPACK.Matrix.Private (Full, Square)
@@ -68,10 +69,10 @@
 newtype Permutation sh = Permutation (Vector (Shape sh) (Element sh))
    deriving (Eq, Show)
 
-format :: (Shape.C sh, Output out) => Permutation sh -> out
-format (Permutation perm) =
+format :: (Shape.C sh, Output out) => Config -> Permutation sh -> out
+format cfg (Permutation perm) =
    let n = Shape.size $ Array.shape perm
-       s0 = Output.text "."
+       s0 = Output.text $ configEmpty cfg
        s1 = Output.text "1"
    in formatAligned $
       map (\k -> map Identity $ replicate (k-1) s0 ++ s1 : replicate (n-k) s0) $
diff --git a/src/Numeric/LAPACK/Private.hs b/src/Numeric/LAPACK/Private.hs
--- a/src/Numeric/LAPACK/Private.hs
+++ b/src/Numeric/LAPACK/Private.hs
@@ -1,8 +1,7 @@
 {-# LANGUAGE TypeFamilies #-}
 module Numeric.LAPACK.Private where
 
-import Numeric.LAPACK.Matrix.Layout.Private
-         (Order(RowMajor, ColumnMajor), transposeFromOrder)
+import Numeric.LAPACK.Matrix.Layout.Private (Order(RowMajor, ColumnMajor))
 import Numeric.LAPACK.Wrapper (Flip(Flip, getFlip))
 
 import qualified Numeric.LAPACK.FFI.Generic as LapackGen
@@ -411,35 +410,6 @@
    when (n == 0 && isZero beta) $
       Marshal.with 0 $ \incbPtr ->
       BlasGen.copy mtPtr betaPtr incbPtr yPtr incyPtr
-
-
-multiplyMatrix ::
-   (Class.Floating a) =>
-   Order -> Order -> Int -> Int -> Int ->
-   ForeignPtr a -> ForeignPtr a -> Ptr a -> IO ()
-multiplyMatrix orderA orderB m k n a b cPtr = evalContT $ do
-   let lda = case orderA of RowMajor -> k; ColumnMajor -> m
-   let ldb = case orderB of RowMajor -> n; ColumnMajor -> k
-   let ldc = m
-   if k==0
-      then liftIO $ fill zero (m*n) cPtr
-      else do
-      transaPtr <- Call.char $ transposeFromOrder orderA
-      transbPtr <- Call.char $ transposeFromOrder orderB
-      mPtr <- Call.cint m
-      nPtr <- Call.cint n
-      kPtr <- Call.cint k
-      alphaPtr <- Call.number one
-      aPtr <- ContT $ withForeignPtr a
-      ldaPtr <- Call.leadingDim lda
-      bPtr <- ContT $ withForeignPtr b
-      ldbPtr <- Call.leadingDim ldb
-      betaPtr <- Call.number zero
-      ldcPtr <- Call.leadingDim ldc
-      liftIO $
-         BlasGen.gemm
-            transaPtr transbPtr mPtr nPtr kPtr alphaPtr aPtr ldaPtr
-            bPtr ldbPtr betaPtr cPtr ldcPtr
 
 
 
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.BLAS.Matrix.Modifier
+import Numeric.Netlib.Modifier
          (Transposition, transposeOrder,
           Conjugation(NonConjugated, Conjugated))
 import Numeric.LAPACK.Matrix.Private (Full)
diff --git a/test-module.list b/test-module.list
--- a/test-module.list
+++ b/test-module.list
@@ -1,3 +1,3 @@
-Numeric.LAPACK.Permutation.Private
-Numeric.LAPACK.Example.DividedDifference
 Numeric.LAPACK.Example.EconomicAllocation
+Numeric.LAPACK.Example.DividedDifference
+Numeric.LAPACK.Permutation.Private
diff --git a/test/DocTest/Main.hs b/test/DocTest/Main.hs
--- a/test/DocTest/Main.hs
+++ b/test/DocTest/Main.hs
@@ -1,14 +1,14 @@
 -- Do not edit! Automatically created with doctest-extract.
 module DocTest.Main where
 
-import qualified DocTest.Numeric.LAPACK.Permutation.Private
-import qualified DocTest.Numeric.LAPACK.Example.DividedDifference
 import qualified DocTest.Numeric.LAPACK.Example.EconomicAllocation
+import qualified DocTest.Numeric.LAPACK.Example.DividedDifference
+import qualified DocTest.Numeric.LAPACK.Permutation.Private
 
 import qualified Test.DocTest.Driver as DocTest
 
 main :: DocTest.T ()
 main = do
-    DocTest.Numeric.LAPACK.Permutation.Private.test
-    DocTest.Numeric.LAPACK.Example.DividedDifference.test
     DocTest.Numeric.LAPACK.Example.EconomicAllocation.test
+    DocTest.Numeric.LAPACK.Example.DividedDifference.test
+    DocTest.Numeric.LAPACK.Permutation.Private.test
diff --git a/test/DocTest/Numeric/LAPACK/Permutation/Private.hs b/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
--- a/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
+++ b/test/DocTest/Numeric/LAPACK/Permutation/Private.hs
@@ -1,11 +1,11 @@
 -- Do not edit! Automatically created with doctest-extract from src/Numeric/LAPACK/Permutation/Private.hs
-{-# LINE 48 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 49 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
 
 module DocTest.Numeric.LAPACK.Permutation.Private where
 
 import qualified Test.DocTest.Driver as DocTest
 
-{-# LINE 49 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 50 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
 import     qualified Test.QuickCheck as QC
 import     Test.Permutation (genPerm, genPivots)
 
@@ -25,21 +25,21 @@
 
 test :: DocTest.T ()
 test = do
- DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:110: "
-{-# LINE 110 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:111: "
+{-# LINE 111 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
  DocTest.property(
-{-# LINE 110 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 111 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
       QC.forAll QC.arbitraryBoundedEnum $ \inv -> QC.forAll (QC.arbitrary >>= genPivots) $ \xs -> xs == Perm.toPivots inv (Perm.fromPivots inv xs)
   )
- DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:193: "
-{-# LINE 193 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:194: "
+{-# LINE 194 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
  DocTest.property(
-{-# LINE 193 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 194 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
       QC.forAll genPerm2 $ \(p0,p1) -> determinant (multiply p0 p1) == determinant p0 <> determinant p1
   )
- DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:222: "
-{-# LINE 222 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+ DocTest.printPrefix "Numeric.LAPACK.Permutation.Private:223: "
+{-# LINE 223 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
  DocTest.property(
-{-# LINE 222 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
+{-# LINE 223 "src/Numeric/LAPACK/Permutation/Private.hs" #-}
       QC.forAll genPerm2 $ \(p0,p1) -> transpose (multiply p0 p1) == multiply (transpose p1) (transpose p0)
   )