diff --git a/resistor-cube.cabal b/resistor-cube.cabal
--- a/resistor-cube.cabal
+++ b/resistor-cube.cabal
@@ -1,5 +1,5 @@
 Name:                resistor-cube
-Version:             0.0.0.4
+Version:             0.0.1
 Synopsis:            Compute total resistance of a cube of resistors
 Description:
   This is an example of how to compute the total resistance
@@ -19,7 +19,7 @@
 Cabal-Version:       >=1.10
 
 Source-Repository this
-  Tag:         0.0.0.4
+  Tag:         0.0.1
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/resistor-cube
 
@@ -30,8 +30,8 @@
 Executable resistor-cube
   Main-is:             Main.hs
   Build-Depends:
-    hmatrix >=0.16 && <0.17,
-    utility-ht >=0.0.11 && <0.1,
+    lapack >=0.2.2 && <0.3,
+    comfort-array >=0.3.1 && <0.4,
     base >=4.5 && <5
   Hs-Source-Dirs:      src
   Default-Language:    Haskell2010
diff --git a/src/Main.hs b/src/Main.hs
--- a/src/Main.hs
+++ b/src/Main.hs
@@ -2,47 +2,51 @@
 Consider a cube of resistors of equal resistance.
 What is the overall resistance from one corner to the opposite one?
 -}
