diff --git a/magico.cabal b/magico.cabal
--- a/magico.cabal
+++ b/magico.cabal
@@ -1,5 +1,5 @@
 Name:                magico
-Version:             0.0.1.2
+Version:             0.0.2
 Synopsis:            Compute solutions for Magico puzzle
 Description:
   Compute solutions for Magico puzzle:
@@ -14,7 +14,7 @@
 Cabal-Version:       >=1.10
 
 Source-Repository this
-  Tag:         0.0.1.2
+  Tag:         0.0.2
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/magico
 
@@ -25,7 +25,8 @@
 Executable magico
   Main-is:             Main.hs
   Build-Depends:
-    hmatrix >=0.16 && <0.17,
+    lapack >=0.2.2 && <0.3,
+    comfort-array >=0.3.1 && <0.4,
     transformers >=0.3 && <0.6,
     utility-ht >=0.0.11 && <0.1,
     base >=4.5 && <5
diff --git a/src/Main.hs b/src/Main.hs
--- a/src/Main.hs
+++ b/src/Main.hs
@@ -2,147 +2,138 @@
 Magico
 <http://www.spectra-verlag.de/>
 -}
+{-# 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 Numeric.LinearAlgebra.HMatrix (null1, nullspace, rank)
-import Numeric.LinearAlgebra ((<\>))
-import Data.Packed.Matrix (Matrix)
-import Data.Packed.Vector (Vector)
+import qualified Numeric.LAPACK.Singular as Singular
+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix as Matrix
+import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix ((===))
+import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Format ((##))
 
-import Text.Printf (printf)
+import qualified Data.Array.Comfort.Boxed as BoxedArray
+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 ((:+:)((:+:)))
 
-import qualified Control.Monad.Trans.State as MS
-import Control.Applicative (Applicative, liftA2, liftA3, pure, (<*>))
+import Foreign.Storable (Storable)
+
 import Control.Functor.HT (outerProduct)
 
 import qualified Data.Foldable as Fold
-import qualified Data.List as List
-import Data.Traversable (Traversable, traverse, sequenceA)
-import Data.Foldable (Foldable, foldMap)
-import Data.Monoid ((<>))
+import Data.Tuple.HT (thd3)
 
 
-data Triple a = Triple a a a deriving (Show)
-data Sums a = Sums a a (Triple a) (Triple a) deriving (Show)
-type Solution a = Triple (Triple a)
+type Triple = BoxedArray.Array (Shape.Enumeration Index)
+type Sums = BoxedArray.Array SumsShape
+type Board = BoxedArray.Array SquareShape
 
+type Matrix height width = Matrix.General height width
+type ExtMatrix = Matrix (SumsShape:+:()) SquareShape
+type SumsVector = Vector SumsShape
+type Solution = Vector SquareShape
+type SquareShape = (Shape.Enumeration Index, Shape.Enumeration Index)
 
-instance Functor Triple where
-   fmap f (Triple x y z) = Triple (f x) (f y) (f z)
 
-instance Applicative Triple where
-   pure x = Triple x x x
-   Triple fx fy fz <*> Triple x y z = Triple (fx x) (fy y) (fz z)
+data Index = I0 | I1 | I2 deriving (Eq, Ord, Enum, Bounded, Show)
+type Index2 = (Index, Index)
 
-instance Foldable Triple where
-   foldMap f (Triple x y z) = f x <> f y <> f z
+squareShape :: SquareShape
+squareShape = (Shape.Enumeration, Shape.Enumeration)
 
-instance Traversable Triple where
-   traverse f (Triple x y z) = liftA3 Triple (f x) (f y) (f z)
 
+type SumsShape =
+      (():+:()) :+: (Shape.Enumeration Index :+: Shape.Enumeration Index)
 
-instance Functor Sums where
-   fmap f (Sums diag0 diag1 horiz vert) =
-      Sums (f diag0) (f diag1) (fmap f horiz) (fmap f vert)
+sumsShape :: SumsShape
+sumsShape = (():+:()) :+: (Shape.Enumeration :+: Shape.Enumeration)
 
-instance Foldable Sums where
-   foldMap f (Sums diag0 diag1 horiz vert) =
-      f diag0 <> f diag1 <> foldMap f horiz <> foldMap f vert
+sums :: (Storable a) => a -> a -> (a,a,a) -> (a,a,a) -> SumsVector a
+sums diag0 diag1 horiz vert =
+   let vec3 (x,y,z) = Vector.fromList Shape.Enumeration [x,y,z]
+       x<+>y = Vector.append x y
+   in (Vector.singleton diag0 <+> Vector.singleton diag1)
+      <+>
+      (vec3 horiz <+> vec3 vert)
 
 
-getFromList :: MS.State [a] a
-getFromList =
-   MS.state $ \(x:xs) -> (x,xs)
-
-solFromList :: [a] -> Solution a
-solFromList =
-   MS.evalState (sequenceA $ pure $ sequenceA $ pure getFromList)
-
-data Index = I0 | I1 | I2 deriving (Eq, Ord, Enum, Show)
-type Index2 = (Index, Index)
-
 indices :: Triple Index
-indices = Triple I0 I1 I2
-
-index :: Index -> Triple a -> a
-index i (Triple x y z) =
-   case i of
-      I0 -> x
-      I1 -> y
-      I2 -> z
-
-index2 :: Index2 -> Triple (Triple a) -> a
-index2 (i,j) = index j . index i
+indices = BoxedArray.indices Shape.Enumeration
 
-flattenIndex2 :: Index2 -> Int
-flattenIndex2 (i,j) = 3 * fromEnum i + fromEnum j
+triple :: a -> a -> a -> Triple a
+triple x y z = BoxedArray.fromList Shape.Enumeration [x,y,z]
 
 sumIndices :: Sums (Triple Index2)
 sumIndices =
-   Sums
-      (Triple (I0,I0) (I1,I1) (I2,I2))
-      (Triple (I0,I2) (I1,I1) (I2,I0))
-      (outerProduct (,) indices indices)
-      (outerProduct (flip (,)) indices indices)
-
-sumOnBoard :: (Num a) => Solution a -> Triple Index2 -> a
-sumOnBoard sol is = Fold.sum $ fmap (flip index2 sol) is
+   BoxedArray.fromList sumsShape $
+      triple (I0,I0) (I1,I1) (I2,I2) :
+      triple (I0,I2) (I1,I1) (I2,I0) :
+      outerProduct (,) (Shape.indices Shape.Enumeration) indices ++
+      outerProduct (flip (,)) (Shape.indices Shape.Enumeration) indices ++
+      []
 
-sums :: (Num a) => Solution a -> Sums a
-sums sol = fmap (sumOnBoard sol) sumIndices
+sumsFromSolution :: (Num a, Storable a) => Solution a -> Sums a
+sumsFromSolution sol = fmap (Fold.sum . fmap (sol Array.!)) sumIndices
 
-fullMatrix :: Matrix Double
+fullMatrix :: Matrix SumsShape SquareShape Double
 fullMatrix =
-   Matrix.fromLists $
-   outerProduct
-      (\is j -> if Fold.elem j is then 1 else 0)
-      (Fold.toList sumIndices)
-      (liftA2 (,) (Fold.toList indices) (Fold.toList indices))
-
-removeAt :: Int -> [a] -> (a, [a])
-removeAt k xs =
-   case splitAt k xs of
-      (_, []) -> error "removeAt: index too large"
-      (leftXs, pivot:rightXs) -> (pivot, leftXs++rightXs)
-
-insertAt :: Int -> a -> [a] -> [a]
-insertAt k x xs =
-   case splitAt k xs of
-      (leftXs, rightXs) -> leftXs ++ x : rightXs
-
+   Matrix.fromRowArray squareShape $
+   fmap
+      (\is -> Array.sample squareShape (\j -> if Fold.elem j is then 1 else 0))
+      sumIndices
 
-isIntegerVector :: Vector Double -> Bool
-isIntegerVector =
-   all (\x -> abs (x - fromInteger (round x)) < 1e-7) . Vector.toList
+isIntegerVector :: Solution Double -> Bool
+isIntegerVector x =
+   (Vector.normInf $ Vector.sub x $ Array.map (fromInteger . round) x) < 1e-7
 
--- | it must be @NC.sumElements offset == 0@
-integerSolutions :: Int -> Vector Double -> Vector Double -> [Vector Double]
+-- | it must be @Vector.sum offset == 0@
+integerSolutions ::
+   Int -> Solution Double -> Solution Double -> [Solution Double]
 integerSolutions maxN offset start =
-   let maxI = NC.maxIndex offset
-       startI = NC.atIndex start maxI
-       offsetI = NC.atIndex offset maxI
+   let (maxI,offsetI) = Vector.argAbs1Maximum offset
+       startI = start!maxI
    in  filter isIntegerVector $
        map
           (\x ->
              let c = (fromIntegral x - startI) / offsetI
-             in  NC.add start $ NC.scale c offset)
+             in  Vector.mac c offset start)
           [0 .. maxN]
 
-splittedMatrix :: Index2 -> ((Int, Vector Double), Matrix Double)
-splittedMatrix givenI2 =
-   let splitPos = flattenIndex2 givenI2
-       (givenCol, remCols) = removeAt splitPos $ Matrix.toColumns fullMatrix
-   in  ((splitPos, givenCol), Matrix.fromColumns remCols)
+extendedMatrix :: Index2 -> ExtMatrix Double
+extendedMatrix givenI2 =
+   fullMatrix
+   ===
+   Matrix.singleRow MatrixShape.RowMajor (Vector.unit squareShape givenI2)
 
-ranks :: Triple (Triple Int)
-ranks = outerProduct (curry $ rank . snd . splittedMatrix) indices indices
 
-masks :: Triple (Triple (Matrix Double))
-masks = outerProduct (curry $ nullspace . snd . splittedMatrix) indices indices
+rank :: ExtMatrix Double -> Int
+rank = fst . Singular.pseudoInverseRCond 1e-7
 
+ranks :: Board Int
+ranks = BoxedArray.map (rank . extendedMatrix) $ BoxedArray.indices squareShape
+
+nullspace ::
+   ExtMatrix Double -> Matrix.General Matrix.ZeroInt SquareShape Double
+nullspace mat =
+   let (_vr,sigma,vl) = Singular.decompose mat
+       rnk = length $ takeWhile (>1e-7) $ Vector.toList sigma
+   in Matrix.dropRows rnk $
+      Matrix.mapHeight (Shape.ZeroBased . Shape.size) $ Matrix.fromFull vl
+
+masks :: Board (Matrix.General Matrix.ZeroInt SquareShape Double)
+masks =
+   BoxedArray.map (nullspace . extendedMatrix) $ BoxedArray.indices squareShape
+
+lastRow :: Square.Square SquareShape Double -> Solution Double
+lastRow mat = Matrix.takeRow mat (last $ Shape.indices $ Matrix.height mat)
+
+null1 :: ExtMatrix Double -> Solution Double
+null1 = lastRow . thd3 . Singular.decompose
+
 {-
 Vectors of the null space look like
 
@@ -162,59 +153,63 @@
 In this case, 'solve' misses some solutions,
 because it checks only one dimension, not two.
 -}
-solve :: (Index2, Int) -> Sums Int -> [Solution Int]
+solve :: (Index2, Int) -> SumsVector Int -> [Solution Int]
 solve (givenI2, given) ss =
-   let vec = Vector.fromList $ Fold.toList $ fmap fromIntegral ss
-       ((splitPos, givenCol), mat) = splittedMatrix givenI2
-   in  map (solFromList . insertAt splitPos given) $
-       filter (all (>=0)) $
-       map (map round . Vector.toList) $
-       integerSolutions (Fold.maximum ss) (null1 mat) $
-          mat <\> NC.sub vec (NC.scale (fromIntegral given) givenCol)
+   let vec =
+          Vector.append
+             (Array.map fromIntegral ss)
+             (Vector.singleton $ fromIntegral given)
+       mat = extendedMatrix givenI2
+   in  filter (all (>=0) . Vector.toList) $
+       map (Array.map round) $
+       integerSolutions (Fold.maximum $ Vector.toList ss) (null1 mat) $
+       Matrix.unliftColumn MatrixShape.ColumnMajor
+          (snd . Singular.leastSquaresMinimumNormRCond 1e-7 mat) vec
 
 printSolution :: Solution Int -> IO ()
 printSolution sol = do
-   Fold.forM_ sol $
-      putStrLn . List.intercalate " " . Fold.toList . fmap (printf "%2d")
+   Matrix.fromRowMajor (Array.map (fromIntegral :: Int -> Float) sol) ## "%2.f"
    putStrLn ""
 
 
 exampleA1 :: [Solution Int]
-exampleA1 =
-   solve ((I0,I2), 2) (Sums 5 8 (Triple 5 8 8) (Triple 9 6 6))
+exampleA1 = solve ((I0,I2), 2) (sums 5 8 (5,8,8) (9,6,6))
 
 exampleA2 :: [Solution Int]
-exampleA2 =
-   solve ((I2,I1), 3) (Sums 9 6 (Triple 8 9 7) (Triple 6 9 9))
+exampleA2 = solve ((I2,I1), 3) (sums 9 6 (8,9,7) (6,9,9))
 
 exampleA3 :: [Solution Int]
-exampleA3 =
-   solve ((I2,I1), 4) (Sums 9 4 (Triple 6 8 8) (Triple 9 7 6))
+exampleA3 = solve ((I2,I1), 4) (sums 9 4 (6,8,8) (9,7,6))
 
 exampleA4 :: [Solution Int]
-exampleA4 =
-   solve ((I1,I2), 4) (Sums 9 6 (Triple 8 9 7) (Triple 9 9 6))
+exampleA4 = solve ((I1,I2), 4) (sums 9 6 (8,9,7) (9,9,6))
 
 exampleA7 :: [Solution Int]
-exampleA7 =
-   solve ((I0,I0), 5) (Sums 8 4 (Triple 9 6 6) (Triple 9 9 3))
+exampleA7 = solve ((I0,I0), 5) (sums 8 4 (9,6,6) (9,9,3))
 
 exampleA20 :: [Solution Int]
-exampleA20 =
-   solve ((I1,I2), 4) (Sums 7 7 (Triple 9 6 8) (Triple 8 6 9))
+exampleA20 = solve ((I1,I2), 4) (sums 7 7 (9,6,8) (8,6,9))
 
 exampleB1 :: [Solution Int]
-exampleB1 =
-   solve ((I0,I1), 4) (Sums 11 11 (Triple 10 7 9) (Triple 9 9 8))
+exampleB1 = solve ((I0,I1), 4) (sums 11 11 (10,7,9) (9,9,8))
 
 exampleD7 :: [Solution Int]
-exampleD7 =
-   solve ((I0,I0), 20) (Sums 44 34 (Triple 43 40 36) (Triple 37 49 33))
+exampleD7 = solve ((I0,I0), 20) (sums 44 34 (43,40,36) (37,49,33))
 
 exampleD20 :: [Solution Int]
-exampleD20 =
-   solve ((I2,I2), 9) (Sums 37 43 (Triple 46 38 39) (Triple 49 35 39))
+exampleD20 = solve ((I2,I2), 9) (sums 37 43 (46,38,39) (49,35,39))
 
+
+analysis :: IO ()
+analysis = do
+   Fold.forM_ masks $ \maskBasis -> do
+      Fold.forM_ (Matrix.toRows maskBasis) $ \mask -> do
+         Matrix.fromRowMajor mask ## "%6.3f"
+         putStrLn ""
+      putStrLn ""
+   Matrix.fromRowMajor
+      (Array.map (fromIntegral :: Int -> Float) (Array.fromBoxed ranks))
+         ## "%2.f"
 
 main :: IO ()
 main = mapM_ printSolution exampleD7