{-# LANGUAGE FlexibleContexts #-}

module Numeric.LinearAlgebra.Repa
  ( Numeric
  , Field
  , Product
  , Vector
  , H.Matrix
  , RandDist(..)
  , Seed
  , HShape(..)
  , LSDiv
  -- * Dot product
  , dot
  , dotS
  , dotSIO
  , dotP
  , dotPIO
  -- * Dense matrix-vector product.
  , app
  , appS
  , appSIO
  , appP
  , appPIO
  -- * Dense matrix-matrix product.
  , mul
  , mulS
  , mulSIO
  , mulP
  , mulPIO
  -- * Vector outer product.
  , outer
  , outerS
  , outerSIO
  , outerP
  , outerPIO
  -- * Kronecker product.
  , kronecker
  , kroneckerS
  , kroneckerSIO
  , kroneckerP
  , kroneckerPIO
  -- * Cross product.
  , cross
  , crossS
  , crossSIO
  -- * Sum of elements.
  , sumElements
  , sumElementsS
  , sumElementsSIO
  , sumElementsP
  , sumElementsPIO
  -- * Product of elements.
  , prodElements
  , prodElementsS
  , prodElementsSIO
  , prodElementsP
  , prodElementsPIO
  -- * Linear systems.
  , (<\>)
  , solve
  , solveS
  , solveSIO
  , solveP
  , solvePIO
  , linearSolve
  , linearSolveS
  , linearSolveSIO
  , linearSolveP
  , linearSolvePIO
  , linearSolveLS
  , linearSolveLS_S
  , linearSolveLS_SIO
  , linearSolveLS_P
  , linearSolveLS_PIO
  , linearSolveSVD
  , linearSolveSVD_S
  , linearSolveSVD_SIO
  , linearSolveSVD_P
  , linearSolveSVD_PIO
  , luSolve
  , luSolveS
  , luSolveSIO
  , luSolveP
  , luSolvePIO
  , cholSolve
  , cholSolveS
  , cholSolveSIO
  , cholSolveP
  , cholSolvePIO
  -- * Inverse and pseudoinverse
  , inv
  , invS
  , invSIO
  , invP
  , invPIO
  , pinv
  , pinvS
  , pinvSIO
  , pinvP
  , pinvPIO
  , pinvTol
  , pinvTolS
  , pinvTolSIO
  , pinvTolP
  , pinvTolPIO
  -- * Determinant and rank
  , rcond
  , rcondS
  , rcondSIO
  , rcondP
  , rcondPIO
  , rank
  , rankS
  , rankSIO
  , rankP
  , rankPIO
  , det
  , detS
  , detSIO
  , detP
  , detPIO
  , invlndet
  , invlndetS
  , invlndetSIO
  , invlndetP
  , invlndetPIO
  -- * Norms
  , norm_Frob
  , norm_FrobS
  , norm_FrobSIO
  , norm_FrobP
  , norm_FrobPIO
  , norm_nuclear
  , norm_nuclearS
  , norm_nuclearSIO
  , norm_nuclearP
  , norm_nuclearPIO
  -- * Nullspace and range
  , orth
  , orthS
  , orthSIO
  , orthP
  , orthPIO
  , nullspace
  , nullspaceS
  , nullspaceSIO
  , nullspaceP
  , nullspacePIO
  , null1
  , null1S
  , null1SIO
  , null1P
  , null1PIO
  , null1sym
  , null1symS
  , null1symSIO
  , null1symP
  , null1symPIO
  -- * SVD
  , svd
  , svdS
  , svdSIO
  , svdP
  , svdPIO
  , thinSVD
  , thinSVD_S
  , thinSVD_SIO
  , thinSVD_P
  , thinSVD_PIO
  , compactSVD
  , compactSVD_S
  , compactSVD_SIO
  , compactSVD_P
  , compactSVD_PIO
  , singularValues
  , singularValuesS
  , singularValuesSIO
  , singularValuesP
  , singularValuesPIO
  , leftSV
  , leftSV_S
  , leftSV_SIO
  , leftSV_P
  , leftSV_PIO
  , rightSV
  , rightSV_S
  , rightSV_SIO
  , rightSV_P
  , rightSV_PIO
  -- * Eigensystems
  , eig
  , eigS
  , eigSIO
  , 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
  -- * QR
  , qr
  , qrS
  , qrSIO
  , qrP
  , qrPIO
  , rq
  , rqS
  , rqSIO
  , rqP
  , rqPIO
  , qrRaw
  , qrRawS
  , qrRawSIO
  , qrRawP
  , qrRawPIO
  , qrgr
  -- * Cholesky
  , chol
  , cholS
  , cholSIO
  , cholP
  , cholPIO
  , chol'
  , chol'S
  , chol'SIO
  , chol'P
  , chol'PIO
  -- * Hessenberg
  , hess
  , hessS
  , hessSIO
  , hessP
  , hessPIO
  -- * Schur
  , schur
  , schurS
  , schurSIO
  , schurP
  , schurPIO
  -- * LU
  , lu
  , luS
  , luSIO
  , luP
  , luPIO
  , luPacked
  , luPackedS
  , luPackedSIO
  , luPackedP
  , luPackedPIO
  -- * Matrix functions
  , expm
  , expmS
  , expmSIO
  , expmP
  , expmPIO
  , sqrtm
  , sqrtmS
  , sqrtmSIO
  , sqrtmP
  , sqrtmPIO
  , matFunc
  , matFuncS
  , matFuncSIO
  , matFuncP
  , matFuncPIO
  -- *Correlation and convolution
  , corr
  , corrS
  , corrSIO
  , corrP
  , corrPIO
  , conv
  , convS
  , convSIO
  , convP
  , convPIO
  , corrMin
  , corrMinS
  , corrMinSIO
  , corrMinP
  , corrMinPIO
  , corr2
  , corr2S
  , corr2SIO
  , corr2P
  , corr2PIO
  , conv2
  , conv2S
  , conv2SIO
  , conv2P
  , conv2PIO
  -- *Random vectors and matrices
  , randomVector
  , randomMatrix
  , gaussianSample
  , uniformSample
  -- *Misc
  , meanCov
  , rowOuters
  ) where

import Numeric.LinearAlgebra.Repa.Conversion

import Data.Array.Repa hiding (rank)
import Data.Array.Repa.Repr.ForeignPtr
import qualified Numeric.LinearAlgebra.HMatrix as H
import Numeric.LinearAlgebra.HMatrix (Complex, Numeric, Field, LSDiv, Normed, Product, Vector, RealElement, RandDist(..), Seed)

-- Dot product

dot :: Numeric t => Array F DIM1 t -> Array F DIM1 t -> t
-- ^Vector dot product.
dot v u = repa2hv v `H.dot` repa2hv u

dotS :: Numeric t => Array D DIM1 t -> Array D DIM1 t -> t
-- ^Vector dot product. Arguments computed sequentially.
dotS v u = repa2hvS v `H.dot` repa2hvS u

dotSIO :: Numeric t => Array D DIM1 t -> Array D DIM1 t -> IO t
-- ^Vector dot product. Arguments computed sequentially inside the IO monad.
dotSIO v u = H.dot <$> repa2hvSIO v <*> repa2hvSIO u

dotP :: (Numeric t, Monad m) => Array D DIM1 t -> Array D DIM1 t -> m t
-- ^Vector dot product. Arguments computed in parallel.
dotP v u = H.dot <$> repa2hvP v <*> repa2hvP u

dotPIO :: Numeric t => Array D DIM1 t -> Array D DIM1 t -> IO t
-- ^Vector dot product. Arguments computed in parallel inside the IO monad.
dotPIO v u = H.dot <$> repa2hvPIO v <*> repa2hvPIO u

-- Dense matrix-vector product

app :: Numeric t => Array F DIM2 t -> Array F DIM1 t -> Array F DIM1 t
-- ^Dense matrix-vector product.
app m v = hv2repa $ repa2hm m `H.app` repa2hv v

appS :: Numeric t => Array D DIM2 t -> Array D DIM1 t -> Array F DIM1 t
-- ^Dense matrix-vector product. Arguments computed sequentially.
appS m v = hv2repa $ repa2hmS m `H.app` repa2hvS v

appSIO :: Numeric t => Array D DIM2 t -> Array D DIM1 t -> IO (Array F DIM1 t)
-- ^Dense matrix-vector product. Arguments computed sequentially inside the IO monad.
appSIO m v = hv2repa <$> (H.app <$> repa2hmSIO m <*> repa2hvSIO v)

appP :: (Numeric t, Monad m) => Array D DIM2 t -> Array D DIM1 t -> m (Array F DIM1 t)
-- ^Dense matrix-vector product. Arguments computed in parallel.
appP m v = hv2repa <$> (H.app <$> repa2hmP m <*> repa2hvP v)

appPIO :: Numeric t => Array D DIM2 t -> Array D DIM1 t -> IO (Array F DIM1 t)
-- ^Dense matrix-vector product. Arguments computed in parallel inside the IO monad.
appPIO m v = hv2repa <$> (H.app <$> repa2hmPIO m <*> repa2hvPIO v)

-- Dense matrix-matrix product

mul :: Numeric t => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^Dense matrix-matrix product.
mul m n = hm2repa $ repa2hm m `H.mul` repa2hm n

mulS :: Numeric t => Array D DIM2 t -> Array D DIM2 t -> Array F DIM2 t
-- ^Dense matrix-matrix product. Arguments computed sequentially.
mulS m n = hm2repa $ repa2hmS m `H.mul` repa2hmS n

mulSIO :: Numeric t => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
-- ^Dense matrix-matrix product. Arguments computed sequentially inside the IO monad
mulSIO m n = hm2repa <$> (H.mul <$> repa2hmSIO m <*> repa2hmSIO n)

mulP :: (Numeric t, Monad m) => Array D DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
-- ^Dense matrix-matrix product. Arguments computed in parallel.
mulP m n = hm2repa <$> (H.mul <$> repa2hmP m <*> repa2hmP n)

mulPIO :: Numeric t => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
-- ^Dense matrix-matrix product. Arguments computed in parallel inside the IO monad
mulPIO m n = hm2repa <$> (H.mul <$> repa2hmPIO m <*> repa2hmPIO n)

-- Outer product of two vectors

outer :: (Product t, Numeric t) => Array F DIM1 t -> Array F DIM1 t -> Array F DIM2 t
-- |Outer product of two vectors.
outer v u = hm2repa $ repa2hv v `H.outer` repa2hv u

outerS :: (Product t, Numeric t) => Array D DIM1 t -> Array D DIM1 t -> Array F DIM2 t
-- |Outer product of two vectors. Arguments computed sequentially.
outerS v u = hm2repa $ repa2hvS v `H.outer` repa2hvS u

outerSIO :: (Product t, Numeric t) => Array D DIM1 t -> Array D DIM1 t -> IO (Array F DIM2 t)
-- |Outer product of two vectors. Arguments computed sequentially inside the IO monad.
outerSIO v u = hm2repa <$> (H.outer <$> repa2hvSIO v <*> repa2hvSIO u)

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)

