diff --git a/repa-linear-algebra.cabal b/repa-linear-algebra.cabal
--- a/repa-linear-algebra.cabal
+++ b/repa-linear-algebra.cabal
@@ -1,5 +1,5 @@
 name:                repa-linear-algebra
-version:             0.2.0.0
+version:             0.3.0.0
 synopsis:            HMatrix operations for Repa.
 description:      HMatrix Vector and Matrix conversions to and from REPA Array F DIM1/2 (Complex) Double, together with a port of linear algebraic functions.
 license:             BSD3
@@ -17,20 +17,20 @@
    type: git
    location: https://github.com/marcinmrotek/repa-linear-algebra.git
 
-Flag Devel
-   Description: Development mode (-Werror).
-   Default:        False
-   Manual:        True
+flag devel
+   description: Development mode (-Werror).
+   default:        False
+   manual:        True
 
 library
   exposed-modules:     Numeric.LinearAlgebra.Repa
                        Numeric.LinearAlgebra.Repa.Conversion
   build-depends:       base    >= 4.8  && <4.9
-                     , hmatrix >= 0.16 && < 0.17
+                     , hmatrix >= 0.17 && < 0.18
                      , repa    >= 3.4  && < 3.5
                      , vector  >= 0.10 && < 0.12
   hs-source-dirs:      src
   ghc-options:       -Wall
-  if flag(Devel)
+  if flag(devel)
     ghc-options:     -Werror
   default-language:    Haskell2010
