hmatrix 0.19.0.0 → 0.20.0.0

raw patch · 14 files changed

+235/−44 lines, 14 filesdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: vector

API changes (from Hackage documentation)

Files

• hmatrix.cabal +3/−3
• src/Internal/Algorithms.hs +28/−4
• src/Internal/C/lapack-aux.c +87/−0
• src/Internal/Container.hs +18/−13
• src/Internal/Convolution.hs +4/−2
• src/Internal/LAPACK.hs +63/−1
• src/Internal/Matrix.hs +2/−1
• src/Internal/Sparse.hs +3/−1
• src/Internal/Static.hs +8/−8
• src/Internal/Util.hs +4/−2
• src/Internal/Vector.hs +8/−4
• src/Internal/Vectorized.hs +2/−1
• src/Numeric/LinearAlgebra.hs +4/−3
• src/Numeric/LinearAlgebra/Static.hs +1/−1

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.19.0.0+Version:            0.20.0.0 License:            BSD3 License-file:       LICENSE Author:             Alberto Ruiz@@ -38,7 +38,7 @@  library -    Build-Depends:      base >= 4.8 && < 5,+    Build-Depends:      base >= 4.9 && < 5,                         binary,                         array,                         deepseq,@@ -47,7 +47,7 @@                         bytestring,                         storable-complex,                         semigroups,-                        vector >= 0.8+                        vector >= 0.11      hs-source-dirs:     src 