-- Kronecker product of two matrices

kronecker :: (Product t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^Kronecker product of two matrices.
kronecker m n = hm2repa $ repa2hm m `H.kronecker` repa2hm n

kroneckerS :: (Product t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> Array F DIM2 t
-- ^Kronecker product of two matrices. Arguments computed sequentially.
kroneckerS m n = hm2repa $ repa2hmS m `H.kronecker` repa2hmS n

kroneckerSIO :: (Product t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
-- ^Kronecker product of two matrices. Arguments computed sequentially inside the IO monad.
kroneckerSIO m n = hm2repa <$> (H.kronecker <$> repa2hmSIO m <*> repa2hmSIO n)

kroneckerP :: (Product t, Numeric t, Monad m) => Array D DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
-- ^Kronecker product of two matrices. Arguments computed in parallel.
kroneckerP m n = hm2repa <$> (H.kronecker <$> repa2hmP m <*> repa2hmP n)

kroneckerPIO :: (Product t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
-- ^Kronecker product of two matrices. Arguments computed in parallel inside the IO monad.
kroneckerPIO m n = hm2repa <$> (H.kronecker <$> repa2hmPIO m <*> repa2hmPIO n)

-- Cross product

cross :: Array F DIM1 Double -> Array F DIM1 Double -> Array F DIM1 Double
-- ^Vector cross product.
cross v u = hv2repa $ repa2hv v `H.cross` repa2hv u

crossS :: Array D DIM1 Double -> Array D DIM1 Double -> Array F DIM1 Double
-- ^Vector cross product. Arguments computed sequentially.
crossS v u = hv2repa $ repa2hvS v `H.cross` repa2hvS u

crossSIO :: Array D DIM1 Double -> Array D DIM1 Double -> IO (Array F DIM1 Double)
-- ^Vector cross product. Arguments computed sequentially inside the IO monad.
crossSIO v u = hv2repa <$> (H.cross <$> repa2hvSIO v <*> repa2hvSIO u)

-- Sum of elements

sumElements :: (Numeric t, HShape sh, Container (HType sh) t) => Array F sh t -> t
-- ^Sum elements of a matrix or a vector.
sumElements = H.sumElements . fromRepa

sumElementsS :: (Numeric t, HShape sh, Container (HType sh) t) => Array D sh t -> t
-- ^Sum elements of a matrix or a vector. Argument computed sequentially.
sumElementsS = H.sumElements . fromRepaS

sumElementsSIO :: (Numeric t, HShape sh, Container (HType sh) t) => Array D sh t -> IO t
-- ^Sum elements of a matrix or a vector. Argument computed sequentially in the IO monad.
sumElementsSIO = fmap H.sumElements . fromRepaSIO

sumElementsP :: (Numeric t, HShape sh, Container (HType sh) t, Monad m) => Array D sh t -> m t
-- ^Sum elements of a matrix or a vector. Argument computed in parallel.
sumElementsP = fmap H.sumElements . fromRepaP

sumElementsPIO :: (Numeric t, HShape sh, Container (HType sh) t) => Array D sh t -> IO t
-- ^Sum elements of a matrix or a vector. Argument computed in parallel in the IO monad.
sumElementsPIO = fmap H.sumElements . fromRepaPIO

-- Product of elements

prodElements :: (Numeric t, HShape sh, Container (HType sh) t) => Array F sh t -> t
-- ^Multiply elements of a matrix or a vector.
prodElements = H.prodElements . fromRepa

prodElementsS :: (Numeric t, HShape sh, Container (HType sh) t) => Array D sh t -> t
-- ^Multiply elements of a matrix or a vector. Argument computed sequentially.
prodElementsS = H.prodElements . fromRepaS

prodElementsSIO :: (Numeric t, HShape sh, Container (HType sh) t) => Array D sh t -> IO t
-- ^Multiply elements of a matrix or a vector. Argument computed sequentially inside the IO monad.
prodElementsSIO = fmap H.prodElements . fromRepaSIO

prodElementsP :: (Numeric t, HShape sh, Container (HType sh) t, Monad m) => Array D sh t -> m t
-- ^Multiply elements of a matrix or a vector. Argument computed in parallel.
prodElementsP = fmap H.prodElements . fromRepaP

prodElementsPIO :: (Numeric t, HShape sh, Container (HType sh) t) => Array D sh t -> IO t
-- ^Multiply elements of a matrix or a vector. Argument computed in parallel inside the IO monad.
prodElementsPIO = fmap H.prodElements . fromRepaPIO

-- Linear systems.

(<\>) :: (Field t, Numeric t, HShape sh, LSDiv (HType sh)) => Array F DIM2 t -> Array F sh t -> Array F sh t
-- ^Infix alias for 'solve'.
(<\>) = solve

solve :: (Field t, Numeric t, HShape sh, LSDiv (HType sh)) => Array F DIM2 t -> Array F sh t -> Array F sh t
-- ^Least squares solution of a linear system, similar to the \ operator of Matlab/Octave (based on linearSolveSD).
solve m n = toRepa $ repa2hm m H.<\> fromRepa n

solveS :: (Field t, Numeric t, HShape sh, LSDiv (HType sh)) => Array D DIM2 t -> Array D sh t -> Array F sh t
-- ^Least squares solution of a linear system, similar to the \ operator of Matlab/Octave (based on linearSolveSD). Arguments are computed sequentially.
solveS m n = toRepa $ repa2hmS m H.<\> fromRepaS n

solveSIO :: (Field t, Numeric t, HShape sh, LSDiv (HType sh)) => Array D DIM2 t -> Array D sh t -> IO (Array F sh t)
-- ^Least squares solution of a linear system, similar to the \ operator of Matlab/Octave (based on linearSolveSD). Arguments are computed sequentially inside the IO monad.
solveSIO m n = toRepa <$> ((H.<\>) <$> repa2hmSIO m <*> fromRepaSIO n)

solveP :: (Field t, Numeric t, HShape sh, LSDiv (HType sh), Monad m) => Array D DIM2 t -> Array D sh t -> m (Array F sh t)
-- ^Least squares solution of a linear system, similar to the \ operator of Matlab/Octave (based on linearSolveSD). Arguments are computed in parallel.
solveP m n = toRepa <$> ((H.<\>) <$> repa2hmP m <*> fromRepaP n)

solvePIO :: (Field t, Numeric t, HShape sh, LSDiv (HType sh)) => Array D DIM2 t -> Array D sh t -> IO (Array F sh t)
-- ^Least squares solution of a linear system, similar to the \ operator of Matlab/Octave (based on linearSolveSD). Arguments are computed in parallel inside the IO monad.
solvePIO m n = toRepa <$> ((H.<\>) <$> repa2hmPIO m <*> fromRepaPIO n)


linearSolve :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Maybe (Array F DIM2 t)
-- ^Solve a linear system (for square coefficient matrix and several right hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
linearSolve m n = hm2repa <$> H.linearSolve (repa2hm m) (repa2hm n)

linearSolveS :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> Maybe (Array F DIM2 t)
-- ^Solve a linear system using the LU decomposition. Arguments computed sequentially.
linearSolveS m n = hm2repa <$> H.linearSolve (repa2hmS m) (repa2hmS n)

linearSolveP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> Array D DIM2 t -> m (Maybe (Array F DIM2 t))
-- ^Solve a linear system using the LU decomposition. Arguments computed in parallel.
linearSolveP m n = (hm2repa <$>) <$> (H.linearSolve <$> repa2hmP m <*> repa2hmP n)

linearSolveSIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Maybe (Array F DIM2 t))
-- ^Solve a linear system using the LU decomposition. Arguments computed sequentially inside the IO monad.
linearSolveSIO m n = (hm2repa <$>) <$> (H.linearSolve <$> repa2hmSIO m <*> repa2hmSIO n)

linearSolvePIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Maybe (Array F DIM2 t))
-- ^Solve a linear system using the LU decomposition. Arguments computed in parallel inside the IO monad.
linearSolvePIO m n = (hm2repa <$>) <$> (H.linearSolve <$> repa2hmPIO m <*> repa2hmPIO n)


linearSolveLS :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^Least squared error solution of an overcompensated system, or the minimum norm solution of an undercompensated system. For rank-deficient systems use 'linearSolveSVD'.
linearSolveLS m n = hm2repa $ H.linearSolveLS (repa2hm m) (repa2hm n)

linearSolveLS_S :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> Array F DIM2 t
linearSolveLS_S m n = hm2repa $ H.linearSolveLS (repa2hmS m) (repa2hmS n)

linearSolveLS_SIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
linearSolveLS_SIO m n = hm2repa <$> (H.linearSolveLS <$> repa2hmSIO m <*> repa2hmSIO n)

linearSolveLS_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
linearSolveLS_P m n = hm2repa <$> (H.linearSolveLS <$> repa2hmP m <*> repa2hmP n)

linearSolveLS_PIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
linearSolveLS_PIO m n = hm2repa <$> (H.linearSolveLS <$> repa2hmPIO m <*> repa2hmPIO n)


linearSolveSVD :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^Minimum norm solution of a general linear least squares problem Ax=b using the SVD. Admits rank-deficient systems but is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular values which are less than eps times the largest singular value.
linearSolveSVD m n = hm2repa $ H.linearSolveSVD (repa2hm m) (repa2hm n)

linearSolveSVD_S :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> Array F DIM2 t
linearSolveSVD_S m n = hm2repa $ H.linearSolveSVD (repa2hmS m) (repa2hmS n)

linearSolveSVD_SIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
linearSolveSVD_SIO m n = hm2repa <$> (H.linearSolveSVD <$> repa2hmSIO m <*> repa2hmSIO n)

linearSolveSVD_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
linearSolveSVD_P m n = hm2repa <$> (H.linearSolveSVD <$> repa2hmP m <*> repa2hmP n)

linearSolveSVD_PIO :: (Field t, Numeric t) => Array D DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
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
-- ^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)

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)

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

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

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


cholSolve :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^Solve a symmetric or Herimitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
cholSolve ch m = hm2repa $ H.cholSolve (repa2hm ch) (repa2hm m)

cholSolveS :: (Field t, Numeric t) => Array F DIM2 t -> Array D DIM2 t -> Array F DIM2 t
cholSolveS ch m = hm2repa $ H.cholSolve (repa2hm ch) (repa2hmS m)

cholSolveSIO :: (Field t, Numeric t) => Array F DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
cholSolveSIO ch m = hm2repa . H.cholSolve (repa2hm ch) <$> repa2hmSIO m

cholSolveP :: (Field t, Numeric t, Monad m) => Array F DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
cholSolveP ch m = hm2repa . H.cholSolve (repa2hm ch) <$> repa2hmP m

cholSolvePIO :: (Field t, Numeric t) => Array F DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
cholSolvePIO ch m = hm2repa . H.cholSolve (repa2hm ch) <$> repa2hmPIO m

-- Inverse and pseudoinverse

inv :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
-- ^Inverse of a square matrix.
inv = hm2repa . H.inv . repa2hm

invS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
invS = hm2repa . H.inv . repa2hmS

invSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
invSIO = fmap (hm2repa . H.inv) . repa2hmSIO

invP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
invP = fmap (hm2repa . H.inv) . repa2hmP

invPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
invPIO = fmap (hm2repa . H.inv) . repa2hmPIO

pinv :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
-- ^Pseudoinverse of a general matrix, with default tolerance ('pinvTol' 1, similar to GNU-Octave)
pinv = hm2repa . H.pinv . repa2hm

pinvS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
pinvS = hm2repa . H.pinv . repa2hmS

pinvSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
pinvSIO = fmap (hm2repa . H.pinv) . repa2hmSIO

pinvP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
pinvP = fmap (hm2repa . H.pinv) . repa2hmP

pinvPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
pinvPIO = fmap (hm2repa . H.pinv) . repa2hmPIO

pinvTol :: (Field t, Numeric t) => Double -> Array F DIM2 t -> Array F DIM2 t
-- ^pinvTol r computes the pseudoinverse of a matrix with tolerance tol=r*g*eps*(max rows cols), where g is the greatest singular value.
pinvTol r = hm2repa . H.pinvTol r . repa2hm

pinvTolS :: (Field t, Numeric t) => Double -> Array D DIM2 t -> Array F DIM2 t
pinvTolS r = hm2repa . H.pinvTol r . repa2hmS

pinvTolSIO :: (Field t, Numeric t) => Double -> Array D DIM2 t -> IO (Array F DIM2 t)
pinvTolSIO r= fmap (hm2repa . H.pinvTol r) . repa2hmSIO

pinvTolP :: (Field t, Numeric t, Monad m) => Double -> Array D DIM2 t -> m (Array F DIM2 t)
pinvTolP r = fmap (hm2repa . H.pinvTol r) . repa2hmP

pinvTolPIO :: (Field t, Numeric t) => Double -> Array D DIM2 t -> IO (Array F DIM2 t)
pinvTolPIO r = fmap (hm2repa . H.pinvTol r) . repa2hmPIO

-- Determinant and rank

rcond :: (Field t, Numeric t) => Array F DIM2 t -> Double
-- ^Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
rcond = H.rcond . repa2hm

rcondS :: (Field t, Numeric t) => Array D DIM2 t -> Double
rcondS = H.rcond . repa2hmS

rcondSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO Double
rcondSIO = fmap H.rcond . repa2hmSIO

rcondP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m Double
rcondP = fmap H.rcond . repa2hmP

rcondPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO Double
rcondPIO = fmap H.rcond . repa2hmPIO

rank :: (Field t, Numeric t) => Array F DIM2 t -> Int
-- ^Number of linearly independent rows or columns. See also 'ranksv'.
rank = H.rank . repa2hm

rankS :: (Field t, Numeric t) => Array D DIM2 t -> Int
rankS = H.rank . repa2hmS

rankSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO Int
rankSIO = fmap H.rank . repa2hmSIO

rankP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m Int
rankP = fmap H.rank . repa2hmP

rankPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO Int
rankPIO = fmap H.rank . repa2hmPIO

det :: (Field t, Numeric t) => Array F DIM2 t -> t
-- ^Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
det = H.det . repa2hm

detS :: (Field t, Numeric t) => Array D DIM2 t -> t
detS = H.det . repa2hmS

detSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO t
detSIO = fmap H.det . repa2hmSIO

detP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m t
detP = fmap H.det . repa2hmP

detPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO t
detPIO = fmap H.det . repa2hmPIO

invlndet :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, (t, t)) -- ^(inverse, (log abs det, sign or phase of det))
-- ^Joint computation of inverse and logarithm of determinant of a square matrix.
invlndet m = let (h, r) = H.invlndet $ repa2hm m in (hm2repa h, r)

