diff --git a/linear-circuit.cabal b/linear-circuit.cabal
--- a/linear-circuit.cabal
+++ b/linear-circuit.cabal
@@ -1,5 +1,5 @@
 Name:                linear-circuit
-Version:             0.0.0.1
+Version:             0.1
 Synopsis:            Compute resistance of linear electrical circuits
 Description:
   Compute resistance of linear electrical circuits.
@@ -15,7 +15,7 @@
 Cabal-Version:       >=1.10
 
 Source-Repository this
-  Tag:         0.0.0.1
+  Tag:         0.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/linear-circuit
 
@@ -27,9 +27,12 @@
   Exposed-Modules:
     Math.LinearCircuit
   Build-Depends:
+    lapack >=0.2.2 && <0.3,
+    netlib-ffi >=0.1.1 && <0.2,
+    comfort-array >=0.3.1 && <0.4,
     comfort-graph >=0.0 && <0.1,
-    hmatrix >=0.16 && <0.17,
     containers >=0.4 && <0.7,
+    transformers >=0.4 && <0.6,
     utility-ht >=0.0.11 && <0.1,
     base >=4.5 && <5
   Hs-Source-Dirs:      src
diff --git a/src/Math/LinearCircuit.hs b/src/Math/LinearCircuit.hs
--- a/src/Math/LinearCircuit.hs
+++ b/src/Math/LinearCircuit.hs
@@ -1,96 +1,76 @@
+{-# LANGUAGE TypeOperators #-}
 module Math.LinearCircuit (resistance) where
 
 import qualified Data.Graph.Comfort as Graph
 import Data.Graph.Comfort (Graph)
 
-import qualified Numeric.Container as NC
-import qualified Numeric.LinearAlgebra.HMatrix as HMatrix
-import qualified Data.Packed.Matrix as Matrix
-import qualified Data.Packed.Vector as Vector
-import Numeric.LinearAlgebra.HMatrix (Field, (<\>))
-import Data.Packed.Matrix (Matrix)
-import Data.Packed.Vector (Vector)
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Map as Map
-import qualified Data.List as List
-import Data.Monoid (mconcat)
+import qualified Data.Array.Comfort.Storable as Array
+import Data.Array.Comfort.Storable ((!))
+import Data.Array.Comfort.Shape ((:+:))
 
-import Control.Functor.HT (outerProduct)
+import Control.Monad.Trans.Identity (IdentityT(IdentityT))
+
+import qualified Data.Map as Map
+import qualified Data.Set as Set
+import Data.Set (Set)
 import Data.Bool.HT (if')
 
 
+type Wrap = IdentityT
+
 voltageMatrix ::
-   (Graph.Edge edge, Ord node, Field a) =>
-   Graph edge node a nodeLabel -> Matrix a
+   (Graph.Edge edge, Ord node, Class.Floating a) =>
+   Graph edge node a nodeLabel ->
+   Matrix.General (Set (Wrap edge node)) (Set node) a
 voltageMatrix gr =
-   Matrix.fromLists $
-   outerProduct
-      (\e n ->
-         if' (Graph.from e == n) 1 $
-         if' (Graph.to   e == n) (-1) $
-         0)
-      (Graph.edges gr)
-      (Graph.nodes gr)
-
-{- |
-It is almost currentMatrix = trans voltageMatrix,
-except that a row is deleted in currentMatrix.
--}
-currentMatrix ::
-   (Graph.Edge edge, Ord node, Field a) =>
-   Graph edge node a nodeLabel -> node -> node -> Matrix a
-currentMatrix gr _src dst =
-   Matrix.fromLists $
-   outerProduct
-      (\n e ->
-         if' (Graph.from e == n) 1 $
-         if' (Graph.to   e == n) (-1) $
-         0)
-      (List.delete dst $ Graph.nodes gr)
-      (Graph.edges gr)
+   Matrix.fromRowMajor $
+   Array.sample
+      (Set.mapMonotonic IdentityT $ Graph.edgeSet gr, Graph.nodeSet gr) $
+         \(IdentityT e, n) ->
+      if' (Graph.from e == n) 1 $
+      if' (Graph.to   e == n) (-1) $
+      0
 
 resistanceMatrix ::
-   (Graph.Edge edge, Ord node, Field a) =>
-   Graph edge node a nodeLabel -> Matrix a
+   (Graph.Edge edge, Ord node, Class.Floating a) =>
+   Graph edge node a nodeLabel ->
+   Triangular.Symmetric (Set (Wrap edge node)) a
 resistanceMatrix gr =
-   HMatrix.diag $ Vector.fromList $
-   Map.elems $ Graph.edgeLabels gr
+   Triangular.diagonal MatrixShape.RowMajor $
+   Array.fromMap $ Map.mapKeysMonotonic IdentityT $ Graph.edgeLabels gr
 
 fullMatrix ::
-   (Graph.Edge edge, Ord node, Field a) =>
-   Graph edge node a nodeLabel -> node -> node -> Matrix a
-fullMatrix gr src dst =
-   let currents = currentMatrix gr src dst
-       voltages = voltageMatrix gr
-   in  Matrix.fromBlocks
-          [[NC.konst 0 (1, Matrix.cols currents),
-               Matrix.asRow $ Vector.fromList $
-               map (\n -> if n==src then 1 else 0) $ Graph.nodes gr],
-           [resistanceMatrix gr, voltages],
-           [currents, NC.konst 0 (Matrix.rows currents, Matrix.cols voltages)]]
-
-rhs ::
-   (Graph.Edge edge, Ord node, Field a) =>
-   Graph edge node a nodeLabel -> node -> node -> Vector a
-rhs gr src dst =
-   mconcat
-      [NC.konst 0 1,
-       NC.konst 0 (length (Graph.edges gr)),
-       Vector.fromList $
-       map (\n -> if n==src then 1 else 0) $
-       List.delete dst $ Graph.nodes gr]
-
-
-solution ::
-   (Graph.Edge edge, Ord node, Field a) =>
-   Graph edge node a nodeLabel -> node -> node -> Vector a
-solution gr src dst =
-   fullMatrix gr src dst <\> rhs gr src dst
+   (Graph.Edge edge, Ord node, Class.Floating a) =>
+   Graph edge node a nodeLabel -> node ->
+   Triangular.Symmetric (():+:(Set (Wrap edge node):+:Set node)) a
+fullMatrix gr src =
+   let vm = voltageMatrix gr
+       mainMatrix =
+         Triangular.stackSymmetric
+            (resistanceMatrix gr)
+            vm
+            (Vector.constant
+               (MatrixShape.symmetric MatrixShape.RowMajor $ Matrix.width vm) 0)
+   in Triangular.stackSymmetric
+         (Triangular.symmetricFromList MatrixShape.RowMajor () [0])
+         (Matrix.singleRow MatrixShape.RowMajor $
+            Vector.unit (Triangular.size mainMatrix) (Right src))
+         mainMatrix
 
 resistance ::
-   (Graph.Edge edge, Ord node, Field a) =>
+   (Graph.Edge edge, Ord node, Class.Floating a) =>
    Graph edge node a nodeLabel -> node -> node -> a
 resistance gr src dst =
-   solution gr src dst
-   `NC.atIndex`
-   (length (Graph.edges gr) + length (takeWhile (dst/=) $ Graph.nodes gr))
+   let m = fullMatrix gr src
+       ix = Right (Right dst)
+   in negate $
+      Matrix.unliftColumn MatrixShape.ColumnMajor
+         (Triangular.solve m)
+         (Vector.unit (Triangular.size m) ix)
+            ! ix