diff --git a/src/Numeric/LinearAlgebra/Repa.hs b/src/Numeric/LinearAlgebra/Repa.hs
--- a/src/Numeric/LinearAlgebra/Repa.hs
+++ b/src/Numeric/LinearAlgebra/Repa.hs
@@ -1,15 +1,39 @@
+{-|
+Module      : Numeric.LinearAlgebra.Repa
+License     : BSD3
+Maintainer  : marcin.jan.mrotek@gmail.com
+Stability   : experimental
+
+HMatrix linear algebra functions wrapped to accept Repa arrays.
+
+* Unqualified functions use raw 'F' arrays.
+* "-S" functions precompute 'D' arrays sequentially.
+* "-SIO" functions precompute 'D' arrays sequentially in the 'IO' monad.
+* "-P" functions precompute 'D' arrays in parralel in any monad.
+* "-PIO" functions precompute 'D' arrays in parralel in the 'IO' monad.
+
+-}
+
 {-# LANGUAGE FlexibleContexts #-}
 
 module Numeric.LinearAlgebra.Repa
-  ( Numeric
+  (
+  -- * Typeclasses
+    Numeric
   , Field
   , Product
+  , LSDiv
+  -- * Shape-polymorphic conversion
+  , HShape(..)
+  -- * Data types
   , Vector
   , H.Matrix
   , RandDist(..)
   , Seed
-  , HShape(..)
-  , LSDiv
+  -- * Special types 
+  , H.Herm
+  , H.LU
+  , H.LDL
   -- * Dot product
   , dot
   , dotS
@@ -153,10 +177,6 @@
   , null1P
   , null1PIO
   , null1sym
-  , null1symS
-  , null1symSIO
-  , null1symP
-  , null1symPIO
   -- * SVD
   , svd
   , svdS
@@ -195,35 +215,13 @@
   , eigP
   , eigPIO
   , eigSH
-  , eigSH_S
-  , eigSH_SIO
-  , eigSH_P
-  , eigSH_PIO
-  , eigSH'
-  , eigSH'S
-  , eigSH'SIO
-  , eigSH'P
-  , eigSH'PIO
   , eigenvalues
   , eigenvaluesS
   , eigenvaluesSIO
   , eigenvaluesP
   , eigenvaluesPIO
   , eigenvaluesSH
-  , eigenvaluesSH_S
-  , eigenvaluesSH_SIO
-  , eigenvaluesSH_P
-  , eigenvaluesSH_PIO
-  , eigenvaluesSH'
-  , eigenvaluesSH'S
-  , eigenvaluesSH'SIO
-  , eigenvaluesSH'P
-  , eigenvaluesSH'PIO
-  , geigSH'
-  , geigSH'S
-  , geigSH'SIO
-  , geigSH'P
-  , geigSH'PIO
+  , geigSH
   -- * QR
   , qr
   , qrS
@@ -243,15 +241,7 @@
   , qrgr
   -- * Cholesky
   , chol
-  , cholS
-  , cholSIO
-  , cholP
-  , cholPIO
-  , chol'
-  , chol'S
-  , chol'SIO
-  , chol'P
-  , chol'PIO
+  , mbChol
   -- * Hessenberg
   , hess
   , hessS
@@ -275,6 +265,14 @@
   , luPackedSIO
   , luPackedP
   , luPackedPIO
+  , luFact
+  -- * Symmetric indefinite
+  , ldlSolve
+  , ldlSolveS
+  , ldlSolveSIO
+  , ldlSolveP
+  , ldlSolvePIO
+  , ldlPacked
   -- * Matrix functions
   , expm
   , expmS
@@ -325,6 +323,22 @@
   -- *Misc
   , meanCov
   , rowOuters
+  , sym
+  , symS
+  , symSIO
+  , symP
+  , symPIO
+  , mTm
+  , mTmS
+  , mTmSIO
+  , mTmP
+  , mTmPIO
+  , trustSym
+  , trustSymS
+  , trustSymSIO
+  , trustSymP
+  , trustSymPIO
+  , unSym
   ) where
 
 import Numeric.LinearAlgebra.Repa.Conversion
@@ -417,6 +431,7 @@
 outerP :: (Product t, Numeric t, Monad m) => Array D DIM1 t -> Array D DIM1 t -> m (Array F DIM2 t)
 -- |Outer product of two vectors. Arguments computed in parallel.
 outerP v u = hm2repa <$> (H.outer <$> repa2hvP v <*> repa2hvP u)
+
 outerPIO :: (Product t, Numeric t) => Array D DIM1 t -> Array D DIM1 t -> IO (Array F DIM2 t)
 -- |Outer product of two vectors. Arguments computed in parallel inside the IO monad.
 outerPIO v u = hm2repa <$> (H.outer <$> repa2hvPIO v <*> repa2hvPIO u)
@@ -583,21 +598,21 @@
 linearSolveSVD_PIO m n = hm2repa <$> (H.linearSolveLS <$> repa2hmPIO m <*> repa2hmPIO n)
 
 
-luSolve :: (Field t, Numeric t) => PackedLU t -> Array F DIM2 t -> Array F DIM2 t
+luSolve :: (Field t, Numeric t) => H.LU t -> Array F DIM2 t -> Array F DIM2 t
 -- ^Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
-luSolve (PackedLU lu' l) m = hm2repa $ H.luSolve (lu', l) (repa2hm m)
+luSolve lu' m = hm2repa $ H.luSolve lu' (repa2hm m)
 
-luSolveS :: (Field t, Numeric t) => PackedLU t -> Array D DIM2 t -> Array F DIM2 t
-luSolveS (PackedLU lu' l) m = hm2repa $ H.luSolve (lu', l) (repa2hmS m)
+luSolveS :: (Field t, Numeric t) => H.LU t -> Array D DIM2 t -> Array F DIM2 t
+luSolveS lu' m = hm2repa $ H.luSolve lu' (repa2hmS m)
 
-luSolveSIO :: (Field t, Numeric t) => PackedLU t -> Array D DIM2 t -> IO (Array F DIM2 t)
-luSolveSIO (PackedLU lu' l) m = hm2repa . H.luSolve (lu', l) <$> repa2hmSIO m
+luSolveSIO :: (Field t, Numeric t) => H.LU t -> Array D DIM2 t -> IO (Array F DIM2 t)
+luSolveSIO lu' m = hm2repa . H.luSolve lu' <$> repa2hmSIO m
 
-luSolveP :: (Field t, Numeric t, Monad m) => PackedLU t -> Array D DIM2 t -> m (Array F DIM2 t)
-luSolveP (PackedLU lu' l) m = hm2repa . H.luSolve (lu', l) <$> repa2hmP m
+luSolveP :: (Field t, Numeric t, Monad m) => H.LU t -> Array D DIM2 t -> m (Array F DIM2 t)
+luSolveP lu' m = hm2repa . H.luSolve lu' <$> repa2hmP m
 
-luSolvePIO :: (Field t, Numeric t) => PackedLU t -> Array D DIM2 t -> IO (Array F DIM2 t)
-luSolvePIO (PackedLU lu' l) m = hm2repa . H.luSolve (lu', l) <$> repa2hmPIO m
+luSolvePIO :: (Field t, Numeric t) => H.LU t -> Array D DIM2 t -> IO (Array F DIM2 t)
+luSolvePIO lu' m = hm2repa . H.luSolve lu' <$> repa2hmPIO m
 
 
 cholSolve :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
@@ -821,21 +836,9 @@
 null1PIO :: Array D DIM2 Double -> IO (Array F DIM1 Double)
 null1PIO = fmap (hv2repa . H.null1) . repa2hmPIO
 
-null1sym :: Array F DIM2 Double -> Array F DIM1 Double
+null1sym :: H.Herm Double -> Array F DIM1 Double
 -- ^Solution of an overconstrained homogenous symmetric linear system.
-null1sym = hv2repa . H.null1sym . repa2hm
-
-null1symS :: Array D DIM2 Double -> Array F DIM1 Double
-null1symS = hv2repa . H.null1sym . repa2hmS
-
-null1symSIO :: Array D DIM2 Double -> IO (Array F DIM1 Double)
-null1symSIO = fmap (hv2repa . H.null1sym) . repa2hmSIO
-
-null1symP :: Monad m => Array D DIM2 Double -> m (Array F DIM1 Double)
-null1symP = fmap (hv2repa . H.null1sym) . repa2hmP
-
-null1symPIO :: Array D DIM2 Double -> IO (Array F DIM1 Double)
-null1symPIO = fmap (hv2repa . H.null1sym) . repa2hmPIO
+null1sym = hv2repa . H.null1sym
 
 -- SVD
 
@@ -989,49 +992,9 @@
   (s,v) <- H.eig <$> repa2hmPIO m
   return (hv2repa s, hm2repa v)
 
-eigSH :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
+eigSH :: (Field t, Numeric t) => H.Herm t -> (Array F DIM1 Double, Array F DIM2 t)
 -- ^Eigenvalues and eigenvectors (as columns) of a complex hermitian or a real symmetric matrix, in descending order.
-eigSH m = let (s,v) = H.eigSH $ repa2hm m in (hv2repa s, hm2repa v)
-
-eigSH_S :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
-eigSH_S m = let (s,v) = H.eigSH $ repa2hmS m in (hv2repa s, hm2repa v)
-
-eigSH_SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
-eigSH_SIO m = do
-  (s,v) <- H.eigSH <$> repa2hmSIO m
-  return (hv2repa s, hm2repa v)
-
-eigSH_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 Double, Array F DIM2 t)
-eigSH_P m = do
-  (s,v) <- H.eigSH <$> repa2hmP m
-  return (hv2repa s, hm2repa v)
-
-eigSH_PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
-eigSH_PIO m = do
-  (s,v) <- H.eigSH <$> repa2hmPIO m
-  return (hv2repa s, hm2repa v)
-
-eigSH' :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
--- ^Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
-eigSH' m = let (s,v) = H.eigSH' $ repa2hm m in (hv2repa s, hm2repa v)
-
-eigSH'S :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
-eigSH'S m = let (s,v) = H.eigSH' $ repa2hmS m in (hv2repa s, hm2repa v)
-
-eigSH'SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
-eigSH'SIO m = do
-  (s,v) <- H.eigSH' <$> repa2hmSIO m
-  return (hv2repa s, hm2repa v)
-
-eigSH'P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 Double, Array F DIM2 t)
-eigSH'P m = do
-  (s,v) <- H.eigSH' <$> repa2hmP m
-  return (hv2repa s, hm2repa v)
-
-eigSH'PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
-eigSH'PIO m = do
-  (s,v) <- H.eigSH' <$> repa2hmPIO m
-  return (hv2repa s, hm2repa v)
+eigSH h = let (s,v) = H.eigSH h in (hv2repa s, hm2repa v)
 
 eigenvalues :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM1 (Complex Double)
 -- ^Eigenvalues (not ordered) of a general square matrix.
@@ -1049,59 +1012,13 @@
 eigenvaluesPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 (Complex Double))
 eigenvaluesPIO = fmap (hv2repa . H.eigenvalues) . repa2hmPIO
 