@@ -70,8 +70,10 @@     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t     linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t     eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))+    geig'        :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double))     eigSH''      :: Matrix t -> (Vector Double, Matrix t)     eigOnly      :: Matrix t -> Vector (Complex Double)+    geigOnly     :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)     eigOnlySH    :: Matrix t -> Vector Double     cholSH'      :: Matrix t -> Matrix t     mbCholSH'    :: Matrix t -> Maybe (Matrix t)@@ -96,7 +98,9 @@     linearSolveSVD' = linearSolveSVDR Nothing     eig' = eigR     eigSH'' = eigS+    geig' = eigG     eigOnly = eigOnlyR+    geigOnly = eigOnlyG     eigOnlySH = eigOnlyS     cholSH' = cholS     mbCholSH' = mbCholS@@ -126,7 +130,9 @@     linearSolveLS' = linearSolveLSC     linearSolveSVD' = linearSolveSVDC Nothing     eig' = eigC+    geig' = eigGC     eigOnly = eigOnlyC+    geigOnly = eigOnlyGC     eigSH'' = eigH     eigOnlySH = eigOnlyH     cholSH' = cholH@@ -170,7 +176,8 @@ -0.690  -0.352   0.433  -0.233   0.398  >>> s-fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]+[35.18264833189422,1.4769076999800903,1.089145439970417e-15]+it :: Vector Double  >>> disp 3 v 3x3@@ -224,7 +231,8 @@ -0.690  -0.352   0.433  >>> s-fromList [35.18264833189422,1.4769076999800903,1.089145439970417e-15]+[35.18264833189422,1.4769076999800903,1.089145439970417e-15]+it :: Vector Double  >>> disp 3 v 3x3@@ -283,7 +291,8 @@ -0.690  -0.352  >>> s-fromList [35.18264833189422,1.4769076999800903]+[35.18264833189422,1.476907699980091]+it :: Vector Double  >>> disp 3 u 5x2@@ -509,10 +518,25 @@ eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) eig = {-# SCC "eig" #-} eig' +-- | Generalized eigenvalues (not ordered) and eigenvectors (as columns) of a pair of nonsymmetric matrices.+-- Eigenvalues are represented as pairs of alpha, beta, where eigenvalue = alpha / beta. Alpha is always+-- complex, but betas has the same type as the input matrix.+--+-- If @(alphas, betas, v) = geig a b@, then @a \<> v == b \<> v \<> diag (alphas / betas)@+--+-- Note that beta can be 0 and that has reasonable interpretation.+geig :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double))+geig = {-# SCC "geig" #-} geig'+ -- | Eigenvalues (not ordered) of a general square matrix. eigenvalues :: Field t => Matrix t -> Vector (Complex Double) eigenvalues = {-# SCC "eigenvalues" #-} eigOnly +-- | Generalized eigenvalues of a pair of matrices. Represented as pairs of alpha, beta,+-- where eigenvalue is alpha / beta as in 'geig'.+geigenvalues :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)+geigenvalues = {-# SCC "geigenvalues" #-} geigOnly+ -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) eigSH' = {-# SCC "eigSH'" #-} eigSH''@@ -535,7 +559,7 @@ >>> let (l, v) = eigSH a  >>> l-fromList [11.344814282762075,0.17091518882717918,-0.5157294715892575]+[11.344814282762075,0.17091518882717918,-0.5157294715892575]  >>> disp 3 $ v <> diag l <> tr v 3x3

@@ -415,6 +415,93 @@     OK } +//////////////////// generalized real eigensystem ////////////++int dggev_(char *jobvl, char *jobvr, integer *n,+    doublereal *a, integer *lda, doublereal *b, integer *ldb,+    doublereal *alphar, doublereal *alphai, doublereal *beta,+    doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr,+    doublereal *work,+	integer *lwork, integer *info);++int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) {+    integer n = ar;+    REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);+    REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);+    char jobvl = vlp==NULL?'N':'V';+    REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);+    char jobvr = vrp==NULL?'N':'V';+    DEBUGMSG("eig_l_G");+    integer lwork = -1;+    integer res;+    // ask for optimal lwork+    double ans;+    dggev_  (&jobvl,&jobvr,+             &n,+             ap,&n,bp,&n,+             (double*)alphap, (double*)alphap+n, betap,+             vlp, &n, vrp, &n,+             &ans, &lwork,+             &res);+    lwork = ceil(ans);+    double * work = (double*)malloc(lwork*sizeof(double));+    CHECK(!work,MEM);+    dggev_  (&jobvl,&jobvr,+             &n,+             ap,&n,bp,&n,+             (double*)alphap, (double*)alphap+n, betap,+             vlp, &n, vrp, &n,+             work, &lwork,+             &res);+    CHECK(res,res);+    free(work);+    OK+}++//////////////////// generalized complex eigensystem ////////////++int zggev_(char *jobvl, char *jobvr, integer *n,+    doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,+    doublecomplex *alphar, doublecomplex *beta,+    doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr,+    doublecomplex *work, integer *lwork,+    doublereal *rwork, integer *info);++int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) {+    integer n = ar;+    REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);+    REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);+    char jobvl = vlp==NULL?'N':'V';+    REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);+    char jobvr = vrp==NULL?'N':'V';+    DEBUGMSG("eig_l_GC");+    double *rwork = (double*) malloc(8*n*sizeof(double));+    CHECK(!rwork,MEM);+    integer lwork = -1;+    integer res;+    // ask for optimal lwork+    doublecomplex ans;+    zggev_  (&jobvl,&jobvr,+             &n,+             ap,&n,bp,&n,+             alphap, betap,+             vlp, &n, vrp, &n,+             &ans, &lwork,+             rwork, &res);+    lwork = ceil(ans.r);+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));+    CHECK(!work,MEM);+    zggev_  (&jobvl,&jobvr,+             &n,+             ap,&n,bp,&n,+             alphap, betap,+             vlp, &n, vrp, &n,+             work, &lwork,+             rwork, &res);+    CHECK(res,res);+    free(work);+    OK+}  //////////////////// symmetric real eigensystem //////////// 