invlndetS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, (t, t))
invlndetS m = let (h, r) = H.invlndet $ repa2hmS m in (hm2repa h, r)

invlndetSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, (t, t))
invlndetSIO m = do
  (h, r) <- H.invlndet <$> repa2hmSIO m
  return (hm2repa h, r)

invlndetP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, (t, t))
invlndetP m = do
  (h, r) <- H.invlndet <$> repa2hmP m
  return (hm2repa h, r)

invlndetPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, (t, t))
invlndetPIO m = do
  (h, r) <- H.invlndet <$> repa2hmPIO m
  return (hm2repa h, r)

-- Norms

norm_Frob :: (Normed (Vector t), Element t) => Array F DIM2 t -> Double
norm_Frob = H.norm_Frob . repa2hm

norm_FrobS :: (Normed (Vector t), Element t) => Array D DIM2 t -> Double
norm_FrobS = H.norm_Frob . repa2hmS

norm_FrobSIO :: (Normed (Vector t), Element t) => Array D DIM2 t -> IO Double
norm_FrobSIO = fmap H.norm_Frob . repa2hmSIO

norm_FrobP :: (Normed (Vector t), Element t, Monad m) => Array D DIM2 t -> m Double
norm_FrobP = fmap H.norm_Frob . repa2hmP