-eigenvaluesSH :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM1 Double
+eigenvaluesSH :: (Field t, Numeric t) => H.Herm t -> Array F DIM1 Double
 -- ^Eigenvalues (in descending order) of a complex hermitian or real symmetric matrix.
-eigenvaluesSH = hv2repa . H.eigenvaluesSH . repa2hm
-
-eigenvaluesSH_S :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM1 Double
-eigenvaluesSH_S = hv2repa . H.eigenvaluesSH . repa2hmS
-
-eigenvaluesSH_SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double)
-eigenvaluesSH_SIO = fmap (hv2repa . H.eigenvaluesSH) . repa2hmSIO
-
-eigenvaluesSH_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 Double)
-eigenvaluesSH_P = fmap (hv2repa . H.eigenvaluesSH) . repa2hmP
-
-eigenvaluesSH_PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double)
-eigenvaluesSH_PIO = fmap (hv2repa . H.eigenvaluesSH) . repa2hmPIO
-
-eigenvaluesSH' :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM1 Double
--- ^Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
-eigenvaluesSH' = hv2repa . H.eigenvaluesSH' . repa2hm
-
-eigenvaluesSH'S :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM1 Double
-eigenvaluesSH'S = hv2repa . H.eigenvaluesSH' . repa2hmS
-
-eigenvaluesSH'SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double)
-eigenvaluesSH'SIO = fmap (hv2repa . H.eigenvaluesSH') . repa2hmSIO
-
-eigenvaluesSH'P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 Double)
-eigenvaluesSH'P = fmap (hv2repa . H.eigenvaluesSH') . repa2hmP
-
-eigenvaluesSH'PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double)
-eigenvaluesSH'PIO = fmap (hv2repa . H.eigenvaluesSH') . repa2hmPIO
-
-geigSH' :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
--- ^Generalized symmetric positive definite eigensystem Av = IBv, for A and B symmetric, B positive definite (conditions not checked).
-geigSH' a b = let (s,v) = H.geigSH' (repa2hm a) (repa2hm b) in (hv2repa s, hm2repa v)
-
-geigSH'S :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
-geigSH'S a b = let (s,v) = H.geigSH' (repa2hmS a) (repa2hmS b) in (hv2repa s, hm2repa v)
-
-geigSH'SIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
-geigSH'SIO a b = do
-  (s,v) <- H.geigSH' <$> repa2hmSIO a <*> repa2hmSIO b
-  return (hv2repa s, hm2repa v)
-
-geigSH'P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> Array D DIM2 t -> m (Array F DIM1 Double, Array F DIM2 t)
-geigSH'P a b = do
-  (s,v) <- H.geigSH' <$> repa2hmP a <*> repa2hmP b
-  return (hv2repa s, hm2repa v)
+eigenvaluesSH = hv2repa . H.eigenvaluesSH 
 
