packages feed

repa-linear-algebra 0.2.0.0 → 0.3.0.0

raw patch · 3 files changed

+206/−228 lines, 3 filesdep ~hmatrix

Dependency ranges changed: hmatrix

Files

repa-linear-algebra.cabal view
@@ -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
src/Numeric/LinearAlgebra/Repa.hs view
@@ -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
src/Numeric/LinearAlgebra/Repa/Conversion.hs view
@@ -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.