diff --git a/Data/Matrix.hs b/Data/Matrix.hs
--- a/Data/Matrix.hs
+++ b/Data/Matrix.hs
@@ -1,6 +1,8 @@
+
 -- | Matrix datatype and operations.
 --
 --   Every provided example has been tested.
+--   Run @cabal test@ for further tests.
 module Data.Matrix (
     -- * Matrix type
     Matrix , prettyMatrix
@@ -19,18 +21,21 @@
   , getElem , (!) , safeGet
   , getRow  , getCol
   , getDiag
+  , getMatrixAsVector
     -- * Manipulating matrices
   , setElem
-  , transpose , extendTo
+  , transpose , setSize , extendTo
   , mapRow , mapCol
     -- * Submatrices
     -- ** Splitting blocks
   , submatrix
   , minorMatrix
   , splitBlocks
-    -- ** Joining blocks
+   -- ** Joining blocks
   , (<|>) , (<->)
   , joinBlocks
+    -- * Matrix operations
+  , elementwise
     -- * Matrix multiplication
     -- ** About matrix multiplication
     -- $mult
@@ -46,8 +51,8 @@
   , switchRows
   , switchCols
     -- * Decompositions
-  , luDecomp
-  , luDecomp'
+  , luDecomp , luDecompUnsafe
+  , luDecomp', luDecompUnsafe'
   , cholDecomp
     -- * Properties
   , trace , diagProd
@@ -64,7 +69,7 @@
 import Data.Traversable
 -- Data
 import           Control.Monad.Primitive (PrimMonad, PrimState)
-import           Data.List               (maximumBy)
+import           Data.List               (maximumBy,foldl1')
 import           Data.Ord                (comparing)
 import qualified Data.Vector             as V
 import qualified Data.Vector.Mutable     as MV
@@ -87,7 +92,7 @@
 data Matrix a = M {
    nrows :: {-# UNPACK #-} !Int -- ^ Number of rows.
  , ncols :: {-# UNPACK #-} !Int -- ^ Number of columns.
- , mvect :: (V.Vector a) -- ^ Content of the matrix as a plain vector.
+ , mvect :: V.Vector a          -- ^ Content of the matrix as a plain vector.
    } deriving Eq
 
 -- | Just a cool way to output the size of a matrix.
@@ -119,8 +124,12 @@
 ---- FUNCTOR INSTANCE
 
 instance Functor Matrix where
+ {-# INLINE fmap #-}
  fmap f (M n m v) = M n m $ V.map f v
 
+-------------------------------------------------------
+-------------------------------------------------------
+
 -- | /O(rows*cols)/. Map a function over a row.
 --   Example:
 --
@@ -218,7 +227,7 @@
          -> [a] -- ^ List of elements
          -> Matrix a
 {-# INLINE fromList #-}
-fromList n m = M n m . V.fromList
+fromList n m = M n m . V.fromListN (n*m)
 
 -- | Create a matrix from an non-empty list of non-empty lists.
 --   /Each list must have the same number of elements/.
@@ -280,7 +289,7 @@
         -> Matrix a -- ^ Matrix
         -> a
 {-# INLINE getElem #-}
-getElem i j (M _ m v) = v V.! encode m (i,j)
+getElem i j (M _ m v) = V.unsafeIndex v $ encode m (i,j)
 
 -- | Short alias for 'getElem'.
 {-# INLINE (!) #-}
@@ -290,12 +299,12 @@
 -- | Safe variant of 'getElem'.
 safeGet :: Int -> Int -> Matrix a -> Maybe a
 safeGet i j a@(M n m _)
- | i > n || j > m = Nothing
+ | i > n || j > m || i < 1 || j < 1 = Nothing
  | otherwise = Just $ getElem i j a
 
--- | /O(cols)/. Get a row of a matrix as a vector.
+-- | /O(1)/. Get a row of a matrix as a vector.
 getRow :: Int -> Matrix a -> V.Vector a
-getRow i (M _ m v) = V.generate m $ \j -> v V.! encode m (i,j+1)
+getRow i (M _ m v) = V.slice (m*(i-1)) m v
 
 -- | /O(rows)/. Get a column of a matrix as a vector.
 getCol :: Int -> Matrix a -> V.Vector a
@@ -307,6 +316,12 @@
  where
   k = min (nrows m) (ncols m)
 
+-- | /O(1)/. Transform a 'Matrix' to a 'V.Vector' of size /rows*cols/.
+--  This is equivalent to get all the rows of the matrix using 'getRow'
+--  and then append them, but far more efficient.
+getMatrixAsVector :: Matrix a -> V.Vector a
+getMatrixAsVector = mvect
+
 -------------------------------------------------------
 -------------------------------------------------------
 ---- MANIPULATING MATRICES
@@ -330,26 +345,36 @@
 transpose :: Matrix a -> Matrix a
 transpose m = matrix (ncols m) (nrows m) $ \(i,j) -> m ! (j,i)
 
--- | Extend a matrix to a given size adding zeroes.
+-- | Extend a matrix to a given size adding a default element.
 --   If the matrix already has the required size, nothing happens.
 --   The matrix is /never/ reduced in size.
 --   Example:
 --
--- >                          ( 1 2 3 0 0 )
--- >              ( 1 2 3 )   ( 4 5 6 0 0 )
--- >              ( 4 5 6 )   ( 7 8 9 0 0 )
--- > extendTo 4 5 ( 7 8 9 ) = ( 0 0 0 0 0 )
-extendTo :: Num a
-         => Int -- ^ Minimal number of rows.
+-- >                            ( 1 2 3 0 0 )
+-- >                ( 1 2 3 )   ( 4 5 6 0 0 )
+-- >                ( 4 5 6 )   ( 7 8 9 0 0 )
+-- > extendTo 0 4 5 ( 7 8 9 ) = ( 0 0 0 0 0 )
+extendTo :: a   -- ^ Element to add when extending.
+         -> Int -- ^ Minimal number of rows.
          -> Int -- ^ Minimal number of columns.
          -> Matrix a -> Matrix a
-extendTo n m a = a''
- where
-  n'  = n - nrows a
-  a'  = if n' <= 0 then a  else a  <-> zero n' (ncols a)
-  m'  = m - ncols a
-  a'' = if m' <= 0 then a' else a' <|> zero (nrows a') m'
+extendTo e n m a = setSize e (max n $ nrows a) (max m $ ncols a) a
 
+-- | Set the size of a matrix to given parameters. Use a default element
+--   for undefined entries if the matrix has been extended.
+setSize :: a   -- ^ Default element.
+        -> Int -- ^ Number of rows.
+        -> Int -- ^ Number of columns.
+        -> Matrix a
+        -> Matrix a
+setSize e n m a@(M r c _) = M n m $ V.generate (n*m) $
+  \k -> let (i,j) = d k
+        in if i <= r && j <= c
+              then a ! (i,j)
+              else e
+    where
+      d = decode m
+
 -------------------------------------------------------
 -------------------------------------------------------
 ---- WORKING WITH BLOCKS
@@ -411,15 +436,12 @@
             -> Matrix a -- ^ Matrix to split.
             -> (Matrix a,Matrix a
                ,Matrix a,Matrix a) -- ^ (TL,TR,BL,BR)
-{-# INLINE splitBlocks #-}
 splitBlocks i j a@(M n m _) = ( submatrix    1  i 1 j a , submatrix    1  i (j+1) m a
                               , submatrix (i+1) n 1 j a , submatrix (i+1) n (j+1) m a )
 
+
 -- | Join blocks of the form detailed in 'splitBlocks'.
-joinBlocks :: (Matrix a,Matrix a
-              ,Matrix a,Matrix a)
-           ->  Matrix a
-{-# INLINE joinBlocks #-}
+joinBlocks :: (Matrix a,Matrix a,Matrix a,Matrix a) -> Matrix a
 joinBlocks (tl,tr,bl,br) = (tl <|> tr)
                                <->
                            (bl <|> br)
@@ -458,6 +480,15 @@
 
 -------------------------------------------------------
 -------------------------------------------------------
+---- MATRIX OPERATIONS
+
+-- | Perform an operation elementwise. The input matrices are assumed
+--   to have the same dimensions, but this is not checked.
+elementwise :: (a -> b -> c) -> (Matrix a -> Matrix b -> Matrix c)
+elementwise f (M n m v) (M _ _ v') = M n m $ V.zipWith f v v'
+
+-------------------------------------------------------
+-------------------------------------------------------
 ---- MATRIX MULTIPLICATION
 
 {- $mult
@@ -542,12 +573,12 @@
    | otherwise =
        let mx = maximum [n,m,n',m']
            n2  = first (>= mx) $ fmap (2^) [(0 :: Int)..]
-           b1 = extendTo n2 n2 a1
-           b2 = extendTo n2 n2 a2
+           b1 = setSize 0 n2 n2 a1
+           b2 = setSize 0 n2 n2 a2
        in  submatrix 1 n 1 m' $ strassen b1 b2
 
 strmixFactor :: Int
-strmixFactor = 100
+strmixFactor = 2 ^ (6 :: Int)
 
 -- | Strassen's mixed algorithm.
 strassenMixed :: Num a => Matrix a -> Matrix a -> Matrix a
@@ -556,8 +587,8 @@
 strassenMixed a@(M r _ _) b
  | r < strmixFactor = multStd_ a b
  | odd r = let r' = r + 1
-               a' = extendTo r' r' a
-               b' = extendTo r' r' b
+               a' = setSize 0 r' r' a
+               b' = setSize 0 r' r' b
            in  submatrix 1 r 1 r $ strassenMixed a' b'
  | otherwise = joinBlocks (c11,c12,c21,c22)
  where
@@ -590,8 +621,8 @@
    | otherwise =
        let mx = maximum [n,m,n',m']
            n2 = if even mx then mx else mx+1
-           b1 = extendTo n2 n2 a1
-           b2 = extendTo n2 n2 a2
+           b1 = setSize 0 n2 n2 a1
+           b2 = setSize 0 n2 n2 a2
        in  submatrix 1 n 1 m' $ strassenMixed b1 b2
 
 -------------------------------------------------------
@@ -719,10 +750,10 @@
 -- >          ( 1 2 0 )     ( 2 0  2 )   (   1 0 0 )   ( 0 0 1 )
 -- >          ( 0 2 1 )     ( 0 2 -1 )   ( 1/2 1 0 )   ( 1 0 0 )
 -- > luDecomp ( 2 0 2 ) = ( ( 0 0  2 ) , (   0 1 1 ) , ( 0 1 0 ) , 1 )
-luDecomp :: (Ord a, Fractional a) => Matrix a -> (Matrix a,Matrix a,Matrix a,a)
+luDecomp :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a,Matrix a,Matrix a,a)
 luDecomp a = recLUDecomp a i i 1 1 n
  where
-  i = (identity $ nrows a)
+  i = identity $ nrows a
   n = min (nrows a) (ncols a)
 
 recLUDecomp ::  (Ord a, Fractional a)
@@ -732,10 +763,11 @@
             ->  a        -- ^ d
             ->  Int      -- ^ Current row
             ->  Int      -- ^ Total rows
-            -> (Matrix a,Matrix a,Matrix a,a)
+            -> Maybe (Matrix a,Matrix a,Matrix a,a)
 recLUDecomp u l p d k n =
-    if k > n then (u,l,p,d)
-    else recLUDecomp u'' l'' p' d' (k+1) n
+    if k > n then Just (u,l,p,d)
+    else if ukk == 0 then Nothing
+                     else recLUDecomp u'' l'' p' d' (k+1) n
  where
   -- Pivot strategy: maximum value in absolute value below the current row.
   i  = maximumBy (\x y -> compare (abs $ u ! (x,k)) (abs $ u ! (y,k))) [ k .. n ]
@@ -758,6 +790,12 @@
     else let x = (u_ ! (j,k)) / ukk
          in  go (combineRows j (-x) k u_) (setElem x (j,k) l_) (j+1)
 
+-- | Unsafe version of 'luDecomp'. It fails when the input matrix is singular.
+luDecompUnsafe :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, a)
+luDecompUnsafe m = case luDecomp m of
+  Just x -> x
+  _ -> error "luDecompUnsafe of singular matrix."
+
 -- | Matrix LU decomposition with /complete pivoting/.
 --   The result for a matrix /M/ is given in the format /(U,L,P,Q,d,e)/ where:
 --
@@ -765,9 +803,9 @@
 --
 --   * /L/ is an /unit/ lower triangular matrix.
 --
---   * /P,Q/ is a permutation matrix.
+--   * /P,Q/ are permutation matrices.
 --
---   * /d,e/ is the determinant of /P,Q/.
+--   * /d,e/ are the determinants of /P/ and /Q/ respectively.
 --
 --   * /PMQ = LU/.
 --
@@ -784,12 +822,18 @@
 -- >           ( 1 0 )     ( 2 1 )   (   1    0 0 )   ( 0 0 1 )
 -- >           ( 0 2 )     ( 0 2 )   (   0    1 0 )   ( 0 1 0 )   ( 1 0 )
 -- > luDecomp' ( 2 1 ) = ( ( 0 0 ) , ( 1/2 -1/4 1 ) , ( 1 0 0 ) , ( 0 1 ) , -1 , 1 )
-luDecomp' :: (Ord a, Fractional a) => Matrix a -> (Matrix a,Matrix a,Matrix a,Matrix a,a,a)
+luDecomp' :: (Ord a, Fractional a) => Matrix a -> Maybe (Matrix a,Matrix a,Matrix a,Matrix a,a,a)
 luDecomp' a = recLUDecomp' a i i (identity $ ncols a) 1 1 1 n
  where
   i = identity $ nrows a
   n = min (nrows a) (ncols a)
 
+-- | Unsafe version of 'luDecomp''. It fails when the input matrix is singular.
+luDecompUnsafe' :: (Ord a, Fractional a) => Matrix a -> (Matrix a, Matrix a, Matrix a, Matrix a, a, a)
+luDecompUnsafe' m = case luDecomp' m of
+  Just x -> x
+  _ -> error "luDecompUnsafe' of singular matrix."
+
 recLUDecomp' ::  (Ord a, Fractional a)
             =>  Matrix a -- ^ U
             ->  Matrix a -- ^ L
@@ -799,11 +843,13 @@
             ->  a        -- ^ e
             ->  Int      -- ^ Current row
             ->  Int      -- ^ Total rows
-            -> (Matrix a,Matrix a,Matrix a,Matrix a,a,a)
+            ->  Maybe (Matrix a,Matrix a,Matrix a,Matrix a,a,a)
 recLUDecomp' u l p q d e k n =
     if k > n || u'' ! (k, k) == 0
-    then (u,l,p,q,d,e)
-    else recLUDecomp' u'' l'' p' q' d' e' (k+1) n
+    then Just (u,l,p,q,d,e)
+    else if ukk == 0
+            then Nothing
+            else recLUDecomp' u'' l'' p' q' d' e' (k+1) n
  where
   -- Pivot strategy: maximum value in absolute value below the current row & col.
   (i, j) = maximumBy (comparing (\(i0, j0) -> abs $ u ! (i0,j0)))
@@ -892,23 +938,27 @@
 -- DETERMINANT
 
 {-# RULES
-"matrix/detOfProduct"
+"matrix/detLaplaceProduct"
     forall a b. detLaplace (a*b) = detLaplace a * detLaplace b
 
-"matrix/detLUOfProduct"
+"matrix/detLUProduct"
     forall a b. detLU (a*b) = detLU a * detLU b
   #-}
 
 -- | Matrix determinant using Laplace expansion.
 --   If the elements of the 'Matrix' are instance of 'Ord' and 'Fractional'
 --   consider to use 'detLU' in order to obtain better performance.
+--   Function 'detLaplace' is /extremely/ slow.
 detLaplace :: Num a => Matrix a -> a
 detLaplace (M 1 1 v) = V.head v
-detLaplace m =
-    sum [ (-1)^(i-1) * m ! (i,1) * detLaplace (minorMatrix i 1 m) | i <- [1 .. nrows m] ]
+detLaplace m = sum1 [ (-1)^(i-1) * m ! (i,1) * detLaplace (minorMatrix i 1 m) | i <- [1 .. nrows m] ]
+  where
+    sum1 = foldl1' (+)
 
 -- | Matrix determinant using LU decomposition.
+--   It works even when the input matrix is singular.
 detLU :: (Ord a, Fractional a) => Matrix a -> a
-detLU m = d * diagProd u
- where
-  (u,_,_,d) = luDecomp m
+detLU m = case luDecomp m of
+  Just (u,_,_,d) -> d * diagProd u
+  Nothing -> 0
+
diff --git a/license b/license
--- a/license
+++ b/license
@@ -1,4 +1,4 @@
-Copyright (c)2013, Daniel Díaz
+Copyright (c)2014, Daniel Díaz
 
 All rights reserved.
 
diff --git a/matrix.cabal b/matrix.cabal
--- a/matrix.cabal
+++ b/matrix.cabal
@@ -1,5 +1,5 @@
 Name: matrix
-Version: 0.2.4.0
+Version: 0.3.0.0
 Author: Daniel Díaz
 Category: Math
 Build-type: Simple
@@ -29,7 +29,7 @@
   location: git://github.com/Daniel-Diaz/matrix.git
 
 Library
-  Build-depends: base ==4.*
+  Build-depends: base == 4.*
                , vector >= 0.10 && < 0.11
                , deepseq >= 1.3.0.0 && < 1.4
                , primitive >= 0.5 && < 0.6
@@ -40,7 +40,17 @@
   type: exitcode-stdio-1.0
   hs-source-dirs: bench
   main-is: mult.hs
-  build-depends: base ==4.*
+  build-depends: base == 4.*
                , matrix
                , criterion
   ghc-options: -O2
+
+Test-Suite matrix-test
+  type: exitcode-stdio-1.0
+  hs-source-dirs: test
+  main-is: Main.hs
+  build-depends: base == 4.*
+               , matrix
+               , tasty
+               , QuickCheck
+               , tasty-quickcheck
diff --git a/test/Main.hs b/test/Main.hs
new file mode 100644
--- /dev/null
+++ b/test/Main.hs
@@ -0,0 +1,102 @@
+
+import Data.Matrix
+import Data.Ratio
+import Control.Applicative
+import Data.Monoid (mconcat)
+
+import Test.Tasty
+import qualified Test.Tasty.QuickCheck as QC
+import Test.QuickCheck
+
+{- matrix package test set
+
+This program uses QuickCheck to check that the matrix
+functions of the matrix package are working properly.
+
+We use the type Rational to have avoid numerical errors
+that may cause the test to fail while the algorithm is
+correct.
+
+-}
+
+-- | Numbers used in tests.
+type R = Rational
+
+newtype I = I { fromI :: Int }
+
+instance Show I where
+  show (I n) = show n
+
+instance Arbitrary I where
+  arbitrary = I <$> choose (1,9)
+
+instance Arbitrary a => Arbitrary (Matrix a) where
+  arbitrary = do
+    I n <- arbitrary
+    I m <- arbitrary
+    genMatrix' n m
+
+genMatrix' :: Arbitrary a => Int -> Int -> Gen (Matrix a)
+genMatrix' n m = fromList n m <$> vector (n*m)
+
+genMatrix :: Int -> Int -> Gen (Matrix R)
+genMatrix = genMatrix'
+
+
+-- | Square matrices
+newtype Sq = Sq { fromSq :: Matrix R }
+
+instance Show Sq where
+  show (Sq m) = show m
+
+instance Arbitrary Sq where
+  arbitrary = do
+    I n <- arbitrary
+    Sq <$> genMatrix n n
+
+main :: IO ()
+main = defaultMain $ testGroup "matrix tests" [
+    QC.testProperty "identity * m = m * identity = m"
+       $ \(Sq m) -> let n = nrows m in identity n * m == m && m * identity n == m
+  , QC.testProperty "getMatrixAsVector m = mconcat [ getRow i m | i <- [1 .. nrows m]]"
+      $ \m -> getMatrixAsVector (m :: Matrix R) == mconcat [ getRow i m | i <- [1 .. nrows m] ]
+  , QC.testProperty "permMatrix n i j * permMatrix n i j = identity n"
+       $ \(I n) -> forAll (choose (1,n))
+       $ \i     -> forAll (choose (1,n))
+       $ \j     -> permMatrix n i j * permMatrix n i j == identity n
+  , QC.testProperty "setElem (getElem i j m) (i,j) m = m"
+       $ \m -> forAll (choose (1,nrows m))
+       $ \i -> forAll (choose (1,ncols m))
+       $ \j -> setElem (getElem i j m) (i,j) m == (m :: Matrix R)
+  , QC.testProperty "transpose (transpose m) = m"
+       $ \m -> transpose (transpose m) == (m :: Matrix R)
+  , QC.testProperty "getRow i m = getCol i (transpose m)"
+       $ \m -> forAll (choose (1,nrows m))
+       $ \i -> getRow i (m :: Matrix R) == getCol i (transpose m)
+  , QC.testProperty "joinBlocks (splitBlocks i j m) = m"
+       $ \m -> forAll (choose (1,nrows m))
+       $ \i -> forAll (choose (1,ncols m))
+       $ \j -> joinBlocks (splitBlocks i j m) == (m :: Matrix R)
+  , QC.testProperty "(+) = elementwise (+)"
+       $ \m1 -> forAll (genMatrix (nrows m1) (ncols m1))
+       $ \m2 -> m1 + m2 == elementwise (+) m1 m2
+  , QC.testProperty "if (u,l,p,d) = luDecomp m then (p*m = l*u) && (detLaplace p = d)"
+       $ \(Sq m) -> (detLaplace m /= 0) ==>
+             (let (u,l,p,d) = luDecompUnsafe m in p*m == l*u && detLaplace p == d)
+  , QC.testProperty "detLaplace m = detLU m"
+       $ \(Sq m) -> detLaplace m == detLU m
+  , QC.testProperty "if (u,l,p,q,d,e) = luDecomp' m then (p*m*q = l*u) && (detLU p = d) && (detLU q = e)"
+       $ \(Sq m) -> (detLU m /= 0) ==>
+             (let (u,l,p,q,d,e) = luDecompUnsafe' m in p*m*q == l*u && detLU p == d && detLU q == e)
+  , QC.testProperty "detLU (scaleRow k i m) = k * detLU m"
+       $ \(Sq m) k -> forAll (choose (1,nrows m))
+       $ \i -> detLU (scaleRow k i m) == k * detLU m
+  , QC.testProperty "let n = nrows m in detLU (switchRows i j m) = detLU (permMatrix n i j) * detLU m"
+       $ \(Sq m) -> let n = nrows m in forAll (choose (1,n))
+       $ \i      -> forAll (choose (1,n))
+       $ \j      -> detLU (switchRows i j m) == detLU (permMatrix n i j) * detLU m
+  , QC.testProperty "switchCols i j = transpose . switchRows i j . transpose"
+       $ \m -> forAll (choose (1,ncols m))
+       $ \i -> forAll (choose (1,ncols m))
+       $ \j -> switchCols i j (m :: Matrix R) == (transpose $ switchRows i j $ transpose m)
+    ]