diff --git a/Matrix/LU.hs b/Matrix/LU.hs
--- a/Matrix/LU.hs
+++ b/Matrix/LU.hs
@@ -14,8 +14,13 @@
 
 module Matrix.LU (lu, lu_solve, improve, inverse, lu_det, solve, det) where
 
+import qualified Matrix.Matrix as Matrix
+import qualified Matrix.Vector as Vector
+import qualified Data.List as List
 import Data.Array
+import Data.Ord (comparing)
 
+
 -- | LU decomposition via Crout's Algorithm
 
 -- TODO: modify for partial pivoting / permutation matrix
@@ -37,8 +42,8 @@
 
 -- | Solution to Ax=b via LU decomposition
 
--- forward is forumla (2.3.6) in NRIC, but remebering that a11=1
--- backward is forumla (2.3.7) in NRIC
+-- forward is formula (2.3.6) in NRIC, but remembering that a11=1
+-- backward is formula (2.3.7) in NRIC
 
 lu_solve :: Array (Int,Int) Double -- ^ LU(A)
          -> Array Int Double -- ^ b
@@ -107,10 +112,37 @@
 
 -- | determinant using original matrix
 
+{-
+It is based on LU decomposition without singularity check
+and thus returns NaN instead of zero if the matrix is singular.
+-}
+_det :: Array (Int,Int) Double -- ^ A
+    -> Double -- ^ det(A)
+
+_det a = (lu_det . lu) a
+
+{- |
+Determinant computation by implicit LU decomposition with permutations.
+-}
 det :: Array (Int,Int) Double -- ^ A
     -> Double -- ^ det(A)
 
-det a = (lu_det . lu) a
+det a =
+   if rangeSize (bounds a) == 0
+     then 1
+     else
+         let ((m0,n0), (m1,n1)) = bounds a
+             v = Matrix.getColumn n0 a
+             (maxi,maxv) = List.maximumBy (comparing (abs . snd)) $ assocs v
+             reduced =
+                ixmap ((m0,n0), (pred m1, pred n1))
+                   (\(i,j) -> (if i<maxi then i else succ i, succ j)) $
+                Vector.sub a $ Matrix.outer v $
+                Vector.scale (recip maxv) $ Matrix.getRow maxi a
+             sign = if even (rangeSize (m0,maxi)-1) then 1 else -1
+             pivot = a!(maxi,n0)
+         in  if pivot == 0 then 0 else sign * pivot * det reduced
+
 
 -------------------------------------------------------------------------------
 -- tests
diff --git a/Matrix/Matrix.hs b/Matrix/Matrix.hs
--- a/Matrix/Matrix.hs
+++ b/Matrix/Matrix.hs
@@ -14,51 +14,109 @@
 
 module Matrix.Matrix where
 
+import Matrix.Vector (generate)
 import Data.Array
 import Data.Complex
 
 -- | Matrix-matrix multiplication: A x B = C
 
-mm_mult :: (Ix a, Integral a, Num b) => Array (a,a) b -- ^ A
-	-> Array (a,a) b -- ^ B
-	-> Array (a,a) b -- ^ C
+mm_mult :: (Ix i, Ix j, Ix k, Num a) => Array (i,j) a -- ^ A
+	-> Array (j,k) a -- ^ B
+	-> Array (i,k) a -- ^ C
 
-mm_mult a b = if ac /= br
+mm_mult a b = if (ac0,ac1) /= (br0,br1)
 	      then error "mm_mult: inside dimensions inconsistent"
-	      else array bnds [ ((i,j), mult i j) | (i,j) <- range bnds ]
-    where mult i j = sum [ a!(i,k) * b!(k,j) | k <- [1..ac] ]
-	  ((_,_),(ar,ac)) = bounds a
-	  ((_,_),(br,bc)) = bounds b
-	  bnds = ((1,1),(ar,bc))
+	      else generate ((ar0,bc0),(ar1,bc1)) $ \(i,j) ->
+			sum [ a!(i,k) * b!(k,j) | k <- range (ac0,ac1) ]
+    where ((ar0,ac0),(ar1,ac1)) = bounds a
+	  ((br0,bc0),(br1,bc1)) = bounds b
 
 -- | Matrix-vector multiplication: A x b = c
 
-mv_mult :: (Ix a, Integral a, Num b) => Array (a,a) b -- ^ A
-	-> Array a b -- ^ b
-	-> Array a b -- ^ c
+mv_mult :: (Ix i, Ix j, Num a) => Array (i,j) a -- ^ A
+	-> Array j a -- ^ b
+	-> Array i a -- ^ c
 
-mv_mult a b = if ac /= br
+mv_mult a b = if (ac0,ac1) /= bounds b
 	      then error "mv_mult: dimensions inconsistent"