norm_FrobPIO :: (Normed (Vector t), Element t) => Array D DIM2 t -> IO Double
norm_FrobPIO = fmap H.norm_Frob . repa2hmPIO


norm_nuclear :: (Field t, Numeric t) => Array F DIM2 t -> Double
norm_nuclear = H.norm_nuclear . repa2hm

norm_nuclearS :: (Field t, Numeric t) => Array D DIM2 t -> Double
norm_nuclearS = H.norm_nuclear . repa2hmS

norm_nuclearSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO Double
norm_nuclearSIO = fmap H.norm_nuclear . repa2hmSIO

norm_nuclearP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m Double
norm_nuclearP = fmap H.norm_nuclear . repa2hmP

norm_nuclearPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO Double
norm_nuclearPIO = fmap H.norm_nuclear . repa2hmPIO

-- Nullspace and range

orth :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
-- ^An orthonormal basis of the range space of a matrix. See also 'orthSVD'.
orth = hm2repa . H.orth . repa2hm

orthS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
orthS = hm2repa . H.orth . repa2hmS

orthSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
orthSIO = fmap (hm2repa . H.orth) . repa2hmSIO

orthP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
orthP = fmap (hm2repa . H.orth) . repa2hmP

orthPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
orthPIO = fmap (hm2repa . H.orth) . repa2hmPIO