-geigSH'PIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
-geigSH'PIO a b = do
-  (s,v) <- H.geigSH' <$> repa2hmPIO a <*> repa2hmPIO b
-  return (hv2repa s, hm2repa v)
+geigSH :: (Field t, Numeric t) => H.Herm t -> H.Herm t -> (Array F DIM1 Double, Array F DIM2 t)
+-- ^Generalized symmetric positive definite eigensystem Av = IBv, for A and B symmetric, B positive definite.
+geigSH a b = let (s,v) = H.geigSH a b in (hv2repa s, hm2repa v)
 
 -- QR
 
@@ -1149,64 +1066,34 @@
   (r,q) <- H.rq <$> repa2hmPIO m
   return (hm2repa r, hm2repa q)
 
-qrRaw :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM1 t)
-qrRaw m = let (n,v) = H.qrRaw $ repa2hm m in (hm2repa n, hv2repa v)
+qrRaw :: (Field t, Numeric t) => Array F DIM2 t -> H.QR t
+qrRaw m = H.qrRaw $ repa2hm m
 
-qrRawS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM1 t)
-qrRawS m = let (n,v) = H.qrRaw $ repa2hmS m in (hm2repa n, hv2repa v)
+qrRawS :: (Field t, Numeric t) => Array D DIM2 t -> H.QR t
+qrRawS m = H.qrRaw $ repa2hmS m
 
-qrRawSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 t)
-qrRawSIO m = do
-  (n,v) <- H.qrRaw <$> repa2hmSIO m
-  return (hm2repa n, hv2repa v)
+qrRawSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (H.QR t)
+qrRawSIO m = H.qrRaw <$> repa2hmSIO m
 
-qrRawP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM1 t)
-qrRawP m = do
-  (n,v) <- H.qrRaw <$> repa2hmP m
-  return (hm2repa n, hv2repa v)
+qrRawP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (H.QR t)
+qrRawP m = H.qrRaw <$> repa2hmP m
 
-qrRawPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 t)
-qrRawPIO m = do
-  (n,v) <- H.qrRaw <$> repa2hmPIO m
-  return (hm2repa n, hv2repa v)
+qrRawPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (H.QR t)
+qrRawPIO m = H.qrRaw <$> repa2hmPIO m
 
-qrgr :: (Field t, Numeric t) => Int -> (Array F DIM2 t, Array F DIM1 t) -> Array F DIM2 t
+qrgr :: (Field t, Numeric t) => Int -> H.QR t -> Array F DIM2 t
 -- ^Generate a matrix with k othogonal columns from the output of 'qrRaw'.
-qrgr k (m,v) = hm2repa $ H.qrgr k (repa2hm m, repa2hv v)
+qrgr k qr' = hm2repa $ H.qrgr k qr'
 
 -- Cholesky
 
-chol :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
+chol :: Field t => H.Herm t -> Array F DIM2 t
 -- ^Cholesky factorization of a positive definite hermitian or symmetric matrix. c = chol m ==> m == c' * c where c is upper triangular.
-chol = hm2repa . H.chol . repa2hm
-
-cholS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
-cholS = hm2repa . H.chol . repa2hmS
-
-cholSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
-cholSIO = fmap (hm2repa . H.chol) . repa2hmSIO
-
-cholP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
-cholP = fmap (hm2repa . H.chol) . repa2hmP
-
-cholPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
-cholPIO = fmap (hm2repa . H.chol) . repa2hmPIO
-
-chol' :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
--- ^Similar to 'chol' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
-chol' = hm2repa . H.cholSH . repa2hm
-
-chol'S :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
-chol'S = hm2repa . H.cholSH . repa2hmS
-
-chol'SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
-chol'SIO = fmap (hm2repa . H.cholSH) . repa2hmSIO
-
-chol'P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
-chol'P = fmap (hm2repa . H.cholSH) . repa2hmP
+chol = hm2repa . H.chol 
 
-chol'PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
-chol'PIO = fmap (hm2repa . H.cholSH) . repa2hmPIO
+mbChol :: Field t => H.Herm t -> Maybe (Array F DIM2 t)
+-- ^Similar to chol, but instead of an error (e.g., caused by a matrix not positive definite) it returns Nothing.
+mbChol h = hm2repa <$> H.mbChol h
 
 -- Hessenberg
 
@@ -1279,30 +1166,52 @@
   (l,u,p,s) <- H.lu <$> repa2hmPIO m
   return (hm2repa l, hm2repa u, hm2repa p, s)
 