-	      else array bnds [ (i, mult i) | i <- range bnds ]
-    where mult i = sum [ a!(i,k) * b!(k) | k <- [1..ac] ]
-	  ((_,_),(ar,ac)) = bounds a
-	  (_,br) = bounds b
-	  bnds = (1,ar)
+	      else generate (ar0,ar1) $ \i ->
+			sum [ a!(i,k) * bk | (k,bk) <- assocs b ]
+    where ((ar0,ac0),(ar1,ac1)) = bounds a
 
 -- | Transpose of a matrix
 
-m_trans :: (Ix a, Integral a, Num b) => Array (a,a) b -- ^ A
-	-> Array (a,a) b -- ^ A^T
+m_trans :: (Ix i, Ix j, Num a) => Array (i,j) a -- ^ A
+	-> Array (j,i) a -- ^ A^T
 
-m_trans a = array bnds [ ((i,j), a!(j,i)) | (i,j) <- range bnds ]
-    where (_,(m,n)) = bounds a
-	  bnds = ((1,1),(n,m))
+m_trans a = generate ((n0,m0),(n1,m1)) $ \(i,j) -> a!(j,i)
+    where ((m0,n0),(m1,n1)) = bounds a
 
 -- | Hermitian transpose (conjugate transpose) of a matrix
 
-m_hermit :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -- ^ A
-	 -> Array (a,a) (Complex b) -- ^ A^H
+m_hermit :: (Ix i, Ix j, RealFloat a) => Array (i,j) (Complex a) -- ^ A
+	 -> Array (j,i) (Complex a) -- ^ A^H
 