nullspace :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
-- ^An orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'.
nullspace = hm2repa . H.nullspace . repa2hm

nullspaceS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
nullspaceS = hm2repa . H.nullspace . repa2hmS

nullspaceSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
nullspaceSIO = fmap (hm2repa . H.nullspace) . repa2hmSIO

nullspaceP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
nullspaceP = fmap (hm2repa . H.nullspace) . repa2hmP

nullspacePIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
nullspacePIO = fmap (hm2repa . H.nullspace) . repa2hmPIO

null1 :: Array F DIM2 Double -> Array F DIM1 Double
-- ^Solution of an overconstrained homogenous linear system.
null1 = hv2repa . H.null1 . repa2hm

null1S :: Array D DIM2 Double -> Array F DIM1 Double
null1S = hv2repa . H.null1 . repa2hmS

null1SIO :: Array D DIM2 Double -> IO (Array F DIM1 Double)
null1SIO = fmap (hv2repa . H.null1) . repa2hmSIO

null1P :: Monad m => Array D DIM2 Double -> m (Array F DIM1 Double)
null1P = fmap (hv2repa . H.null1) . repa2hmP

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
-- ^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

-- SVD

svd :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
-- ^Full singular value decomposition.
svd m = let (u,s,v) = H.svd $ repa2hm m in (hm2repa u, hv2repa s, hm2repa v)

svdS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
svdS m = let (u,s,v) = H.svd $ repa2hmS m in (hm2repa u, hv2repa s, hm2repa v)

svdSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
svdSIO m = do
  (u,s,v) <- H.svd <$> repa2hmSIO m
  return (hm2repa u, hv2repa s, hm2repa v)

svdP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
svdP m = do
  (u,s,v) <- H.svd <$> repa2hmP m
  return (hm2repa u, hv2repa s, hm2repa v)

svdPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
svdPIO m = do
  (u,s,v) <- H.svd <$> repa2hmPIO m
  return (hm2repa u, hv2repa s, hm2repa v)

thinSVD :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
-- ^A version of 'svd' which returns only the min (rows m) (cols m) singular vectors of m. (u,s,v) = thinSVD m ==> m == u * diag s * tr v
thinSVD m = let (u,s,v) = H.thinSVD $ repa2hm m in (hm2repa u, hv2repa s, hm2repa v)

thinSVD_S :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
thinSVD_S m = let (u,s,v) = H.thinSVD $ repa2hmS m in (hm2repa u, hv2repa s, hm2repa v)

thinSVD_SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
thinSVD_SIO m = do
  (u,s,v) <- H.thinSVD <$> repa2hmSIO m
  return (hm2repa u, hv2repa s, hm2repa v)

thinSVD_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
thinSVD_P m = do
  (u,s,v) <- H.thinSVD <$> repa2hmP m
  return (hm2repa u, hv2repa s, hm2repa v)

thinSVD_PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
thinSVD_PIO m = do
  (u,s,v) <- H.thinSVD <$> repa2hmPIO m
  return (hm2repa u, hv2repa s, hm2repa v)

compactSVD :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
-- ^Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
compactSVD m = let (u,s,v) = H.compactSVD $ repa2hm m in (hm2repa u, hv2repa s, hm2repa v)

compactSVD_S :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
compactSVD_S m = let (u,s,v) = H.compactSVD $ repa2hmS m in (hm2repa u, hv2repa s, hm2repa v)