-data PackedLU t = PackedLU (H.Matrix t) [Int]
-
-luPacked :: (Field t, Numeric t) => Array F DIM2 t -> PackedLU t
+luPacked :: (Field t, Numeric t) => Array F DIM2 t -> H.LU t
 -- ^Obtains the LU decomposition in a packed data structure suitable for 'luSolve'.
-luPacked m = let (lu', is) = H.luPacked $ repa2hm m in PackedLU lu' is
+luPacked m =  H.luPacked $ repa2hm m 
 
-luPackedS :: (Field t, Numeric t) => Array D DIM2 t -> PackedLU t
-luPackedS m = let (lu', is) = H.luPacked $ repa2hmS m in PackedLU lu' is
+luPackedS :: (Field t, Numeric t) => Array D DIM2 t -> H.LU t
+luPackedS m = H.luPacked $ repa2hmS m
 
-luPackedSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (PackedLU t)
-luPackedSIO m = do
-  (lu', is) <- H.luPacked <$> repa2hmSIO m
-  return $ PackedLU lu' is
+luPackedSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (H.LU t)
+luPackedSIO m = H.luPacked <$> repa2hmSIO m
 
-luPackedP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (PackedLU t)
-luPackedP m = do
-  (lu', is) <- H.luPacked <$> repa2hmP m
-  return $ PackedLU lu' is
+luPackedP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (H.LU t)
+luPackedP m = H.luPacked <$> repa2hmP m
 
-luPackedPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (PackedLU t)
-luPackedPIO m = do
-  (lu', is) <- H.luPacked <$> repa2hmPIO m
-  return $ PackedLU lu' is
+luPackedPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (H.LU t)
+luPackedPIO m = H.luPacked <$> repa2hmPIO m
 
+luFact :: Numeric t => H.LU t -> (Array F DIM2 t, Array F DIM2 t, Array F DIM2 t, t)
+-- ^Compute the explicit LU decomposition from the compact one obtained by luPacked.
+luFact lu' = let (l,u,p,s) = H.luFact lu' in (hm2repa l, hm2repa u, hm2repa p, s)
+
+-- Symmetric indefinite
+
+ldlSolve :: Field t => H.LDL t -> Array F DIM2 t -> Array F DIM2 t
+{- ^
+Solution of a linear system (for several right hand sides) from a precomputed LDL factorization obtained by 'ldlPacked'.
+Note: this can be slower than the general solver based on the LU decomposition.
+-}
+
+ldlSolve l = hm2repa . H.ldlSolve l . repa2hm
+
+ldlSolveS :: Field t => H.LDL t -> Array D DIM2 t -> Array F DIM2 t
+ldlSolveS l = hm2repa . H.ldlSolve l . repa2hmS
+
+ldlSolveSIO :: Field t => H.LDL t -> Array D DIM2 t -> IO (Array F DIM2 t)
+ldlSolveSIO l m = hm2repa . H.ldlSolve l <$> repa2hmSIO m
+
+ldlSolveP :: (Field t, Monad m) => H.LDL t -> Array D DIM2 t -> m (Array F DIM2 t)
+ldlSolveP l m = hm2repa . H.ldlSolve l <$> repa2hmP m
+
+ldlSolvePIO :: Field t => H.LDL t -> Array D DIM2 t -> IO (Array F DIM2 t)
+ldlSolvePIO l m = hm2repa . H.ldlSolve l <$> repa2hmPIO m
+
+ldlPacked :: Field t => H.Herm t -> H.LDL t
+-- ^Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'.
+ldlPacked = H.ldlPacked
+
 -- Matrix functions
 
 expm :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
@@ -1464,3 +1373,55 @@
 rowOuters :: Array F DIM2 Double -> Array F DIM2 Double -> Array F DIM2 Double
 -- ^Outer product of the rows of the matrices.
 rowOuters m n = hm2repa $ H.rowOuters (repa2hm m) (repa2hm n)
+
+sym :: Field t => Array F DIM2 t -> H.Herm t
+-- ^Compute the complex Hermitian or real symmetric part of a square matrix ((x + tr x)/2).
+sym = H.sym . repa2hm
+
+symS :: Field t => Array D DIM2 t -> H.Herm t
+symS = H.sym . repa2hmS
+
+symSIO :: Field t => Array D DIM2 t -> IO (H.Herm t)
+symSIO m = H.sym <$> repa2hmSIO m
+
+symP :: (Field t, Monad m) => Array D DIM2 t -> m (H.Herm t)
+symP m = H.sym <$> repa2hmP m
+
+symPIO :: Field t => Array D DIM2 t -> IO (H.Herm t)
+symPIO m = H.sym <$> repa2hmPIO m
+
+mTm :: Field t => Array F DIM2 t -> H.Herm t
+-- ^Compute the contraction tr x <> x of a general matrix.
+mTm = H.mTm . repa2hm
+
+mTmS :: Field t => Array D DIM2 t -> H.Herm t
+mTmS = H.mTm . repa2hmS
+
+mTmSIO :: Field t => Array D DIM2 t -> IO (H.Herm t)
+mTmSIO m = H.mTm <$> repa2hmSIO m
+
+mTmP :: (Field t, Monad m) => Array D DIM2 t -> m (H.Herm t)
+mTmP m = H.mTm <$> repa2hmP m
+
+mTmPIO :: Field t => Array D DIM2 t -> IO (H.Herm t)
+mTmPIO m = H.mTm <$> repa2hmPIO m
+
+trustSym :: Field t => Array F DIM2 t -> H.Herm t
+-- ^At your own risk, declare that a matrix is complex Hermitian or real symmetric for usage in 'chol', 'eigSH', etc. Only a triangular part of the matrix will be used.
+trustSym = H.trustSym . repa2hm
+
+trustSymS :: Field t => Array D DIM2 t -> H.Herm t
+trustSymS = H.trustSym . repa2hmS
+
+trustSymSIO :: Field t => Array D DIM2 t -> IO (H.Herm t)
+trustSymSIO m = H.trustSym <$> repa2hmSIO m 
+
+trustSymP :: (Field t, Monad m) => Array D DIM2 t -> m (H.Herm t)
+trustSymP m = H.trustSym <$> repa2hmP m 
+
+trustSymPIO :: Field t => Array D DIM2 t -> IO (H.Herm t)
+trustSymPIO m = H.trustSym <$> repa2hmPIO m 
+
+unSym :: Numeric t => H.Herm t -> Array F DIM2 t
+-- ^Extract the general matrix from a Herm structure, forgetting its symmetric or Hermitian property.
+unSym = hm2repa . H.unSym
diff --git a/src/Numeric/LinearAlgebra/Repa/Conversion.hs b/src/Numeric/LinearAlgebra/Repa/Conversion.hs
--- a/src/Numeric/LinearAlgebra/Repa/Conversion.hs
+++ b/src/Numeric/LinearAlgebra/Repa/Conversion.hs
@@ -1,3 +1,19 @@
+{-|
+Module      : Numeric.LinearAlgebra.Repa.Conversion
+License     : BSD3
+Maintainer  : marcin.jan.mrotek@gmail.com
+Stability   : experimental
+
+Repa - HMatrix conversion functions.
+
+* Unqualified functions use raw 'F' arrays.
+* "-S" functions precompute 'D' arrays sequentially.
+* "-SIO" functions precompute 'D' arrays sequentially in the 'IO' monad.
+* "-P" functions precompute 'D' arrays in parralel in any monad.
+* "-PIO" functions precompute 'D' arrays in parralel in the 'IO' monad.
+
+-}
+
 {-# LANGUAGE
     FlexibleContexts
   , FlexibleInstances
@@ -103,6 +119,7 @@
   :: ( Storable t
      , Container V.Vector t
      , Element t
+     , Num t
      )
   => H.Matrix t -> Array F DIM2 t
 -- ^O(1). Convert a HMatrix Matrix to a Repa Array.