diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2015, Marcin Mrotek + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Marcin Mrotek nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# repa-linear-algebra +HMatrix operations for Repa. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/repa-linear-algebra.cabal b/repa-linear-algebra.cabal new file mode 100644 --- /dev/null +++ b/repa-linear-algebra.cabal @@ -0,0 +1,38 @@ +name: repa-linear-algebra +version: 0.0.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 +license-file: LICENSE +author: Marcin Mrotek +maintainer: marcin.jan.mrotek@gmail.com +stability: experimental (casts under the hood) +homepage: https://github.com/marcinmrotek/repa-linear-algebra +bug-reports: https://github.com/marcinmrotek/repa-linear-algebra/issues +category: Data +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 +source-repository head + type: git + location: https://github.com/marcinmrotek/repa-linear-algebra.git + +Flag Devel + Description: Development mode (-Werror). + Default: False + Manual: True + +library + exposed-modules: Numeric.LinearAlgebra.Repa + Numeric.LinearAlgebra.Repa.Conversion + -- other-modules: + -- other-extensions: + build-depends: base >=4.8 && <4.9 + , hmatrix + , repa + , vector + hs-source-dirs: src + ghc-options: -Wall + if flag(Devel) + ghc-options: -Werror + default-language: Haskell2010 diff --git a/src/Numeric/LinearAlgebra/Repa.hs b/src/Numeric/LinearAlgebra/Repa.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearAlgebra/Repa.hs @@ -0,0 +1,1464 @@ +{-# LANGUAGE FlexibleContexts #-} + +module Numeric.LinearAlgebra.Repa + ( Numeric + , Field + , Product + , 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) diff --git a/src/Numeric/LinearAlgebra/Repa/Conversion.hs b/src/Numeric/LinearAlgebra/Repa/Conversion.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/LinearAlgebra/Repa/Conversion.hs @@ -0,0 +1,146 @@ +{-# LANGUAGE + FlexibleContexts + , FlexibleInstances + , TypeFamilies + , TypeSynonymInstances + #-} + +module Numeric.LinearAlgebra.Repa.Conversion + ( Container + , Element + , HShape(..) + -- * Vector conversion utilities + , hv2repa + , repa2hv + , repa2hvS + , repa2hvSIO + , repa2hvP + , repa2hvPIO + -- * Matrix conversion utilities + , hm2repa + , repa2hm + , repa2hmS + , repa2hmSIO + , repa2hmP + , repa2hmPIO + ) where + +import Data.Array.Repa +import Data.Array.Repa.Repr.ForeignPtr +import Foreign.Storable +import Foreign.ForeignPtr +import qualified Numeric.LinearAlgebra.HMatrix as H +import Numeric.LinearAlgebra.HMatrix (Container, Element, Numeric) +import qualified Data.Vector.Storable as V + +-- |Shape-polymorphic conversion. +class HShape sh where + type HType sh :: * -> * + toRepa :: Numeric t => HType sh t -> Array F sh t + fromRepa :: Numeric t => Array F sh t -> HType sh t + fromRepaS :: Numeric t => Array D sh t -> HType sh t + fromRepaSIO :: Numeric t => Array D sh t -> IO (HType sh t) + fromRepaP :: (Numeric t, Monad m) => Array D sh t -> m (HType sh t) + fromRepaPIO :: Numeric t => Array D sh t -> IO (HType sh t) + +instance HShape DIM1 where + type HType DIM1 = H.Vector + toRepa = hv2repa + fromRepa = repa2hv + fromRepaS = repa2hvS + fromRepaSIO = repa2hvSIO + fromRepaP = repa2hvP + fromRepaPIO = repa2hvPIO + +instance HShape DIM2 where + type HType DIM2 = H.Matrix + toRepa = hm2repa + fromRepa = repa2hm + fromRepaS = repa2hmS + fromRepaSIO = repa2hmSIO + fromRepaP = repa2hmP + fromRepaPIO = repa2hmPIO + +-- Vector conversion utilities + +hv2repa :: Storable t => H.Vector t -> Array F DIM1 t +-- ^O(1). Convert a HMatrix Vector to a Repa Array. +hv2repa hv = fromForeignPtr (ix1 ln) ptr + where (ptr, ln) = V.unsafeToForeignPtr0 hv + +repa2hv :: Storable t => Array F DIM1 t -> H.Vector t +-- ^O(1). Convert a Repa Array to a HMatrix Vector. +repa2hv r = V.unsafeFromForeignPtr0 (toForeignPtr r) ln + where ln = size $ extent r + +repa2hvS :: Storable t => Array D DIM1 t -> H.Vector t +-- ^Convert a delayed Repa Array to a HMatrix Vector, evaluating it sequentially. +repa2hvS = repa2hv.computeS + +repa2hvSIO :: Storable t => Array D DIM1 t -> IO (H.Vector t) +-- ^O(1). Convert a Repa Array to a HMatrix Vector sequentially inside the IO monad. +repa2hvSIO r = do + ptr <- mallocForeignPtrArray ln + computeIntoS ptr r + return $ V.unsafeFromForeignPtr0 ptr ln + where ln = size $ extent r + +repa2hvP :: (Storable t, Monad m) => Array D DIM1 t -> m (H.Vector t) +-- ^Convert a delayed Repa Array to a HMatrix Vector, evaluating it in parallel. +repa2hvP = fmap repa2hv.computeP + +repa2hvPIO :: Storable t => Array D DIM1 t -> IO (H.Vector t) +-- ^O(1). Convert a Repa Array to a HMatrix Vector in parallel inside the IO monad. +repa2hvPIO r = do + ptr <- mallocForeignPtrArray ln + computeIntoP ptr r + return $ V.unsafeFromForeignPtr0 ptr ln + where ln = size $ extent r + +-- Matrix conversion utilities + +hm2repa + :: ( Storable t + , Container V.Vector t + , Element t + ) + => H.Matrix t -> Array F DIM2 t +-- ^O(1). Convert a HMatrix Matrix to a Repa Array. +hm2repa hm = fromForeignPtr (ix2 r c) ptr + where (ptr, _) = V.unsafeToForeignPtr0 $ H.flatten hm + (r , c) = H.size hm + +repa2hm :: Storable t => Array F DIM2 t -> H.Matrix t +-- ^O(1). Convert a Repa Array to a HMatrix Matrix. +repa2hm r = H.reshape c $ V.unsafeFromForeignPtr0 (toForeignPtr r) ln + where ln = size e + (_:c:[]) = listOfShape e + e = extent r + +repa2hmS :: Storable t => Array D DIM2 t -> H.Matrix t +-- ^Convert a delayed Repa Array to a HMatrix Matrix, evaluating it sequentially. +repa2hmS = repa2hm . computeS + +repa2hmSIO :: Storable t => Array D DIM2 t -> IO (H.Matrix t) +-- ^O(1). Convert a Repa Array to a HMatrix Matrix sequentially inside the IO monad. +repa2hmSIO r = do + ptr <- mallocForeignPtrArray ln + computeIntoS ptr r + return . H.reshape c $ V.unsafeFromForeignPtr0 ptr ln + where ln = size e + (_:c:[]) = listOfShape e + e = extent r + +repa2hmP :: (Storable t, Monad m) => Array D DIM2 t -> m (H.Matrix t) +-- ^Convert a delayed Repa Array to a HMatrix Matrix, evaluating it in parallel. +repa2hmP = fmap repa2hm . computeP + +repa2hmPIO :: Storable t => Array D DIM2 t -> IO (H.Matrix t) +-- ^O(1). Convert a Repa Array to a HMatrix Matrix in parallel inside the IO monad. +repa2hmPIO r = do + ptr <- mallocForeignPtrArray ln + computeIntoP ptr r + return . H.reshape c $ V.unsafeFromForeignPtr0 ptr ln + where ln = size e + (_:c:[]) = listOfShape e + e = extent r