compactSVD_SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
compactSVD_SIO m = do
  (u,s,v) <- H.compactSVD <$> repa2hmSIO m
  return (hm2repa u, hv2repa s, hm2repa v)

compactSVD_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
compactSVD_P m = do
  (u,s,v) <- H.compactSVD <$> repa2hmP m
  return (hm2repa u, hv2repa s, hm2repa v)

compactSVD_PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double, Array F DIM2 t)
compactSVD_PIO m = do
  (u,s,v) <- H.compactSVD <$> repa2hmPIO m
  return (hm2repa u, hv2repa s, hm2repa v)

singularValues :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM1 Double
-- ^Singular values only.
singularValues = hv2repa . H.singularValues . repa2hm

singularValuesS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM1 Double
singularValuesS = hv2repa . H.singularValues . repa2hmS

singularValuesSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double)
singularValuesSIO = fmap (hv2repa . H.singularValues) . repa2hmSIO

singularValuesP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 Double)
singularValuesP = fmap (hv2repa . H.singularValues) . repa2hmP

singularValuesPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double)
singularValuesPIO = fmap (hv2repa . H.singularValues) . repa2hmPIO

leftSV :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM1 Double)
-- ^Singular values and all left singular vectors (as columns).
leftSV m = let (u,s) = H.leftSV $ repa2hm m in (hm2repa u, hv2repa s)

leftSV_S :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM1 Double)
leftSV_S m = let (u,s) = H.leftSV $ repa2hmS m in (hm2repa u, hv2repa s)

leftSV_SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double)
leftSV_SIO m = do
  (u,s) <- H.leftSV <$> repa2hmSIO m
  return (hm2repa u, hv2repa s)

leftSV_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM1 Double)
leftSV_P m = do
  (u,s) <- H.leftSV <$> repa2hmP m
  return (hm2repa u, hv2repa s)

leftSV_PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM1 Double)
leftSV_PIO m = do
  (u,s) <- H.leftSV <$> repa2hmPIO m
  return (hm2repa u, hv2repa s)

rightSV :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
-- ^Singular values and all right singular vectors (as columns).
rightSV m = let (s,v) = H.rightSV $ repa2hm m in (hv2repa s, hm2repa v)

rightSV_S :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM1 Double, Array F DIM2 t)
rightSV_S m = let (s,v) = H.rightSV $ repa2hmS m in (hv2repa s, hm2repa v)

rightSV_SIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
rightSV_SIO m = do
  (s,v) <- H.rightSV <$> repa2hmSIO m
  return (hv2repa s, hm2repa v)

rightSV_P :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 Double, Array F DIM2 t)
rightSV_P m = do
  (s,v) <- H.rightSV <$> repa2hmP m
  return (hv2repa s, hm2repa v)

rightSV_PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 Double, Array F DIM2 t)
rightSV_PIO m = do
  (s,v) <- H.rightSV <$> repa2hmPIO m
  return (hv2repa s, hm2repa v)

-- Eigensystems

eig :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM1 (Complex Double), Array F DIM2 (Complex Double))
-- ^Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix. (s,v) = eig m ==> m * v = v == v <> diag s
eig m = let (s,v) = H.eig $ repa2hm m in (hv2repa s, hm2repa v)

eigS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM1 (Complex Double), Array F DIM2 (Complex Double))
eigS m = let (s,v) = H.eig $ repa2hmS m in (hv2repa s, hm2repa v)

eigSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 (Complex Double), Array F DIM2 (Complex Double))
eigSIO m = do
  (s,v) <- H.eig <$> repa2hmSIO m
  return (hv2repa s, hm2repa v)

eigP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 (Complex Double), Array F DIM2 (Complex Double))
eigP m = do
  (s,v) <- H.eig <$> repa2hmP m
  return (hv2repa s, hm2repa v)

eigPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 (Complex Double), Array F DIM2 (Complex Double))
eigPIO m = do
  (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)
-- ^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)

eigenvalues :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM1 (Complex Double)
-- ^Eigenvalues (not ordered) of a general square matrix.
eigenvalues = hv2repa . H.eigenvalues . repa2hm

eigenvaluesS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM1 (Complex Double)
eigenvaluesS = hv2repa . H.eigenvalues . repa2hmS

eigenvaluesSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM1 (Complex Double))
eigenvaluesSIO = fmap (hv2repa . H.eigenvalues) . repa2hmSIO

eigenvaluesP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM1 (Complex Double))
eigenvaluesP = fmap (hv2repa . H.eigenvalues) . repa2hmP

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
-- ^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)

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)

-- QR

qr :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
-- ^QR factorization. (q,r) = qr m ==> m = q * r where q is unitary and r is upper triangular.
qr m = let (q,r) = H.qr $ repa2hm m in (hm2repa q, hm2repa r)

qrS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
qrS m = let (q,r) = H.qr $ repa2hmS m in (hm2repa q, hm2repa r)

qrSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
qrSIO m = do
  (q,r) <- H.qr <$> repa2hmSIO m
  return (hm2repa q, hm2repa r)

qrP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM2 t)
qrP m = do
  (q,r) <- H.qr <$> repa2hmP m
  return (hm2repa q, hm2repa r)

qrPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
qrPIO m = do
  (q,r) <- H.qr <$> repa2hmPIO m
  return (hm2repa q, hm2repa r)

rq :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
-- ^RQ factorization. (r,q) = rq m ==> m = r * q where q is unitary and r is upper triangular.
rq m = let (r,q) = H.rq $ repa2hm m in (hm2repa r, hm2repa q)

rqS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
rqS m = let (r,q) = H.rq $ repa2hmS m in (hm2repa r, hm2repa q)

rqSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
rqSIO m = do
  (r,q) <- H.rq <$> repa2hmSIO m
  return (hm2repa r, hm2repa q)

rqP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM2 t)
rqP m = do
  (r,q) <- H.rq <$> repa2hmP m
  return (hm2repa r, hm2repa q)

rqPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
rqPIO m = do
  (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)

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)

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)

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)

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)

qrgr :: (Field t, Numeric t) => Int -> (Array F DIM2 t, Array F DIM1 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)

-- Cholesky

chol :: (Field t, Numeric t) => Array F DIM2 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'PIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
chol'PIO = fmap (hm2repa . H.cholSH) . repa2hmPIO

-- Hessenberg