-m_hermit a = array bnds [ ((i,j), conjugate (a!(j,i))) | (i,j) <- range bnds ]
-    where (_,(m,n)) = bounds a
-	  bnds = ((1,1),(n,m))
+m_hermit a = generate ((n0,m0),(n1,m1)) $ \(i,j) -> conjugate (a!(j,i))
+    where ((m0,n0),(m1,n1)) = bounds a
+
+
+columnBounds :: (Ix i, Ix j) => Array (i,j) a -> (i,i)
+columnBounds a =
+   let ((m0,_n0), (m1,_n1)) = bounds a
+   in  (m0,m1)
+
+rowBounds :: (Ix i, Ix j) => Array (i,j) a -> (j,j)
+rowBounds a =
+   let ((_m0,n0), (_m1,n1)) = bounds a
+   in  (n0,n1)
+
+getColumn :: (Ix i, Ix j) => j -> Array (i,j) e -> Array i e
+getColumn j a = ixmap (columnBounds a) (\k -> (k,j)) a
+
+getRow :: (Ix i, Ix j) => i -> Array (i,j) e -> Array j e
+getRow k a = ixmap (rowBounds a) (\j -> (k,j)) a
+
+toColumns :: (Ix i, Ix j) => Array (i,j) a -> [Array i a]
+toColumns a = map (flip getColumn a) $ range $ rowBounds a
+
+toRows :: (Ix i, Ix j) => Array (i,j) a -> [Array j a]
+toRows a = map (flip getRow a) $ range $ columnBounds a
+
+
+{- |
+We need the bounds of the row indices for empty input lists.
+-}
+fromColumns :: (Ix i) => (i,i) -> [Array i a] -> Array (i,Int) a
+fromColumns bnds@(m0,m1) columns =
+   if all ((bnds==) . bounds) columns
+     then array ((m0,0), (m1, length columns - 1)) $ concat $
+          zipWith
+            (\k -> map (\(i,a) -> ((i,k),a)) . assocs)
+            [0..] columns
+     else error "Matrix.fromColumns: column bounds mismatch"
+
+fromRows :: (Ix j) => (j,j) -> [Array j a] -> Array (Int,j) a
+fromRows bnds@(n0,n1) rows =
+   if all ((bnds==) . bounds) rows
+     then array ((0,n0), (length rows - 1, n1)) $ concat $
+          zipWith
+            (\k -> map (\(i,a) -> ((k,i),a)) . assocs)
+            [0..] rows
+     else error "Matrix.fromRows: row bounds mismatch"
+
+
+
+outer :: (Ix i, Ix j, Num a) => Array i a -> Array j a -> Array (i,j) a
+outer x y =
+   let (m0,m1) = bounds x
+       (n0,n1) = bounds y
+   in  array ((m0,n0), (m1,n1)) $ do
+         (i,xi) <- assocs x
+         (j,yj) <- assocs y
+         return ((i,j), xi*yj)
+
+inner :: (Ix i, Num a) => Array i a -> Array i a -> a
+inner x y =
+   if bounds x == bounds y
+     then sum $ zipWith (*) (elems x) (elems y)
+     else error "inner: dimensions mismatch"
diff --git a/Matrix/QR/Givens.hs b/Matrix/QR/Givens.hs
new file mode 100644
--- /dev/null
+++ b/Matrix/QR/Givens.hs
@@ -0,0 +1,160 @@
+module Matrix.QR.Givens  (
+   leastSquares,
+   decompose, solve, det,
+   Rotation, rotateVector,
+   Upper, solveUpper, detUpper,
+   ) where
+
+import qualified Matrix.Sparse as Sparse
+import DSP.Basic (toMaybe, (^!))
+
+import Control.Monad (mfilter)
+
+import qualified Data.Foldable as Fold
+import qualified Data.List as List
+import qualified Data.Map as Map
+import Data.Map (Map)
+import Data.Array
+         (Array, Ix, array, bounds, elems, rangeSize, range, (!), (//), )
+
+
+data Rotation i a = Rotation (i,i) (a,a)
+   deriving Show
+
+data Upper i j a = Upper ((i,j), (i,j)) (Map i (Map j a))
+   deriving Show
+
+{- |
+The decomposition routine is pretty simple.
+It does not try to minimize fill-up by a clever ordering of rotations.
+However, for banded matrices it will work as expected.
+-}
+decompose :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>
+      Sparse.Matrix i j a -- ^ A
+   -> ([Rotation i a], Upper i j a) -- ^ QR(A)
+decompose a =
+   (\(qs,rows) -> (concat qs, Upper (Sparse.bounds a) $ Map.fromList rows)) .
+   unzip .
+   List.unfoldr
+      (\a0 ->
+         let bnds@((m0,_), _) = Sparse.bounds a0
+             (q,a1) = step a0
+         in  toMaybe (not $ emptyRange bnds) $
+               ((q, (m0, Sparse.getRow m0 a1)), submatrix a1))
+    $ a
+
+-- cf. QR.Householder
+emptyRange :: (Ix i) => (i,i) -> Bool
+emptyRange = null . range
+
+-- | assumes that the first column is empty except the top-most element
+submatrix ::
+   (Ord i, Enum i, Ord j, Enum j) => Sparse.Matrix i j a -> Sparse.Matrix i j a
+submatrix a =
+   let ((m0,n0), mn1) = Sparse.bounds a
+   in  Sparse.fromRows ((succ m0, succ n0), mn1) $
+       Map.delete m0 $ Sparse.toRows a
+
+step ::
+   (Ord i, Ord j, RealFloat a) =>
+   Sparse.Matrix i j a -> ([Rotation i a], Sparse.Matrix i j a)
+step a =
+   let bnds@((m0,n0), _) = Sparse.bounds a
+       rows = Sparse.toRows a
+       topRow = Map.findWithDefault Map.empty m0 rows
+   in  (\((xi,xrem), finalRows) ->
+         (Map.elems $ Map.mapMaybe fst finalRows,
+          Sparse.fromRows bnds $ Map.insert m0 (Map.insert n0 xi xrem) $
+            fmap snd finalRows)) $
+       Map.mapAccumWithKey
+         (\(xi,xrem) mi yrow ->
+            let yrem = Map.delete n0 yrow
+                jrot = Just . Rotation (m0,mi)
+            in  case mfilter (0/=) $ Map.lookup n0 yrow of
+                  Nothing -> ((xi,xrem),(Nothing,yrem))
+                  Just yi ->
+                     if xi==0
+                        then ((yi,yrem), (jrot (0,-1), fmap negate xrem))
+                        else
+                           let rot = rotationFromXY (xi,yi)
+                               (rx,ry) = rotateRows rot (xrem, yrem)
+                           in  ((fst $ rotateSingle rot (xi,yi), rx),
+                                (jrot rot, ry)))
+         (Map.findWithDefault 0 n0 topRow, Map.delete n0 topRow)
+         (Map.delete m0 rows)
+
+-- | The argument must not be (0,0).
+rotationFromXY :: (RealFloat a) => (a,a) -> (a,a)
+rotationFromXY (x,y) =
+   if abs x > abs y
+     then let q = y/x; k = recipNorm q in (k,-q*k)
+     else let q = x/y; k = recipNorm q in (q*k,-k)
+
+recipNorm :: Floating a => a -> a
+recipNorm q = recip $ sqrt (1+q^!2)
+
+rotateSingle :: (Num a) => (a,a) -> (a,a) -> (a,a)
+rotateSingle (c,s) (x,y) = (c*x-s*y, s*x+c*y)
+
+rotateRows ::
+   (Ord j, Num a) => (a,a) -> (Map j a, Map j a) -> (Map j a, Map j a)
+rotateRows (c,s) (xs,ys) =
+   let rs =
+         Map.intersectionWith (curry $ rotateSingle (c,s)) xs ys
+         `Map.union`
+         fmap (\x -> ( c*x, s*x)) (Map.difference xs ys)
+         `Map.union`
+         fmap (\y -> (-s*y, c*y)) (Map.difference ys xs)
+
+   in  (fmap fst rs, fmap snd rs)
+
+rotateVector :: (Ix i, Num a) => Rotation i a -> Array i a -> Array i a
+rotateVector (Rotation (i0,i1) cs) x =
+   let (y0,y1) = rotateSingle cs (x!i0,x!i1)
+   in  x // [(i0,y0),(i1,y1)]
+
+
+{- |
+Assumes that 'Upper' matrix is at least as high as wide
+and that it has full rank.
+-}
+solveUpper ::
+   (Ix i, Ix j, Fractional a) => Upper i j a -> Array i a -> Array j a
+solveUpper (Upper ((m0,n0), (m1,n1)) rows0) b =
+   if bounds b == (m0,m1)
+     then
+         array (n0,n1) $ Map.toList $
+         foldr
+            (\(row,bi) xs ->
+               let ((j,a),as) = Map.deleteFindMin row
+               in  Map.insert j
+                     ((bi - Fold.sum (Map.intersectionWith (*) as xs)) / a) xs)
+            Map.empty
+            (zip (Map.elems rows0) (elems b))
+     else error "solveUpper: vertical bounds mismatch"
+
+solve ::
+   (Ix i, Ix j, Fractional a) =>
+   ([Rotation i a], Upper i j a) -> Array i a -> Array j a
+solve (qs, u) b = solveUpper u $ foldl (flip rotateVector) b qs
+
+{- |
+Solve a sparse overconstrained linear problem, i.e. minimize @||Ax-b||@.
+@A@ must have dimensions @m x n@ with @m>=n@ and it must have full-rank.
+None of these conditions is checked.
+-}
+leastSquares ::
+   (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>
+   Sparse.Matrix i j a -> Array i a -> Array j a
+leastSquares = solve . decompose
+
+
+detUpper ::
+   (Ix i, Ix j, Fractional a) => Upper i j a -> a
+detUpper (Upper ((_m0,n0), (_m1,n1)) rows) =
+   if rangeSize (n0,n1) == Map.size rows
+     then product $ map (snd . Map.findMin) $ Map.elems rows
+     else 0
+
+det :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Sparse.Matrix i j a -> a
+det = detUpper . snd . decompose
diff --git a/Matrix/QR/Householder.hs b/Matrix/QR/Householder.hs
new file mode 100644
--- /dev/null
+++ b/Matrix/QR/Householder.hs
@@ -0,0 +1,152 @@
+module Matrix.QR.Householder (
+   leastSquares,
+   decompose, solve, det,
+   Reflection, reflectMatrix, reflectVector,
+   Upper, matrixFromUpper, solveUpper, detUpper,
+   ) where
+
+import Matrix.Matrix (mv_mult, m_trans, getRow, getColumn, inner, outer)
+import Matrix.Vector (sub, scale, norm)
+import DSP.Basic (toMaybe)
+
+import qualified Data.List as List
+import Data.Array
+         (Array, Ix, bounds, elems, range, rangeSize,
+          accum, accumArray, assocs, ixmap, listArray, (!), (//), )
+
+
+decompose :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>
+      Array (i,j) a -- ^ A
+   -> ([Reflection i a], Upper i j a) -- ^ QR(A)
+decompose a =
+   (\(qs,rows) -> (qs, Upper (bounds a) rows)) .
+   unzip .
+   List.unfoldr
+      (\a0 ->
+         let bnds@((m0,_), _) = bounds a0
+         in  toMaybe (not $ emptyRange bnds) $
+               let (q,a1) = step a0
+               in  ((q, getRow m0 a1), submatrix a1))
+    $ a
+
+emptyRange :: (Ix i) => (i,i) -> Bool
+emptyRange = null . range
+
+step ::
+   (Ix i, Ix j, RealFloat a) =>
+   Array (i,j) a -> (Reflection i a, Array (i,j) a)
+step a =
+   let (m0,n0) = fst $ bounds a
+       z = getColumn n0 a
+       sign x = if x<0 then -1 else 1
+       q = reflection $ accum (+) z [(m0, sign(z!m0) * norm z)]
+   in  (q, reflectMatrix q a)
+
+{-
+Submatrices with only Ix constrained indices would not work,
+because we cannot reduce a two-dimensional array by only one element.
+-}
+submatrix :: (Ix i, Enum i, Ix j, Enum j) => Array (i,j) e -> Array (i,j) e
+submatrix a =
+   let ((m0,n0), (m1,n1)) = bounds a
+   in  ixmap ((succ m0, succ n0), (m1,n1)) id a
+
+
+data Upper i j a = Upper ((i,j), (i,j)) [Array j a]
+
+matrixFromUpper :: (Ix i, Ix j, Num a) => Upper i j a -> Array (i,j) a
+matrixFromUpper (Upper bnds@((m0,_n0), (m1,_n1)) rows) =
+   accumArray (const id) 0 bnds $ concat $
+   zipWith (\k -> map (\(j,a) -> ((k,j),a)) . assocs) (range (m0,m1)) rows
+
+
+newtype Reflection i a = Reflection (Array i a)
+
+reflection :: (Ix i, Floating a) => Array i a -> Reflection i a
+reflection v =
+   let normv = norm v
+   in  Reflection $ fmap (/ ((1-signum normv) + normv)) v
+
+reflectMatrixFull ::
+   (Ix i, Ix j, Num a) => Reflection i a -> Array (i,j) a -> Array (i,j) a
+reflectMatrixFull (Reflection v) a =
+   sub a $ scale 2 $ outer v $ mv_mult (m_trans a) v
+
+reflectMatrix ::
+   (Ix i, Ix j, Num a) => Reflection i a -> Array (i,j) a -> Array (i,j) a
+reflectMatrix q@(Reflection v) a =
+   let (k0,k1) = bounds v
+       ((m0,n0), (m1,n1)) = bounds a
+       bnds = ((k0,n0),(k1,n1))
+   in  case (compare k0 m0, compare k1 m1) of
+         (EQ,EQ) -> reflectMatrixFull q a
+         (LT,_) -> error "reflectMatrix: lower reflection dimension too small"
+         (_,GT) -> error "reflectMatrix: upper reflection dimension too big"
+         _ -> replaceSubArray a $ reflectMatrixFull q $ subArray bnds a
+
+
+reflectVectorFull :: (Ix i, Num a) => Reflection i a -> Array i a -> Array i a
+reflectVectorFull (Reflection v) a = sub a $ scale (2 * inner v a) v
+
+reflectVector :: (Ix i, Num a) => Reflection i a -> Array i a -> Array i a
+reflectVector q@(Reflection v) a =
+   let bnds@(k0,k1) = bounds v
+       (m0,m1) = bounds a
+   in  case (compare k0 m0, compare k1 m1) of
+         (EQ,EQ) -> reflectVectorFull q a
+         (LT,_) -> error "reflectVector: lower reflection dimension too small"
+         (_,GT) -> error "reflectVector: upper reflection dimension too big"
+         _ -> replaceSubArray a $ reflectVectorFull q $ subArray bnds a
+
+
+subArray :: (Ix i) => (i,i) -> Array i a -> Array i a
+subArray bnds = ixmap bnds id
+
+replaceSubArray :: (Ix i) => Array i a -> Array i a -> Array i a
+replaceSubArray x y = x // assocs y
+
+
+{- |
+Assumes that 'Upper' matrix is at least as high as wide
+and that it has full rank.
+-}
+solveUpper ::
+   (Ix i, Ix j, Fractional a) => Upper i j a -> Array i a -> Array j a
+solveUpper (Upper ((m0,n0), (m1,n1)) rs0) b =
+   if bounds b == (m0,m1)
+     then
+         listArray (n0,n1) $
+         foldr
+            (\(r,bi) xs ->
+               let (a:as) = elems r
+               in  (bi - sum (zipWith (*) as xs)) / a : xs)
+            []
+            (zip rs0 (elems b))
+     else error "solveUpper: vertical bounds mismatch"
+
+solve ::
+   (Ix i, Ix j, Fractional a) =>
+   ([Reflection i a], Upper i j a) -> Array i a -> Array j a
+solve (qs, u) b = solveUpper u $ foldl (flip reflectVector) b qs
+
+{- |
+Solve an overconstrained linear problem, i.e. minimize @||Ax-b||@.
+@A@ must have dimensions @m x n@ with @m>=n@ and it must have full-rank.
+None of these conditions is checked.
+-}
+leastSquares ::
+   (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>
+   Array (i,j) a -> Array i a -> Array j a
+leastSquares = solve . decompose
+
+
+detUpper :: (Ix i, Ix j, Fractional a) => Upper i j a -> a
+detUpper (Upper ((_m0,n0), (_m1,n1)) rs) =
+   if rangeSize (n0,n1) == length rs
+     then product $ map (head . elems) rs
+     else 0
+
+det :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Array (i,j) a -> a
+det a =
+   let (qs,u) = decompose a
+   in  (if even (length qs) then 1 else -1) * detUpper u
diff --git a/Matrix/Sparse.hs b/Matrix/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/Matrix/Sparse.hs
@@ -0,0 +1,84 @@
+module Matrix.Sparse (
+   Matrix,
+   bounds,
+   fromMap,
+   fromRows,
+   fromColumns,
+   fromDense,
+   toRows,
+   toColumns,
+   toDense,
+   getRow,
+   getColumn,
+   mulVector,
+   ) where
+
+import qualified Matrix.Vector as Vector
+import qualified Data.Foldable as Fold
+import qualified Data.Map as Map
+import qualified Data.Array as Array
+import Data.Map (Map)
+import Data.Array (Array, Ix, accumArray, (!))
+
+
+data Matrix i j a = Matrix ((i,j), (i,j)) (Map i (Map j a))
+   deriving Show
+
+instance Functor (Matrix i j) where
+   fmap f (Matrix bnds m) = Matrix bnds $ fmap (fmap f) m
+
+
+bounds :: Matrix i j a -> ((i,j), (i,j))
+bounds (Matrix bnds _) = bnds
+
+fromMap :: (Ord i, Ord j) => ((i,j), (i,j)) -> Map (i,j) a -> Matrix i j a
+fromMap bnds =
+   Matrix bnds . Map.fromListWith Map.union .
+   map (\((i,j),a) -> (i, Map.singleton j a)) . Map.toList
+
+fromRows ::
+   (Ord i, Ord j) => ((i,j), (i,j)) -> Map i (Map j a) -> Matrix i j a
+fromRows = Matrix
+
+fromColumns ::
+   (Ord i, Ord j) => ((i,j), (i,j)) -> Map j (Map i a) -> Matrix i j a
+fromColumns bnds = Matrix bnds . flipMap
+
+fromDense :: (Ix i, Ix j) => Array (i,j) a -> Matrix i j a
+fromDense a = fromMap (Array.bounds a) $ Map.fromList $ Array.assocs a
+
+
+toRows :: (Ord i, Ord j) => Matrix i j a -> Map i (Map j a)
+toRows (Matrix _bnds rows) = rows
+
+toColumns :: (Ord i, Ord j) => Matrix i j a -> Map j (Map i a)
+toColumns (Matrix _bnds rows) = flipMap rows
+
+toDense :: (Ix i, Ix j, Num a) => Matrix i j a -> Array (i,j) a
+toDense (Matrix bnds a) =
+   accumArray (const id) 0 bnds $ Fold.fold $
+   Map.mapWithKey (\i -> map (\(j,e) -> ((i,j),e)) .  Map.toList) a
+
+
+-- cf. comfort-graph:Graph.Comfort.Map.flip
+flipMap :: (Ord i, Ord j) => Map i (Map j a) -> Map j (Map i a)
+flipMap =
+   Map.unionsWith (Map.unionWith (error $ "Map.flip: duplicate key")) .
+   Map.elems . Map.mapWithKey (fmap . Map.singleton)
+
+
+getRow :: (Ord i, Ord j) => i -> Matrix i j a -> Map j a
+getRow i (Matrix _ rows) = Map.findWithDefault Map.empty i rows
+
+getColumn :: (Ord i, Ord j) => j -> Matrix i j a -> Map i a
+getColumn j (Matrix _ rows) = Map.mapMaybe (Map.lookup j) rows
+
+
+mulVector :: (Ix i, Ix j, Num a) => Matrix i j a -> Array j a -> Array i a
+mulVector a@(Matrix ((m0,n0), (m1,n1)) _) v =
+   if (n0,n1) == Array.bounds v
+     then Vector.generate (m0,m1) $ flip mulRowVector v . flip getRow a
+     else error "Sparse.mulVector: dimensions mismatch"
+
+mulRowVector :: (Ix j, Num a) => Map j a -> Array j a -> a
+mulRowVector row v = Fold.sum $ Map.mapWithKey (\j x -> x * v!j) row
diff --git a/Matrix/Vector.hs b/Matrix/Vector.hs
new file mode 100644
--- /dev/null
+++ b/Matrix/Vector.hs
@@ -0,0 +1,35 @@
+module Matrix.Vector where
+
+import DSP.Basic ((^!))
+import Data.Array
+         (Array, Ix, bounds, elems, range, array, assocs, listArray, (!), )
+
+
+generate :: (Ix i) => (i,i) -> (i -> a) -> Array i a
+generate bnds f = array bnds $ map (\i -> (i, f i)) $ range bnds
+
+fromList :: [a] -> Array Int a
+fromList xs = listArray (0, length xs - 1) xs
+
+toList :: Array Int a -> [a]
+toList = elems
+
+
+norm :: (Ix i, Floating a) => Array i a -> a
+norm = sqrt . sum . elems . fmap (^!2)
+
+scale :: (Ix i, Num a) => a -> Array i a -> Array i a
+scale x = fmap (x*)
+
+
+lift2 :: (Ix i) => (a -> b -> c) -> Array i a -> Array i b -> Array i c
+lift2 f x y =
+   if bounds x == bounds y
+     then array (bounds x) [ (k, f xk (y!k)) | (k, xk) <- assocs x ]
+     else error "Vector.lift2: matrix dimensions mismatch"
+
+add :: (Ix i, Num a) => Array i a -> Array i a -> Array i a
+add = lift2 (+)
+
+sub :: (Ix i, Num a) => Array i a -> Array i a -> Array i a
+sub = lift2 (-)
diff --git a/dsp.cabal b/dsp.cabal
--- a/dsp.cabal
+++ b/dsp.cabal
@@ -1,5 +1,5 @@
 Name:             dsp
-Version:          0.2.3.1
+Version:          0.2.4
 License:          GPL
 License-File:     COPYING
 Copyright:        Matt Donadio, 2003
@@ -11,14 +11,13 @@
 Synopsis:         Haskell Digital Signal Processing
 Description:      Digital Signal Processing, Fourier Transform, Linear Algebra, Interpolation
 Category:         Sound, Math
-
-Tested-With:      GHC==6.4.1, GHC==6.8.2
-Tested-With:      GHC==7.4.2, GHC==7.6.3
+Tested-With:      GHC==7.4.2, GHC==7.6.3, GHC==7.8.4
+Tested-With:      GHC==8.0.2, GHC==8.2.1
 Cabal-Version:    >=1.8
 Build-Type:       Simple
 
 Source-Repository this
-  Tag:         0.2.3.1
+  Tag:         0.2.4
   Type:        darcs
   Location:    http://code.haskell.org/~thielema/dsp/
 
@@ -26,22 +25,16 @@
   Type:        darcs
   Location:    http://code.haskell.org/~thielema/dsp/
 
-Flag splitBase
-  Description: Choose the new smaller, split-up base package.
-
 Flag buildExamples
   Description: Build demo executables
-  Default: True
+  Default: False
 
 Library
-  If flag(splitBase)
-    Build-Depends:
-      array >=0.1 && <0.6,
-      random >=1.0 && <1.2,
-      base >= 2 && <5
-  Else
-    Build-Depends:
-      base >=1.0 && <2
+  Build-Depends:
+    containers >=0.3 && <0.6,
+    array >=0.1 && <0.6,
+    random >=1.0 && <1.2,
+    base >= 2 && <5
 
   GHC-Options:      -Wall
   Exposed-modules:
@@ -75,6 +68,7 @@
     DSP.Filter.IIR.Matchedz
     DSP.Filter.IIR.Prony
     DSP.Filter.IIR.Transform
+    DSP.Filter.IIR.Cookbook
     DSP.Flowgraph
     DSP.Multirate.CIC
     DSP.Multirate.Halfband
@@ -85,8 +79,12 @@
     DSP.Window
     Matrix.Cholesky
     Matrix.LU
+    Matrix.QR.Householder
+    Matrix.QR.Givens
     Matrix.Levinson
     Matrix.Matrix
+    Matrix.Vector
+    Matrix.Sparse
     Matrix.Simplex
     Numeric.Approximation.Chebyshev
     Numeric.Random.Distribution.Binomial
@@ -123,9 +121,21 @@
     Polynomial.Chebyshev
     Polynomial.Maclaurin
     Polynomial.Roots
-    DSP.Filter.IIR.Cookbook
   Other-Modules:
     Numeric.Transform.Fourier.Eigensystem
+
+Test-Suite dsp-test
+  Type: exitcode-stdio-1.0
+  Main-Is: Main.hs
+  Other-Modules: Test.Matrix.QR
+  Hs-Source-Dirs: test
+  GHC-Options: -Wall
+  Build-Depends:
+    QuickCheck >=2.5 && <3,
+    dsp,
+    containers,
+    array,
+    base
 
 Executable dsp-demo-article
   Main-Is: Article.hs
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,15 @@
+module Main where
+
+import qualified Test.Matrix.QR as QR
+
+
+prefix :: String -> [(String, IO ())] -> [(String, IO ())]
+prefix msg =
+   map (\(str,test) -> (msg ++ "." ++ str, test))
+
+main :: IO ()
+main =
+   mapM_ (\(msg,io) -> putStr (msg++": ") >> io) $
+   concat $
+      prefix "Matrix.QR" QR.tests :
+      []
diff --git a/test/Test/Matrix/QR.hs b/test/Test/Matrix/QR.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Matrix/QR.hs
@@ -0,0 +1,164 @@
+module Test.Matrix.QR where
+
+import qualified Matrix.QR.Householder as Householder
+import qualified Matrix.QR.Givens as Givens
+import qualified Matrix.LU as LU
+import qualified Matrix.Vector as Vector
+import qualified Matrix.Sparse as Sparse
+import Matrix.Matrix (m_trans, mm_mult, mv_mult)
+import DSP.Basic ((^!))
+
+import Control.Applicative (liftA2, (<$>))
+
+import qualified Data.Foldable as Fold
+import qualified Data.Map as Map
+import Data.Array
+
+import qualified Test.QuickCheck as QC
+
+
+doubleArray :: (Ix i) => Array i Int -> Array i Double
+doubleArray = fmap fromIntegral
+
+gramian :: Array (Int,Int) Double -> Double
+gramian m = LU.det (m_trans m `mm_mult` m)
+
+fullRank :: Array (Int,Int) Int -> Bool
+fullRank m = round (gramian $ doubleArray m) /= (0::Integer)
+
+arbitraryIntArray :: (Ix i) => (i,i) -> QC.Gen (Array i Int)
+arbitraryIntArray bnds =
+   fmap (listArray bnds) $ QC.vectorOf (rangeSize bnds) $ QC.choose (-10,10)
+
+genMatrix :: QC.Gen (Array (Int,Int) Int)
+genMatrix = do
+   m <- QC.choose (0,5)
+   n <- QC.choose (0,m)
+   arbitraryIntArray ((1,1),(m,n))
+
+genForward :: QC.Gen (Array (Int,Int) Int, Array Int Int)
+genForward = do
+   a <- genMatrix `QC.suchThat` fullRank
+   let ((_m0,n0), (_m1,n1)) = bounds a
+   x <- arbitraryIntArray (n0,n1)
+   return (a,x)
+
+genInverse :: QC.Gen (Array (Int,Int) Int, Array Int Int)
+genInverse = do
+   a <- genMatrix `QC.suchThat` fullRank
+   let ((m0,_n0), (m1,_n1)) = bounds a
+   b <- arbitraryIntArray (m0,m1)
+   return (a,b)
+
+
+arbitraryIntSparse ::
+   (Ix i, Ix j) => ((i, j), (i, j)) -> QC.Gen (Sparse.Matrix i j Int)
+arbitraryIntSparse bnds =
+   fmap
+      (Sparse.fromMap bnds . fmap fst .
+       Map.filter snd . Map.fromList . zip (range bnds)) $
+      QC.vectorOf (rangeSize bnds) $
+      liftA2 (,) (QC.choose (-10,10)) QC.arbitrary
+
+genSparse :: QC.Gen (Sparse.Matrix Int Int Int)
+genSparse = do
+   m <- QC.choose (0,5)
+   n <- QC.choose (0,m)
+   arbitraryIntSparse ((1,1),(m,n))
+
+genSparseInverse :: QC.Gen (Sparse.Matrix Int Int Int, Array Int Int)
+genSparseInverse = do
+   a <- genSparse `QC.suchThat` (fullRank . Sparse.toDense)
+   let ((m0,_n0), (m1,_n1)) = Sparse.bounds a
+   b <- arbitraryIntArray (m0,m1)
+   return (a,b)
+
+genSquare :: QC.Gen (Array (Int,Int) Int)
+genSquare = do
+   m <- QC.choose (0,5)
+   arbitraryIntArray ((1,1),(m,m))
+
+genSparseSquare :: QC.Gen (Sparse.Matrix Int Int Int)
+genSparseSquare = do
+   m <- QC.choose (0,5)
+   arbitraryIntSparse ((1,1),(m,m))
+
+
+
+approx :: (Fractional a, Ord a) => a -> a -> Bool
+approx x y  =  abs (x-y) <= 1e-5 * max 1 (abs x + abs y)
+
+approxAbsVector :: (Fractional a, Ord a, Ix i) => Array i a -> Array i a -> Bool
+approxAbsVector x y = (Fold.foldl max 0 $ fmap abs $ Vector.sub x y) < 1e-5
+
+
+
+solveHouseholder :: QC.Property
+solveHouseholder =
+   QC.forAll genForward $ \(a,x) ->
+      let b = mv_mult a x
+      in  x ==
+          fmap round (Householder.leastSquares (doubleArray a) (doubleArray b))
+
+solveGivens :: QC.Property
+solveGivens =
+   QC.forAll genForward $ \(a,x) ->
+      let b = mv_mult a x
+      in  x ==
+          fmap round
+            (Givens.leastSquares
+               (Sparse.fromDense $ doubleArray a) (doubleArray b))
+
+
+leastSquares :: QC.Property
+leastSquares =
+   QC.forAll genInverse $ \(a,b) ->
+      Householder.leastSquares (doubleArray a) (doubleArray b)
+      `approxAbsVector`
+      Givens.leastSquares (Sparse.fromDense $ doubleArray a) (doubleArray b)
+
+leastSquaresSparse :: QC.Property
+leastSquaresSparse =
+   QC.forAll genSparseInverse $ \(a,b) ->
+      Householder.leastSquares (doubleArray $ Sparse.toDense a) (doubleArray b)
+      `approxAbsVector`
+      Givens.leastSquares (fmap fromIntegral a) (doubleArray b)
+
+
+
+gramianHouseholder :: QC.Property
+gramianHouseholder =
+   QC.forAll (fmap doubleArray genMatrix) $ \a ->
+      gramian a `approx` (Householder.det a ^! 2)
+
+gramianGivens :: QC.Property
+gramianGivens =
+   QC.forAll (fmap fromIntegral <$> genSparse) $ \a ->
+      gramian (Sparse.toDense a)  `approx`  (Givens.det a ^! 2)
+
+detHouseholder :: QC.Property
+detHouseholder =
+   QC.forAll (fmap doubleArray genSquare) $ \a ->
+      LU.det a `approx` Householder.det a
+
+detGivens :: QC.Property
+detGivens =
+   QC.forAll (fmap fromIntegral <$> genSparseSquare) $ \a ->
+      LU.det (Sparse.toDense a)  `approx`  Givens.det a
+
+
+longCheck :: QC.Property -> IO ()
+longCheck =
+   QC.quickCheckWith (QC.stdArgs {QC.maxSuccess=10000})
+
+tests :: [(String, IO ())]
+tests =
+   ("solveHouseholder", longCheck solveHouseholder) :
+   ("solveGivens", longCheck solveGivens) :
+   ("leastSquares", longCheck leastSquares) :
+   ("leastSquaresSparse", longCheck leastSquaresSparse) :
+   ("gramianHouseholder", longCheck gramianHouseholder) :
+   ("gramianGivens", longCheck gramianGivens) :
+   ("detHouseholder", longCheck detHouseholder) :
+   ("detGivens", longCheck detGivens) :
+   []