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.1
+Version:             0.1.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.1
+  Tag:         0.1.0.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/linear-circuit
 
@@ -27,9 +27,9 @@
   Exposed-Modules:
     Math.LinearCircuit
   Build-Depends:
-    lapack >=0.2.2 && <0.3,
+    lapack >=0.3 && <0.4,
     netlib-ffi >=0.1.1 && <0.2,
-    comfort-array >=0.3.1 && <0.4,
+    comfort-array >=0.4 && <0.5,
     comfort-graph >=0.0 && <0.1,
     containers >=0.4 && <0.7,
     transformers >=0.4 && <0.6,
diff --git a/src/Math/LinearCircuit.hs b/src/Math/LinearCircuit.hs
--- a/src/Math/LinearCircuit.hs
+++ b/src/Math/LinearCircuit.hs
@@ -9,59 +9,50 @@
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.Netlib.Class as Class
+import Numeric.LAPACK.Matrix.Triangular ((#%%%#))
+import Numeric.LAPACK.Matrix ((#\|))
 
+import qualified Data.Array.Comfort.Boxed as BoxedArray
 import qualified Data.Array.Comfort.Storable as Array
 import Data.Array.Comfort.Storable ((!))
-import Data.Array.Comfort.Shape ((:+:))
+import Data.Array.Comfort.Shape ((:+:)((:+:)))
 
 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, Class.Floating a) =>
-   Graph edge node a nodeLabel ->
+   Set node -> Set (Wrap edge node) ->
    Matrix.General (Set (Wrap edge node)) (Set node) a
-voltageMatrix 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, Class.Floating a) =>
-   Graph edge node a nodeLabel ->
-   Triangular.Symmetric (Set (Wrap edge node)) a
-resistanceMatrix gr =
-   Triangular.diagonal MatrixShape.RowMajor $
-   Array.fromMap $ Map.mapKeysMonotonic IdentityT $ Graph.edgeLabels gr
+voltageMatrix nodes =
+   Matrix.fromRowArray nodes .
+   fmap
+      (\(IdentityT e) ->
+         Array.fromAssociations 0 nodes
+            [(Graph.from e, 1), (Graph.to e, -1)]) .
+   BoxedArray.indices
 
 fullMatrix ::
    (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
+   let edges = Map.mapKeysMonotonic IdentityT $ Graph.edgeLabels gr
+       nodes = Graph.nodeSet gr
+       order = MatrixShape.RowMajor
+       symmetricZero sh = Matrix.zero $ MatrixShape.symmetric order sh
+       unit = Vector.unit (Map.keysSet edges :+: nodes) (Right src)
+   in (symmetricZero (), Matrix.singleRow order unit)
+      #%%%#
+      (Triangular.diagonal order $ Array.fromMap edges,
+            voltageMatrix nodes $ Map.keysSet edges)
+      #%%%#
+      symmetricZero nodes
 
 resistance ::
    (Graph.Edge edge, Ord node, Class.Floating a) =>
@@ -69,8 +60,4 @@
 resistance gr src dst =
    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
+   in - (m #\| Vector.unit (Triangular.size m) ix) ! ix