hess :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
-- ^Hessenberg factorization. (p,h) == hess m ==> p * h * p' where p is unitary and h is in upper Hessenberg form (zero entries below the first subdiagonal).
hess m = let (p,h) = H.hess $ repa2hm m in (hm2repa p, hm2repa h)

hessS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
hessS m = let (p,h) = H.hess $ repa2hmS m in (hm2repa p, hm2repa h)

hessSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
hessSIO m = do
  (p,h) <- H.hess <$> repa2hmSIO m
  return (hm2repa p, hm2repa h)

hessP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM2 t)
hessP m = do
  (p,h) <- H.hess <$> repa2hmP m
  return (hm2repa p, hm2repa h)

hessPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
hessPIO m = do
  (p,h) <- H.hess <$> repa2hmPIO m
  return (hm2repa p, hm2repa h)

-- Schur
schur :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
-- ^Schur factorization. (u,s) = schur m ==> m == u * s * u' where u is unitary and s is a Schur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is upper triangular in 2x2 blocks.
schur m = let (u,s) = H.schur $ repa2hm m in (hm2repa u, hm2repa s)

schurS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM2 t)
schurS m = let (u,s) = H.schur $ repa2hmS m in (hm2repa u, hm2repa s)

schurSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
schurSIO m = do
  (u,s) <- H.schur <$> repa2hmSIO m
  return (hm2repa u, hm2repa s)

schurP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM2 t)
schurP m = do
  (u,s) <- H.schur <$> repa2hmP m
  return (hm2repa u, hm2repa s)

schurPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t)
schurPIO m = do
  (u,s) <- H.schur <$> repa2hmPIO m
  return (hm2repa u, hm2repa s)

-- LU

lu :: (Field t, Numeric t) => Array F DIM2 t -> (Array F DIM2 t, Array F DIM2 t, Array F DIM2 t, t)
-- ^Explicit LU factorization of a general matrix. (l,u,p,s) = lu m ==> m = p * l * u where l is lower triangular, u is upper triangular, p is a permutation matrix, and s is the signature of the permutation.
lu m = let (l,u,p,s) = H.lu $ repa2hm m in (hm2repa l, hm2repa u, hm2repa p, s)

luS :: (Field t, Numeric t) => Array D DIM2 t -> (Array F DIM2 t, Array F DIM2 t, Array F DIM2 t, t)
luS m = let (l,u,p,s) = H.lu $ repa2hmS m in (hm2repa l, hm2repa u, hm2repa p, s)

luSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t, Array F DIM2 t, t)
luSIO m = do
  (l,u,p,s) <- H.lu <$> repa2hmSIO m
  return (hm2repa l, hm2repa u, hm2repa p, s)

luP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t, Array F DIM2 t, Array F DIM2 t, t)
luP m = do
  (l,u,p,s) <- H.lu <$> repa2hmP m
  return (hm2repa l, hm2repa u, hm2repa p, s)

luPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t, Array F DIM2 t, Array F DIM2 t, t)
luPIO m = do
  (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
-- ^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

luPackedS :: (Field t, Numeric t) => Array D DIM2 t -> PackedLU t
luPackedS m = let (lu', is) = H.luPacked $ repa2hmS m in PackedLU lu' is

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

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

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

-- Matrix functions

expm :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
-- ^Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Val Loan, based on a scaled Pade approximation.
expm = hm2repa . H.expm . repa2hm

expmS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
expmS = hm2repa . H.expm . repa2hmS

expmSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
expmSIO = fmap (hm2repa . H.expm) . repa2hmSIO

expmP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
expmP = fmap (hm2repa . H.expm) . repa2hmP

expmPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
expmPIO = fmap (hm2repa . H.expm) . repa2hmPIO

sqrtm :: (Field t, Numeric t) => Array F DIM2 t -> Array F DIM2 t
-- ^Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. It only works with invertible matrices that have a real solution.
sqrtm = hm2repa . H.sqrtm . repa2hm

sqrtmS :: (Field t, Numeric t) => Array D DIM2 t -> Array F DIM2 t
sqrtmS = hm2repa . H.sqrtm . repa2hmS

sqrtmSIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
sqrtmSIO = fmap (hm2repa . H.sqrtm) . repa2hmSIO

sqrtmP :: (Field t, Numeric t, Monad m) => Array D DIM2 t -> m (Array F DIM2 t)
sqrtmP = fmap (hm2repa . H.sqrtm) . repa2hmP

sqrtmPIO :: (Field t, Numeric t) => Array D DIM2 t -> IO (Array F DIM2 t)
sqrtmPIO = fmap (hm2repa . H.sqrtm) . repa2hmPIO

matFunc :: (Complex Double -> Complex Double) -> Array F DIM2 (Complex Double) -> Array F DIM2 (Complex Double)
-- Generic matrix function for diagonalizable matrices.
matFunc f = hm2repa . H.matFunc f . repa2hm

matFuncS :: (Complex Double -> Complex Double) -> Array D DIM2 (Complex Double) -> Array F DIM2 (Complex Double)
matFuncS f = hm2repa . H.matFunc f . repa2hmS

matFuncSIO :: (Complex Double -> Complex Double) -> Array D DIM2 (Complex Double) -> IO (Array F DIM2 (Complex Double))
matFuncSIO f = fmap (hm2repa . H.matFunc f) . repa2hmSIO

matFuncP :: Monad m => (Complex Double -> Complex Double) -> Array D DIM2 (Complex Double) -> m (Array F DIM2 (Complex Double))
matFuncP f = fmap (hm2repa . H.matFunc f) . repa2hmP

matFuncPIO :: (Complex Double -> Complex Double) -> Array D DIM2 (Complex Double) -> IO (Array F DIM2 (Complex Double))
matFuncPIO f = fmap (hm2repa . H.matFunc f) . repa2hmPIO

-- Correlation and convolution

corr :: (Product t, Numeric t) => Array F DIM1 t -> Array F DIM1 t -> Array F DIM1 t
-- ^Correlation.
corr k = hv2repa . H.corr (repa2hv k) . repa2hv

corrS :: (Product t, Numeric t) => Array F DIM1 t -> Array D DIM1 t -> Array F DIM1 t
corrS k = hv2repa . H.corr (repa2hv k) . repa2hvS

corrSIO :: (Product t, Numeric t) => Array F DIM1 t -> Array D DIM1 t -> IO (Array F DIM1 t)
corrSIO k = fmap (hv2repa . H.corr (repa2hv k)) . repa2hvSIO

corrP :: (Product t, Numeric t, Monad m) => Array F DIM1 t -> Array D DIM1 t -> m (Array F DIM1 t)
corrP k = fmap (hv2repa . H.corr (repa2hv k)) . repa2hvP

corrPIO :: (Product t, Numeric t) => Array F DIM1 t -> Array D DIM1 t -> IO (Array F DIM1 t)
corrPIO k = fmap (hv2repa . H.corr (repa2hv k)) . repa2hvPIO

conv :: (Product t, Numeric t) => Array F DIM1 t -> Array F DIM1 t -> Array F DIM1 t
-- ^Convolution - 'corr' with reversed kernel and padded input, equivalent to polynomial multiplication.
conv k = hv2repa . H.conv (repa2hv k) . repa2hv

convS :: (Product t, Numeric t) => Array F DIM1 t -> Array D DIM1 t -> Array F DIM1 t
convS k = hv2repa . H.conv (repa2hv k) . repa2hvS

convSIO :: (Product t, Numeric t) => Array F DIM1 t -> Array D DIM1 t -> IO (Array F DIM1 t)
convSIO k = fmap (hv2repa . H.conv (repa2hv k)) . repa2hvSIO

convP :: (Product t, Numeric t, Monad m) => Array F DIM1 t -> Array D DIM1 t -> m (Array F DIM1 t)
convP k = fmap (hv2repa . H.conv (repa2hv k)) . repa2hvP

convPIO :: (Product t, Numeric t) => Array F DIM1 t -> Array D DIM1 t -> IO (Array F DIM1 t)
convPIO k = fmap (hv2repa . H.conv (repa2hv k)) . repa2hvPIO

corrMin :: (Product t, Numeric t, RealElement t) => Array F DIM1 t -> Array F DIM1 t -> Array F DIM1 t
-- ^Similar to 'corr' but using 'min' instead of (*).
corrMin k = hv2repa . H.corrMin (repa2hv k) . repa2hv

corrMinS :: (Product t, Numeric t, RealElement t) => Array F DIM1 t -> Array D DIM1 t -> Array F DIM1 t
corrMinS k = hv2repa . H.corrMin (repa2hv k) . repa2hvS

corrMinSIO :: (Product t, Numeric t, RealElement t) => Array F DIM1 t -> Array D DIM1 t -> IO (Array F DIM1 t)
corrMinSIO k = fmap (hv2repa . H.corrMin (repa2hv k)) . repa2hvSIO

corrMinP :: (Product t, Numeric t, RealElement t, Monad m) => Array F DIM1 t -> Array D DIM1 t -> m (Array F DIM1 t)
corrMinP k = fmap (hv2repa . H.corrMin (repa2hv k)) . repa2hvP

corrMinPIO :: (Product t, Numeric t, RealElement t) => Array F DIM1 t -> Array D DIM1 t -> IO (Array F DIM1 t)
corrMinPIO k = fmap (hv2repa . H.corrMin (repa2hv k)) . repa2hvPIO

corr2 :: (Product t, Numeric t) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^2D correlation (without padding).
corr2 k = hm2repa . H.corr2 (repa2hm k) . repa2hm

corr2S :: (Product t, Numeric t) => Array F DIM2 t -> Array D DIM2 t -> Array F DIM2 t
corr2S k = hm2repa . H.corr2 (repa2hm k) . repa2hmS

corr2SIO :: (Product t, Numeric t) => Array F DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
corr2SIO k = fmap (hm2repa . H.corr2 (repa2hm k)) . repa2hmSIO

corr2P :: (Product t, Numeric t, Monad m) => Array F DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
corr2P k = fmap (hm2repa . H.corr2 (repa2hm k)) . repa2hmP

corr2PIO :: (Product t, Numeric t) => Array F DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
corr2PIO k = fmap (hm2repa . H.corr2 (repa2hm k)) . repa2hmPIO

conv2 :: (Product t, Numeric t, Num (Vector t)) => Array F DIM2 t -> Array F DIM2 t -> Array F DIM2 t
-- ^2D convolution.
conv2 k = hm2repa . H.conv2 (repa2hm k) . repa2hm

conv2S :: (Product t, Numeric t, Num (Vector t)) => Array F DIM2 t -> Array D DIM2 t -> Array F DIM2 t
conv2S k = hm2repa . H.conv2 (repa2hm k) . repa2hmS

conv2SIO :: (Product t, Numeric t, Num (Vector t)) => Array F DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
conv2SIO k = fmap (hm2repa . H.conv2 (repa2hm k)) . repa2hmSIO

conv2P :: (Product t, Numeric t, Num (Vector t), Monad m) => Array F DIM2 t -> Array D DIM2 t -> m (Array F DIM2 t)
conv2P k = fmap (hm2repa . H.conv2 (repa2hm k)) . repa2hmP

conv2PIO :: (Product t, Numeric t, Num (Vector t)) => Array F DIM2 t -> Array D DIM2 t -> IO (Array F DIM2 t)
conv2PIO k = fmap (hm2repa . H.conv2 (repa2hm k)) . repa2hmPIO

-- Random arrays

randomVector :: Seed -> RandDist -> Int -> Array F DIM1 Double
-- ^Pseudorandom vector of na given size. Usee 'randomIO' to get a random seed.
randomVector s d n = hv2repa $ H.randomVector s d n

randomMatrix :: RandDist -> Int -> Int -> IO (Array F DIM2 Double)
randomMatrix d a b = fmap hm2repa
  $ case d of
       Uniform    -> H.rand a b
       Gaussian -> H.randn a b

gaussianSample :: Seed -> Int -> Array F DIM1 Double -> Array F DIM2 Double -> Array F DIM2 Double
-- ^A matrix whose rows are pseudorandom samples from a multivariate Gaussian distribution.
gaussianSample s r mean cov = hm2repa $ H.gaussianSample s r (repa2hv mean) (repa2hm cov)

uniformSample :: Seed -> Int -> [(Double,Double)] -> Array F DIM2 Double
-- ^A matrix whose rows are pseudorandom samples from a multivariate uniform distribution.
uniformSample s r rng = hm2repa $ H.uniformSample s r rng

-- misc

meanCov :: Array F DIM2 Double -> (Array F DIM1 Double, Array F DIM2 Double)
-- ^Compute mean vector and a covariance matrix of the rows of a matrix.
meanCov m = let (v,c) = H.meanCov $ repa2hm m in (hv2repa v, hm2repa c)

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)