@@ -40,10 +40,12 @@ {- | Creates a real vector containing a range of values:  >>> linspace 5 (-3,7::Double)-fromList [-3.0,-0.5,2.0,4.5,7.0]@+[-3.0,-0.5,2.0,4.5,7.0]+it :: Vector Double ->>> linspace 5 (8,2+i) :: Vector (Complex Double)-fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]+>>> linspace 5 (8,3:+2) :: Vector (Complex Double)+[8.0 :+ 0.0,6.75 :+ 0.5,5.5 :+ 1.0,4.25 :+ 1.5,3.0 :+ 2.0]+it :: Vector (Complex Double)  Logarithmic spacing can be defined as follows: @@ -87,7 +89,8 @@ >>> let v = vector [10,20,30]  >>> m #> v-fromList [140.0,320.0]+[140.0,320.0]+it :: Vector Numeric.LinearAlgebra.Data.R  -} infixr 8 #>@@ -136,10 +139,12 @@  >>> let x = a <\> v >>> x-fromList [3.0799999999999996,5.159999999999999]+[3.0799999999999996,5.159999999999999]+it :: Vector Numeric.LinearAlgebra.Data.R  >>> a #> x-fromList [13.399999999999999,26.799999999999997,1.0]+[13.399999999999999,26.799999999999997,0.9999999999999991]+it :: Vector Numeric.LinearAlgebra.Data.R  It also admits multiple right-hand sides stored as columns in a matrix. @@ -167,7 +172,8 @@   where     -- |     -- >>> build 5 (**2) :: Vector Double-    -- fromList [0.0,1.0,4.0,9.0,16.0]+    -- [0.0,1.0,4.0,9.0,16.0]+    -- it :: Vector Double     --     -- Hilbert matrix of order N:     --@@ -204,12 +210,11 @@  {- | Compute mean vector and covariance matrix of the rows of a matrix. ->>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])-(fromList [4.010341078059521,5.0197204699640405],-(2><2)- [     1.9862461923890056, -1.0127225830525157e-2- , -1.0127225830525157e-2,     3.0373954915729318 ])-+>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (trustSym $ diagl [2,3])+([3.9933155655086696,5.061409102770331],Herm (2><2)+ [    1.9963242906624408, -4.227815571404954e-2+ , -4.227815571404954e-2,    3.2003833097832857 ])+it :: (Vector Double, Herm Double) -} meanCov :: Matrix Double -> (Vector Double, Herm Double) meanCov x = (med,cov) where

@@ -41,7 +41,8 @@ {- ^ correlation  >>> corr (fromList[1,2,3]) (fromList [1..10])-fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]+[14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]+it :: (Enum t, Product t, Container Vector t) => Vector t  -} corr ker v@@ -54,7 +55,8 @@ {- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product)  >>> conv (fromList[1,1]) (fromList [-1,1])-fromList [-1.0,0.0,1.0]+[-1.0,0.0,1.0]+it :: (Product t, Container Vector t) => Vector t  -} conv ker v

