diff --git a/example/Fit.hs b/example/Fit.hs
--- a/example/Fit.hs
+++ b/example/Fit.hs
@@ -1,31 +1,50 @@
+{-# LANGUAGE Rank2Types #-}
 module Main where
 
+import qualified Plot2DExtra
+
 import qualified Numeric.Interpolation.NodeList as Nodes
 import qualified Numeric.Interpolation.Piecewise as Piecewise
 import qualified Numeric.Interpolation.Type as Type
 
-import qualified Data.Packed.ST as PackST
-import qualified Data.Packed.Matrix as Matrix
-import qualified Data.Packed.Vector as Vector
-import Data.Packed.Matrix (Matrix)
-import Data.Packed.Vector (Vector)
+import qualified Numeric.LAPACK.Singular as Singular
+import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.BandedHermitianPositiveDefinite
+                                                       as BandedSPD
+import qualified Numeric.LAPACK.Matrix.BandedHermitian as BandedHermitian
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
 
-import qualified Numeric.LinearAlgebra.Banded as Banded
-import qualified Numeric.Container as Container
-import Numeric.Container ((<\>))
+import qualified Type.Data.Num.Unary as Unary
+import Type.Base.Proxy (Proxy(Proxy))
 
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+
 import qualified Graphics.Gnuplot.Advanced as GP
+import qualified Graphics.Gnuplot.Frame as Frame
+import qualified Graphics.Gnuplot.Frame.OptionSet as Opts
 import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
 import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
+import qualified Graphics.Gnuplot.LineSpecification as LineSpec
 
 import System.Random (randomRs, mkStdGen)
 import Control.Monad.HT (void)
-import Control.Monad (when, zipWithM_, forM_)
+import Control.Applicative ((<$>))
 
 import qualified Data.Foldable as Fold
+import Data.Tuple.HT (mapFst)
 import Data.Monoid ((<>))
 
 
+type BandedHermitianMatrix k =
+       BandedHermitian.BandedHermitian k Matrix.ZeroInt
+type Matrix = Matrix.General Matrix.ZeroInt Matrix.ZeroInt
+type Vector = Vector.Vector Matrix.ZeroInt
+type ShortVector k = Vector.Vector (ShapeStatic.ZeroBased k)
+
+
 noisy :: [(Double, Double)]
 noisy =
    take 100 $
@@ -37,19 +56,24 @@
 basisMatrixFull ::
    Type.T Double Double ny -> [Double] -> [Double] -> Matrix Double
 basisMatrixFull typ xs txs0 =
-   let txs = Vector.fromList txs0
-   in  Matrix.fromColumns $
-       map (flip Container.cmap txs . Piecewise.interpolateConstantExt typ) $
+   let txs = Vector.autoFromList txs0
+   in  Matrix.fromColumns (Array.shape txs) $
+       map (flip Array.map txs . Piecewise.interpolateConstantExt typ) $
        Type.basisFunctions typ xs
 
+zipRowsWith :: (a -> b -> c) -> [a] -> [[b]] -> [c]
+zipRowsWith f as bs = concat $ zipWith (map . f) as bs
+
+-- ToDo: generalize shapes
 basisMatrixSparse ::
    Type.T Double Double ny -> [Double] -> [Double] -> Matrix Double
-basisMatrixSparse typ xs txs = PackST.runSTMatrix $ do
-   mat <- PackST.newMatrix 0 (length txs) (length $ Type.basisFunctions typ xs)
-   zipWithM_
-      (\k -> mapM_ (uncurry (PackST.writeMatrix mat k))) [0..] $
-      map (Type.sampleBasisFunctions typ xs) txs
-   return mat
+basisMatrixSparse typ xs txs =
+   Array.fromAssociations
+      (MatrixShape.general MatrixShape.RowMajor
+         (Shape.ZeroBased $ length txs)
+         (Shape.ZeroBased $ length $ Type.basisFunctions typ xs)) 0 $
+   zipRowsWith (\k (j,x) -> ((k,j),x)) [0..] $
+   map (Type.sampleBasisFunctions typ xs) txs
 
 fit ::
    Type.T Double Double ny ->
@@ -58,39 +82,56 @@
    let (txs, tys) = unzip target
        matrix = basisMatrixSparse typ xs txs
    in  Type.coefficientsToInterpolator typ xs $
-       Vector.toList $ matrix <\> Vector.fromList tys
+       Vector.toList $
+       Matrix.unliftColumn MatrixShape.ColumnMajor
+         (snd . Singular.leastSquaresMinimumNormRCond 1e-5 matrix) $
+       Vector.autoFromList tys
 
 matrixDiff ::
-   Type.T Double Double ny ->
-   [Double] -> [(Double, Double)] -> Double
+   Type.T Double Double ny -> [Double] -> [(Double, Double)] -> Double
 matrixDiff typ xs target =
    let txs = map fst target
-   in  Container.maxElement $ Container.cmap abs $
-       Container.sub
+   in  snd $ Vector.argAbs1Maximum $
+       Matrix.sub
           (basisMatrixFull typ xs txs)
           (basisMatrixSparse typ xs txs)
 
 
 mulSparseMatrixVector ::
    Int -> [[(Int, Double)]] -> [Double] -> Vector Double
-mulSparseMatrixVector size samples tys = PackST.runSTVector $ do
-   vec <- PackST.newVector 0 size
-   forM_ (zip samples tys) $ \(row,ty) ->
-      forM_ row $ \(k,y) ->
-         PackST.modifyVector vec k (+y*ty)
-   return vec
+mulSparseMatrixVector size samples tys =
+   Array.accumulate (+)
+      (Vector.constant (Matrix.zeroInt size) 0)
+      (zipRowsWith (\ty (k,y) -> (k, y*ty)) tys samples)
 
+{- ToDo:
+use index type (ZeroInt, Enumeration Order)
+with Order = Absolute | Derivative
+Problem: not all interpolation types use derivatives
+-}
+shortVector ::
+   (Unary.Natural k) =>
+   Proxy k -> [(Int, Double)] -> (Int, ShortVector k Double)
+shortVector width xs =
+   let i0 = minimum $ map fst xs
+   in (i0,
+--       Array.fromAssociations (ShapeStatic.ZeroBased Proxy) 0 $
+       Array.reshape (ShapeStatic.ZeroBased Proxy) $
+       Array.fromAssociations (Matrix.zeroInt $ Unary.integralFromProxy width) 0 $
+       map (mapFst (subtract i0)) xs)
+
 bandedGramian ::
-   Int -> Int -> [[(Int, Double)]] -> Banded.SymmetricMatrix Double
+   (Unary.Natural k) =>
+   Int -> Proxy (Unary.Succ k) ->
+   [[(Int, Double)]] -> BandedHermitianMatrix k Double
 bandedGramian size width samples =
-      Banded.SymmetricMatrix $ PackST.runSTMatrix $ do
-   mat <- PackST.newMatrix 0 size width
-   forM_ samples $ \row ->
-      forM_ row $ \(k,yk) ->
-      forM_ row $ \(j,yj) ->
-         when (k<=j) $ PackST.modifyMatrix mat k (j-k) (+yk*yj)
-   return mat
+   BandedHermitian.sumRank1 MatrixShape.ColumnMajor (Matrix.zeroInt size) $
+   map ((,) 1 . shortVector width) samples
 
+reifyPositive ::
+   Integer -> (forall s. Unary.Natural s => Proxy (Unary.Succ s) -> w) -> w
+reifyPositive n f = Unary.reifyNatural (n-1) (f . Unary.succ)
+
 fitBanded ::
    Type.T Double Double ny ->
    [Double] -> [(Double, Double)] -> Nodes.T Double ny
@@ -98,11 +139,12 @@
    let size = length $ Type.basisFunctions typ xs
        (txs, tys) = unzip target
        samples = map (Type.sampleBasisFunctions typ xs) txs
-       matrix =
-          Banded.choleskyDecompose $
-          bandedGramian size (Type.basisOverlap typ) samples
-   in  Type.coefficientsToInterpolator typ xs $ Vector.toList $
-       Banded.choleskySolve matrix $ mulSparseMatrixVector size samples tys
+   in reifyPositive (toInteger $ Type.basisOverlap typ) 
+         (\width ->
+            Type.coefficientsToInterpolator typ xs $ Vector.toList $
+            Matrix.unliftColumn MatrixShape.ColumnMajor
+               (BandedSPD.solve (bandedGramian size width samples)) $
+            mulSparseMatrixVector size samples tys)
 
 bandedDiff ::
    (ny -> ny -> Double) ->
@@ -122,44 +164,62 @@
 
 
 plotBasisFunctions ::
-   Type.T Double Double ny -> [Double] -> Plot2D.T Double Double
-plotBasisFunctions nodeType xs =
+   (String, Type.T Double Double ny) -> [Double] ->
+   Frame.T (Graph2D.T Double Double)
+plotBasisFunctions (title, nodeType) xs =
    let abscissa = Plot2D.linearScale 1000 (minimum xs, maximum xs)
-   in  Plot2D.functions Graph2D.lines abscissa $
+   in  Frame.cons (Opts.title title $ Opts.key False Opts.deflt) $
+       Plot2D.functions Graph2D.lines abscissa $
        map (Piecewise.interpolateConstantExt nodeType) $
        Type.basisFunctions nodeType xs
 
 
+typeLinear :: (String, Type.T Double Double Double)
+typeLinear        = ("linear",        Type.linear)
+typeHermite1, typeCubicLinear, typeCubicParabola ::
+   (String, Type.T Double Double (Double, Double))
+typeHermite1      = ("hermite1",      Type.hermite1)
+typeCubicLinear   = ("cubicLinear",   Type.cubicLinear)
+typeCubicParabola = ("cubicParabola", Type.cubicParabola)
+
+
 main :: IO ()
 main = do
    let xs = [0, 1, 3, 4, 6, 7]
        exs = (-1) : xs ++ [8]
-   void $ GP.plotDefault $ plotBasisFunctions Type.linear xs
-   void $ GP.plotDefault $ plotBasisFunctions Type.hermite1 xs
-   void $ GP.plotDefault $ plotBasisFunctions Type.cubicLinear exs
-   void $ GP.plotDefault $ plotBasisFunctions Type.cubicParabola exs
+   void $ GP.plotDefault $ plotBasisFunctions typeLinear xs
+   void $ GP.plotDefault $ plotBasisFunctions typeHermite1 xs
+   void $ GP.plotDefault $ plotBasisFunctions typeCubicLinear exs
+   void $ GP.plotDefault $ plotBasisFunctions typeCubicParabola exs
    let linearNodes = fit Type.linear xs noisy
        hermite1Nodes = fit Type.hermite1 xs noisy
        cubicLinearNodes = fit Type.cubicLinear exs noisy
        cubicParabolaNodes = fit Type.cubicParabola exs noisy
    void $ GP.plotDefault $
-      Plot2D.list Graph2D.points noisy
+      (Graph2D.lineSpec (LineSpec.title "target" LineSpec.deflt)
+         <$> Plot2D.list Graph2D.points noisy)
       <>
-      (Plot2D.functions Graph2D.lines (Plot2D.linearScale 1000 (-2,10)) $
-       Piecewise.interpolateConstantExt Type.linear linearNodes :
-       Piecewise.interpolateConstantExt Type.hermite1 hermite1Nodes :
-       Piecewise.interpolateConstantExt Type.cubicLinear cubicLinearNodes :
-       Piecewise.interpolateConstantExt Type.cubicParabola cubicParabolaNodes :
-       [])
+      (let interpolate (name,typ) nodes =
+             (name, Piecewise.interpolateConstantExt typ nodes)
+       in Plot2DExtra.functions (Plot2D.linearScale 1000 (-2,10)) $
+            interpolate typeLinear linearNodes :
+            interpolate typeHermite1 hermite1Nodes :
+            interpolate typeCubicLinear cubicLinearNodes :
+            interpolate typeCubicParabola cubicParabolaNodes :
+            [])
 
    putStrLn "differences between matrices should be almost zero:"
-   putStrLn $ "linear: " ++ show (matrixDiff Type.linear xs noisy)
-   putStrLn $ "hermite1: " ++ show (matrixDiff Type.hermite1 xs noisy)
-   putStrLn $ "cubicLinear: " ++ show (matrixDiff Type.cubicLinear exs noisy)
-   putStrLn $ "cubicParabola: " ++ show (matrixDiff Type.cubicParabola exs noisy)
+   let printMatrixDiff (name,typ) ps =
+         putStrLn $ name ++ ": " ++ show (matrixDiff typ ps noisy)
+   printMatrixDiff typeLinear xs
+   printMatrixDiff typeHermite1 xs
+   printMatrixDiff typeCubicLinear exs
+   printMatrixDiff typeCubicParabola exs
 
    putStrLn "differences between samples should be almost zero:"
-   putStrLn $ "linear: " ++ show (bandedDiff absDiffSingle Type.linear xs noisy)
-   putStrLn $ "hermite1: " ++ show (bandedDiff absDiffPair Type.hermite1 xs noisy)
-   putStrLn $ "cubicLinear: " ++ show (bandedDiff absDiffPair Type.cubicLinear exs noisy)
-   putStrLn $ "cubicParabola: " ++ show (bandedDiff absDiffPair Type.cubicParabola exs noisy)
+   let printBandedDiff diff (name,typ) ps =
+         putStrLn $ name ++ ": " ++ show (bandedDiff diff typ ps noisy)
+   printBandedDiff absDiffSingle typeLinear xs
+   printBandedDiff absDiffPair typeHermite1 xs
+   printBandedDiff absDiffPair typeCubicLinear exs
+   printBandedDiff absDiffPair typeCubicParabola exs
diff --git a/example/Plot.hs b/example/Plot.hs
--- a/example/Plot.hs
+++ b/example/Plot.hs
@@ -1,5 +1,7 @@
 module Main where
 
+import qualified Plot2DExtra
+
 import qualified Numeric.Interpolation.NodeList as Nodes
 import qualified Numeric.Interpolation.Piecewise as Piecewise
 import qualified Numeric.Interpolation.Basis as Basis
@@ -9,7 +11,6 @@
 import qualified Graphics.Gnuplot.Advanced as GP
 
 import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
-import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
 
 import Control.Monad.HT (void)
 
@@ -24,10 +25,11 @@
        cubicLinearNodes = Basis.coefficientsToCubicLinear xs $ map sin xs
        cubicParabolaNodes = Basis.coefficientsToCubicParabola xs $ map sin xs
    void $ GP.plotDefault $
-      Plot2D.functions Graph2D.lines
-         (Plot2D.linearScale 1000 (-2,15))
-         [Piecewise.interpolateConstantExt Type.linear linearNodes,
-          Piecewise.interpolateConstantExt Type.hermite1 hermite1Nodes,
-          Piecewise.interpolateConstantExt Type.hermite1 cubicLinearNodes,
-          Piecewise.interpolateConstantExt Type.hermite1 cubicParabolaNodes,
-          sin]
+      Plot2DExtra.functions (Plot2D.linearScale 1000 (-2,15)) $
+         let interpolate typ nodes = Piecewise.interpolateConstantExt typ nodes
+         in ("linear",        interpolate Type.linear linearNodes) :
+            ("hermite1",      interpolate Type.hermite1 hermite1Nodes) :
+            ("cubicLinear",   interpolate Type.hermite1 cubicLinearNodes) :
+            ("cubicParabola", interpolate Type.hermite1 cubicParabolaNodes) :
+            ("sin", sin) :
+            []
diff --git a/example/Plot2DExtra.hs b/example/Plot2DExtra.hs
new file mode 100644
--- /dev/null
+++ b/example/Plot2DExtra.hs
@@ -0,0 +1,17 @@
+module Plot2DExtra where
+
+import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
+import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
+import qualified Graphics.Gnuplot.LineSpecification as LineSpec
+import qualified Graphics.Gnuplot.Value.Tuple as Tuple
+import qualified Graphics.Gnuplot.Value.Atom as Atom
+
+import Data.Tuple.HT (mapFst)
+
+
+functions ::
+   (Atom.C x, Atom.C y, Tuple.C x, Tuple.C y) =>
+   [x] -> [(String, x -> y)] -> Plot2D.T x y
+functions scale =
+   Plot2D.functionsWithLineSpec Graph2D.lines scale .
+   map (mapFst (flip LineSpec.title LineSpec.deflt))
diff --git a/interpolation.cabal b/interpolation.cabal
--- a/interpolation.cabal
+++ b/interpolation.cabal
@@ -1,5 +1,5 @@
 Name:             interpolation
-Version:          0.1.0.3
+Version:          0.1.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann
@@ -26,7 +26,7 @@
   * @example/Fit.hs@:
     Demonstrates how to use the basis functions
     for fitting an interpolation function to a given function
-    using a linear least squares solver like @<\>@ from @hmatrix@.
+    using a linear least squares solver like @<\>@ from @lapack@.
     We use a distorted sinus as target.
   .
   The package needs only Haskell 98.
@@ -42,12 +42,12 @@
   description: Build example executables
   default:     False
 
-Flag hmatrix
-  description: Build examples that depend on hmatrix
+Flag lapack
+  description: Build examples that depend on lapack
   default:     True
 
 Source-Repository this
-  Tag:         0.1.0.3
+  Tag:         0.1.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/interpolation/
 
@@ -78,12 +78,13 @@
 
 Executable interpolation-plot
   Main-Is:             Plot.hs
+  Other-Modules:       Plot2DExtra
   Hs-Source-Dirs:      example
   GHC-Options:         -Wall
   If flag(buildExamples)
     Build-Depends:
       interpolation,
-      gnuplot >=0.5.2 && <0.6,
+      gnuplot >=0.5.6 && <0.6,
       utility-ht >=0.0.9 && <0.1,
       base >=4.5 && <5
   Else
@@ -91,15 +92,17 @@
 
 Executable interpolation-fit
   Main-Is:             Fit.hs
+  Other-Modules:       Plot2DExtra
   Hs-Source-Dirs:      example
   GHC-Options:         -Wall
-  If flag(buildExamples) && flag(hmatrix)
+  If flag(buildExamples) && flag(lapack)
     Build-Depends:
       interpolation,
-      hmatrix-banded >=0.0 && <0.1,
-      hmatrix >=0.15 && <0.16,
+      lapack >=0.2.2 && <0.3,
+      tfp >=1.0 && <1.1,
+      comfort-array >=0.3 && <0.4,
       random >=1.0 && <1.2,
-      gnuplot >=0.5.2 && <0.6,
+      gnuplot >=0.5.6 && <0.6,
       utility-ht >=0.0.9 && <0.1,
       base >=4.5 && <5
   Else