+{-# LANGUAGE TypeOperators #-}
 module Main where
 
-import qualified Numeric.Container as NC
-import qualified Data.Packed.Matrix as Matrix
-import qualified Data.Packed.Vector as Vector
-import qualified Numeric.LinearAlgebra.HMatrix as HMatrix
-import Data.Packed.Matrix (Matrix)
+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 Numeric.LAPACK.Format ((##))
 
-import Control.Applicative (liftA3)
-import Control.Functor.HT (outerProduct)
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable ((!))
+import Data.Array.Comfort.Shape ((:+:)((:+:)))
 
 
-{-
-Set resistance of a primitive resistor to 1.
-This way, voltage equals current.
--}
 data Coord = C0 | C1 deriving (Eq, Ord, Show, Enum, Bounded)
 data Dim = D0 | D1 | D2 deriving (Eq, Ord, Show, Enum, Bounded)
-data Corner = Corner Coord Coord Coord deriving (Eq, Ord, Show)
-data Edge = Edge Dim Coord Coord deriving (Eq, Ord, Show)
+type Corner = (Coord,Coord,Coord)
+type Edge = (Dim,Coord,Coord)
 
-allCoords :: [Coord]
-allCoords = [minBound .. maxBound]
+type ShapeEnum = Shape.Enumeration
+type CornerShape = (ShapeEnum Coord, ShapeEnum Coord, ShapeEnum Coord)
+type EdgeShape = (ShapeEnum Dim, ShapeEnum Coord, ShapeEnum Coord)
 
-allDims :: [Dim]
-allDims = [minBound .. maxBound]
+type Matrix height width = Matrix.General height width
 
-flattenCornerIndex :: Corner -> Int
-flattenCornerIndex (Corner x y z) =
-   (fromEnum x * 2 + fromEnum y) * 2 + fromEnum z
 
-flattenEdgeIndex :: Edge -> Int
-flattenEdgeIndex (Edge d x y) =
-   (fromEnum d * 2 + fromEnum x) * 2 + fromEnum y
+cornerShape :: CornerShape
+cornerShape = (Shape.Enumeration, Shape.Enumeration, Shape.Enumeration)
 
+edgeShape :: EdgeShape
+edgeShape = (Shape.Enumeration, Shape.Enumeration, Shape.Enumeration)
 
-voltageMatrix :: Matrix Double
+
+{-
+We also need a currentMatrix that implements Kirchhoff's current law
+but it turns out to equal 'transpose voltageMatrix'.
+This makes fullMatrix symmetric.
+-}
+voltageMatrix :: Matrix EdgeShape CornerShape Double
 voltageMatrix =
-   Matrix.fromLists $
-   outerProduct
-      (\(Edge ed ex ey) c ->
+   Matrix.fromRowMajor $
+   Array.sample
+      (edgeShape, cornerShape)
+      (\((ed,ex,ey), c) ->
          let ((cx, cy), cz) = selectCornerCoords ed c
          in  if ex==cx && ey==cy
                then
@@ -50,60 +54,44 @@
                      C0 ->  1
                      C1 -> -1
                else 0)
-      (liftA3 Edge allDims allCoords allCoords)
-      (liftA3 Corner allCoords allCoords allCoords)
 
 
 selectCornerCoords :: Dim -> Corner -> ((Coord, Coord), Coord)
-selectCornerCoords ed (Corner cx cy cz) =
+selectCornerCoords ed (cx,cy,cz) =
    case ed of
       D0 -> ((cy, cz), cx)
       D1 -> ((cx, cz), cy)
       D2 -> ((cx, cy), cz)
 
 sourceCorner, destCorner :: Corner
-sourceCorner = Corner C0 C0 C0
-destCorner = Corner C1 C1 C1
+sourceCorner = (C0,C0,C0)
+destCorner = (C1,C1,C1)
 
--- almost transposed voltageMatrix with some removed rows
-currentMatrix :: Matrix Double
-currentMatrix =
-   Matrix.fromLists $
-   outerProduct
-      (\c@(Corner _ _ _) (Edge ed ex ey) ->
-         let ((cx, cy), cz) = selectCornerCoords ed c
-         in  if ex==cx && ey==cy
-               then
-                  case cz of
-                     C0 -> 1
-                     C1 -> -1
-               else 0)
-      (filter (/= sourceCorner) $
-       filter (/= destCorner) $
-       liftA3 Corner allCoords allCoords allCoords)
-      (liftA3 Edge allDims allCoords allCoords)
+resistances :: Vector.Vector EdgeShape Double
+resistances = Vector.constant edgeShape 1
 
 
-fullMatrix :: Matrix Double
+fullMatrix :: Triangular.Symmetric (():+:(EdgeShape:+:CornerShape)) Double
 fullMatrix =
-   Matrix.fromBlocks
-      [[NC.konst 0 (1,12), Matrix.asRow $ Vector.fromList $ 1 : replicate 7 0],
-       [HMatrix.ident 12, voltageMatrix],
-       [currentMatrix, NC.konst 0 (6,8)]]
+   Triangular.stackSymmetric
+      (Triangular.symmetricFromList MatrixShape.RowMajor () [0])
+      (Matrix.singleRow MatrixShape.RowMajor $
+         Vector.unit (edgeShape:+:cornerShape) (Right sourceCorner)) $
+   Triangular.stackSymmetric
+      (Triangular.diagonal MatrixShape.RowMajor resistances)
+      voltageMatrix
+      (Vector.constant
+         (MatrixShape.symmetric MatrixShape.RowMajor cornerShape) 0)
 
+
 main :: IO ()
 main = do
-   print fullMatrix
-   let [currentVec, potentialVec] =
-         Vector.takesV [12,8] $ HMatrix.null1 fullMatrix
-   let totalCurrent =
-         sum $
-         map
-            (\d -> NC.atIndex currentVec $ flattenEdgeIndex $ Edge d C0 C0)
-            allDims
-   let cornerPot c = NC.atIndex potentialVec (flattenCornerIndex c)
-   let totalVoltage = cornerPot destCorner - cornerPot sourceCorner
-   print $ totalVoltage / totalCurrent
+   fullMatrix ## "%2.f"
+   let solutionVec =
+         Matrix.unliftColumn MatrixShape.ColumnMajor
+            (Triangular.solve fullMatrix) $
+         Vector.unit (():+:(edgeShape:+:cornerShape)) (Right (Right destCorner))
+   print $ - solutionVec ! Right (Right destCorner)
 
 {-
 result: total resistance is 5/(2+2+2)