@@ -18,6 +18,8 @@  module Internal.LAPACK where +import Data.Bifunctor (first)+ import Internal.Devel import Internal.Vector import Internal.Matrix hiding ((#), (#!))@@ -234,7 +236,9 @@ -----------------------------------------------------------------------------  foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok+foreign import ccall unsafe "eig_l_G" dggev :: R ::> R ::> C :> R :> R ::> R ::> Ok foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok+foreign import ccall unsafe "eig_l_GC" zggev :: C ::> C ::> C :> C :> C ::> C ::> Ok foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok @@ -298,13 +302,69 @@     | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) fixeig _ _ = error "fixeig with impossible inputs" +-- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0.+-- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1)+-- can be different. Therefore old 'fixeig' would fail for 'eigG'.+fixeigG  []  _  = []+fixeigG [_] [v] = [comp' v]+fixeigG ((_:+ai1) : an : as) (v1:v2:vs)+    | abs ai1 > 1e-13 = toComplex' (v1, v2) : toComplex' (v1, mapVector negate v2) : fixeigG as vs+    | otherwise = comp' v1 : fixeigG (an:as) (v2:vs)+fixeigG _ _ = error "fixeigG with impossible inputs"  -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. -- The eigenvalues are not sorted. eigOnlyR :: Matrix Double -> Vector (Complex Double) eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" +-- | Generalized eigenvalues and right eigenvectors of a pair of real matrices, using LAPACK's /dggev/.+-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.+eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double))+eigG a b = (alpha', beta, v'')+  where+    (alpha, beta, v) = eigGaux dggev a b "eigG"+    alpha' = fixeig1 alpha+    v' = toRows $ trans v+    v'' = fromColumns $ fixeigG (toList alpha') v' +eigGaux f ma mb st = unsafePerformIO $ do+    a <- copy ColumnMajor ma+    b <- copy ColumnMajor mb+    alpha <- createVector r+    beta <- createVector r+    vr <- createMatrix ColumnMajor r r++    (a # b # alpha # beta #! vr) g #| st++    return (alpha, beta, vr)+  where+    r = rows ma+    g ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta = f ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta 0 0 0 0 nullPtr ++eigGOnlyAux f ma mb st = unsafePerformIO $ do+    a <- copy ColumnMajor ma+    b <- copy ColumnMajor mb+    alpha <- createVector r+    beta <- createVector r++    (a # b # alpha #! beta) g #| st++    return (alpha, beta)+  where+    r = rows ma+    g ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta = f ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta 0 0 0 0 nullPtr 0 0 0 0 nullPtr++-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/.+-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.+eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))+eigGC a b = eigGaux zggev a b "eigGC"++eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)+eigOnlyG a b = first fixeig1 $ eigGOnlyAux dggev a b "eigOnlyG"++eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))+eigOnlyGC a b = eigGOnlyAux zggev a b "eigOnlyGC"+ -----------------------------------------------------------------------------  eigSHAux f st m = unsafePerformIO $ do@@ -446,8 +506,10 @@     | ndl  == nd - 1 &&       ndu  == nd - 1 &&       nd   == r = unsafePerformIO . g $ do+        dl' <- head . toRows <$> copy ColumnMajor (fromRows [dl])+        du' <- head . toRows <$> copy ColumnMajor (fromRows [du])         s <- copy ColumnMajor b-        (dl # d # du #! s) f #| st+        (dl' # d # du' #! s) f #| st         return s     | otherwise = error $ st ++ " of nonsquare matrix"   where

@@ -160,7 +160,8 @@ {- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.  >>> flatten (ident 3)-fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]+[1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]+it :: (Num t, Element t) => Vector t  -} flatten :: Element t => Matrix t -> Vector t

@@ -172,8 +172,10 @@ {- | general matrix - vector product  >>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)]+m :: GMatrix >>> m !#> vector [1..2000]-fromList [1000.0,4000.0]+[1000.0,4000.0]+it :: Vector Double  -} infixr 8 !#>

@@ -324,34 +324,34 @@ instance KnownNat n => Show (R n)   where     show s@(R (Dim v))-      | singleV v = "("++show (v!0)++" :: R "++show d++")"-      | otherwise   = "(vector"++ drop 8 (show v)++" :: R "++show d++")"+      | singleV v = "(" ++ show (v!0) ++ " :: R " ++ show d ++ ")"+      | otherwise = "(vector " ++ show v ++ " :: R " ++ show d ++")"       where         d = size s  instance KnownNat n => Show (C n)   where     show s@(C (Dim v))-      | singleV v = "("++show (v!0)++" :: C "++show d++")"-      | otherwise   = "(vector"++ drop 8 (show v)++" :: C "++show d++")"+      | singleV v = "(" ++ show (v!0) ++ " :: C " ++ show d ++ ")"+      | otherwise = "(vector " ++ show v ++ " :: C " ++ show d ++")"       where         d = size s  instance (KnownNat m, KnownNat n) => Show (L m n)   where-    show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (drop 9 $ show y) m' n'+    show (isDiag -> Just (z,y,(m',n'))) = printf "(diag %s %s :: L %d %d)" (show z) (show y) m' n'     show s@(L (Dim (Dim x)))        | singleM x = printf "(%s :: L %d %d)" (show (x `atIndex` (0,0))) m' n'-       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: L "++show m'++" "++show n'++")"+       | otherwise = "(matrix" ++ dropWhile (/='\n') (show x) ++ " :: L " ++ show m' ++ " " ++ show n' ++ ")"       where         (m',n') = size s  instance (KnownNat m, KnownNat n) => Show (M m n)   where-    show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (drop 9 $ show y) m' n'+    show (isDiagC -> Just (z,y,(m',n'))) = printf "(diag %s %s :: M %d %d)" (show z) (show y) m' n'     show s@(M (Dim (Dim x)))        | singleM x = printf "(%s :: M %d %d)" (show (x `atIndex` (0,0))) m' n'-       | otherwise = "(matrix"++ dropWhile (/='\n') (show x)++" :: M "++show m'++" "++show n'++")"+       | otherwise = "(matrix" ++ dropWhile (/='\n') (show x) ++ " :: M " ++ show m' ++ " " ++ show n' ++ ")"       where         (m',n') = size s 

@@ -99,7 +99,8 @@ {- | Create a real vector.  >>> vector [1..5]-fromList [1.0,2.0,3.0,4.0,5.0]+[1.0,2.0,3.0,4.0,5.0]+it :: Vector R  -} vector :: [R] -> Vector R@@ -378,7 +379,8 @@ On a matrix it gets the k-th row as a vector:  >>> matrix 5 [1..15] ! 1-fromList [6.0,7.0,8.0,9.0,10.0]+[6.0,7.0,8.0,9.0,10.0]+it :: Vector Double  >>> matrix 5 [1..15] ! 1 ! 3 9.0

@@ -116,7 +116,8 @@      be used, for instance, with infinite lists.  >>> 5 |> [1..]-fromList [1.0,2.0,3.0,4.0,5.0]+[1.0,2.0,3.0,4.0,5.0]+it :: (Enum a, Num a, Foreign.Storable.Storable a) => Vector a  -} (|>) :: (Storable a) => Int -> [a] -> Vector a@@ -135,7 +136,8 @@ {- | takes a number of consecutive elements from a Vector  >>> subVector 2 3 (fromList [1..10])-fromList [3.0,4.0,5.0]+[3.0,4.0,5.0]+it :: (Enum t, Num t, Foreign.Storable.Storable t) => Vector t  -} subVector :: Storable t => Int       -- ^ index of the starting element@@ -169,7 +171,8 @@ {- | concatenate a list of vectors  >>> vjoin [fromList [1..5::Double], konst 1 3]-fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]+[1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]+it :: Vector Double  -} vjoin :: Storable t => [Vector t] -> Vector t@@ -191,7 +194,8 @@ {- | Extract consecutive subvectors of the given sizes.  >>> takesV [3,4] (linspace 10 (1,10::Double))-[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]+[[1.0,2.0,3.0],[4.0,5.0,6.0,7.0]]+it :: [Vector Double]  -} takesV :: Storable t => [Int] -> Vector t -> [Vector t]

@@ -417,7 +417,8 @@  -- | -- >>> range 5--- fromList [0,1,2,3,4]+-- [0,1,2,3,4]+-- it :: Vector I -- range :: Int -> Vector I range n = unsafePerformIO $ do

@@ -29,7 +29,8 @@     -- as the Hadamard product or the Schur product):     --     -- >>>  vector [1,2,3] * vector [3,0,-2]-    -- fromList [3.0,0.0,-6.0]+    -- [3.0,0.0,-6.0]+    -- it :: Vector R     --     -- >>> matrix 3 [1..9] * ident 3     -- (3><3)@@ -130,8 +131,8 @@     leftSV, rightSV,      -- * Eigendecomposition-    eig, eigSH,-    eigenvalues, eigenvaluesSH,+    eig, geig, eigSH,+    eigenvalues, geigenvalues, eigenvaluesSH,     geigSH,      -- * QR

@@ -92,7 +92,7 @@   infixl 4 &-(&) :: forall n . (KnownNat n, 1 <= n)+(&) :: forall n . KnownNat n     => R n -> ℝ -> R (n+1) u & x = u # (konst x :: R 1)