hmatrix 0.1.0.0 → 0.1.1.0

raw patch · 10 files changed

+192/−94 lines, 10 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• examples/experiments/Static.hs +7/−15
• examples/experiments/listlike.hs +0/−8
• examples/tests.hs +44/−9
• hmatrix.cabal +6/−6
• lib/Data/Packed/Internal/Matrix.hs +4/−38
• lib/Numeric/LinearAlgebra/Algorithms.hs +51/−7
• lib/Numeric/LinearAlgebra/Instances.hs +11/−0
• lib/Numeric/LinearAlgebra/LAPACK.hs +20/−1
• lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +46/−0
• lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h +3/−10

@@ -6,23 +6,15 @@ import Numeric.LinearAlgebra import Foreign import Language.Haskell.TH.Syntax-import Data.Packed.Internal(Vector(..),Matrix(..))  instance Lift Double where   lift x = return (LitE (RationalL (toRational x))) -instance Lift (Ptr Double) where-    lift p = [e| p |]--instance Lift (ForeignPtr Double) where-    lift p = [e| p |]--instance (Lift a, Storable a, Lift (Ptr a), Lift (ForeignPtr a)) => Lift (Vector a ) where-    lift (V n fp) = [e| V $(lift n) $(lift fp) |]+instance Lift (Vector a ) where+    lift v = [e| v |] -instance (Lift (Vector a)) => Lift (Matrix a) where-    lift (MC r c v) = [e| MC $(lift r) $(lift c) $(lift v) |]-    lift (MF r c v) = [e| MF $(lift r) $(lift c) $(lift v) |]+instance Lift (Matrix a) where+    lift m = [e| m |]  tdim :: Int -> ExpQ tdim 0 = [| Z |]@@ -52,7 +44,7 @@ createv :: Storable t => d -> Vector t -> SVec d t createv d v = SVec v ---vec'' v = [|createv ($(tdim (dim v))) v|]+vec'' v = [|createv ($(tdim (dim v))) v|]  vec' :: [Double] -> ExpQ vec' d = [| createl ($(tdim (length d))) d |]@@ -71,8 +63,8 @@ vec d = mat (length d) 1 d  -mat' :: Matrix Double -> ExpQ-mat' m = [| createm ($(tdim (rows m))) ($(tdim (cols m))) m |]+--mat' :: Matrix Double -> ExpQ+--mat' m = [| createm ($(tdim (rows m))) ($(tdim (cols m))) m |]  covec :: [Double] -> ExpQ covec d = mat 1 (length d) d

@@ -3,15 +3,7 @@ import qualified Data.ListLike as LL import Numeric.LinearAlgebra import Data.Monoid-import Data.Packed.Internal.Vector import Foreign--instance (Storable a) => Monoid (Vector a) where-    mempty = V { dim = 0, fptr = undefined, ptr = undefined }-    mappend a b = mconcat [a,b]-    mconcat = j . filter ((>0).dim)-        where j [] = mempty-              j l  = join l  instance Storable a => LL.FoldableLL (Vector a) a where     foldl f x v = foldl f x (toList v)

@@ -2,20 +2,20 @@  module Main where -import Data.Packed.Internal((>|<), multiply', multiplyG, MatrixOrder(..),debug,fmat) import Numeric.GSL hiding (sin,cos,exp,choose) import Numeric.LinearAlgebra-import Numeric.LinearAlgebra.Linear(Linear) import Numeric.LinearAlgebra.LAPACK-import Numeric.GSL.Matrix(svdg) import qualified Numeric.GSL.Matrix as GSL import Test.QuickCheck hiding (test) import Test.HUnit hiding ((~:),test) import System.Random(randomRs,mkStdGen) import System.Info-import Data.List(foldl1')+import Data.List(foldl1', transpose) import System(getArgs)+import Debug.Trace(trace) +debug x = trace (show x) x+ type RM = Matrix Double type CM = Matrix (Complex Double) @@ -187,8 +187,23 @@  ------------------------------------------------------- -detTest = det m == 26 && det mc == 38 :+ (-3)+feye n = flipud (ident n) :: Matrix Double ++luTest1 m = m |~| p <> l <> u+    where (l,u,p,_) = lu m++detTest1 = det m == 26+        && det mc == 38 :+ (-3)+        && det (feye 2) == -1++detTest2 m = s d1 |~| s d2+    where d1 = det m+          d2 = det' m * det q+          det' m = product $ toList $ takeDiag r+          (q,r) = qr m+          s x = fromList [x]+ invTest m = degenerate m || m <> inv m |~| ident (rows m)  pinvTest m =  m <> p <> m |~| m@@ -340,9 +355,21 @@ asFortran m = (rows m >|< cols m) $ toList (flatten $ trans  m) asC m = (rows m >< cols m) $ toList (flatten m) -mulC a b = multiply' RowMajor a b-mulF a b = multiply' ColumnMajor a b+mulC a b = a <> b+mulF a b = trans $ trans b <> trans a +-------------------------------------------------------------------------++multiplyG a b = reshape (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b)+    where multiplyL a b = [[dotL x y | y <- transpose b] | x <- a]+          dotL a b = sum (zipWith (*) a b)++r >|< c = f where+    f l | dim v == r*c = reshapeF r v+        | otherwise    = error "(>|<)"+        where v = fromList l+    reshapeF r = trans . reshape r+ ---------------------------------------------------------------------  rot :: Double -> Matrix Double@@ -382,6 +409,14 @@     quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == trans (mulF (trans m2) (trans m1 :: CM))     quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: RM)     quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == multiplyG m1 (m2 :: CM)+    putStrLn "--------- lu ---------"+    quickCheck (luTest1 :: RM->Bool)+    quickCheck (luTest1 :: CM->Bool)+    quickCheck (detTest2 . sqm  :: SqM Double -> Bool)+    quickCheck (detTest2 . sqm  :: SqM (Complex Double) -> Bool)+    runTestTT $ TestList+     [ test "det1" detTest1+     ]     putStrLn "--------- svd ---------"     quickCheck (svdTest svdR)     quickCheck (svdTest svdRdd)@@ -389,7 +424,7 @@     quickCheck (svdTest' svdR)     quickCheck (svdTest' svdRdd)     quickCheck (svdTest' svdC)-    quickCheck (svdTest' svdg)+    quickCheck (svdTest' GSL.svdg)     putStrLn "--------- eig ---------"     quickCheck (eigTest  . sqm :: SqM Double -> Bool)     quickCheck (eigTest  . sqm :: SqM (Complex Double) -> Bool)@@ -450,7 +485,6 @@      , exponentialTest      , integrateTest      , polySolveTest-     , test "det" detTest      ]  bigtests = do@@ -460,6 +494,7 @@      , test "eigH" $ eigTestSH bigmatc      , test "eigR" $ eigTest   bigmat      , test "eigC" $ eigTest   bigmatc+     , test "det"  $ det (feye 1000) == 1 && det (feye 1002) == -1      ]  main = do

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.1.0.0+Version:            0.1.1.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -37,11 +37,7 @@     Extensions:         ForeignFunctionInterface      hs-source-dirs:     lib-    Exposed-modules:    Data.Packed.Internal,-                        Data.Packed.Internal.Common,-                        Data.Packed.Internal.Vector-                        Data.Packed.Internal.Matrix,-                        Data.Packed,+    Exposed-modules:    Data.Packed,                         Data.Packed.Vector,                         Data.Packed.Matrix,                         Numeric.GSL.Vector,@@ -87,6 +83,10 @@                         Numeric.LinearAlgebra.Interface,                         Numeric.LinearAlgebra.Algorithms,                         Graphics.Plot+    other-modules:      Data.Packed.Internal,+                        Data.Packed.Internal.Common,+                        Data.Packed.Internal.Vector,+                        Data.Packed.Internal.Matrix     C-sources:          lib/Data/Packed/Internal/auxi.c,                         lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,                         lib/Numeric/GSL/gsl-aux.c

@@ -254,9 +254,9 @@         r2 = dim d `div` c2         noneed = r1 == 1 || c1 == 1 -foreign import ccall safe "auxi.h transR"+foreign import ccall unsafe "auxi.h transR"     ctransR :: TMM -- Double ::> Double ::> IO Int-foreign import ccall safe "auxi.h transC"+foreign import ccall unsafe "auxi.h transC"     ctransC :: TCMCM -- Complex Double ::> Complex Double ::> IO Int  ------------------------------------------------------------------@@ -277,24 +277,19 @@     return r  multiplyR = multiplyAux cmultiplyR-foreign import ccall safe "auxi.h multiplyR"+foreign import ccall unsafe "auxi.h multiplyR"     cmultiplyR :: Int -> Int -> Int -> Ptr Double                -> Int -> Int -> Int -> Ptr Double                -> Int -> Int -> Ptr Double                -> IO Int  multiplyC = multiplyAux cmultiplyC-foreign import ccall safe "auxi.h multiplyC"+foreign import ccall unsafe "auxi.h multiplyC"     cmultiplyC :: Int -> Int -> Int -> Ptr (Complex Double)                -> Int -> Int -> Int -> Ptr (Complex Double)                -> Int -> Int -> Ptr (Complex Double)                -> IO Int -multiply' :: (Element a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a-multiply' RowMajor a b    = multiplyD a b-multiply' ColumnMajor a b = trans $ multiplyD (trans b) (trans a)-- -- | matrix product multiply :: (Element a) => Matrix a -> Matrix a -> Matrix a multiply = multiplyD@@ -402,32 +397,3 @@     --free charname  -- TO DO: free the auxiliary CString     return res foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM------------------------------------------------------------------------------- Generic definitions--{--transL m = matrixFromVector RowMajor (rows m) $ transdata (cols m) (cdat m) (rows m)--subMatrixG (r0,c0) (rt,ct) x = matrixFromVector RowMajor ct $ fromList $ concat $ map (subList c0 ct) (subList r0 rt (toLists x))-    where subList s n = take n . drop s--diagG v = matrixFromVector RowMajor c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..c]]-    where c = dim v-          l = toList v-          delta i j | i==j      = 1-                    | otherwise = 0--}--transdataG c1 d _ = fromList . concat . transpose . partit c1 . toList $ d--dotL a b = sum (zipWith (*) a b)--multiplyG a b = matrixFromVector RowMajor (cols b) $ fromList $ concat $ multiplyL (toLists a) (toLists b)--multiplyL a b | ok = [[dotL x y | y <- transpose b] | x <- a]-              | otherwise = error "inconsistent dimensions in contraction "-    where ok = case common length a of-                   Nothing -> False-                   Just c  -> c == length b

@@ -37,6 +37,8 @@     hess, -- ** Schur     schur,+-- ** LU+    lu, -- * Matrix functions     expm,     sqrtm,@@ -52,11 +54,11 @@ -- * Util     haussholder,     unpackQR, unpackHess,-    Field(linearSolveSVD,lu,eigSH',cholSH)+    Field(linearSolveSVD,eigSH',cholSH) ) where  -import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)+import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj, (//)) import Data.Packed import qualified Numeric.GSL.Matrix as GSL import Numeric.GSL.Vector@@ -64,12 +66,13 @@ import Complex import Numeric.LinearAlgebra.Linear import Data.List(foldl1')+import Data.Array  -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. class (Normed (Matrix t), Linear Matrix t) => Field t where     -- | Singular value decomposition using lapack's dgesvd or zgesvd.     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)-    lu          :: Matrix t -> (Matrix t, Matrix t, [Int], t)+    luPacked    :: Matrix t -> (Matrix t, [Int])     -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv and zgesv.     --  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".     linearSolve :: Matrix t -> Matrix t -> Matrix t@@ -106,7 +109,7 @@  instance Field Double where     svd = svdR-    lu  = GSL.luR+    luPacked = luR     linearSolve = linearSolveR     linearSolveSVD = linearSolveSVDR Nothing     ctrans = trans@@ -119,7 +122,7 @@  instance Field (Complex Double) where     svd = svdC-    lu  = GSL.luC+    luPacked = luC     linearSolve = linearSolveC     linearSolveSVD = linearSolveSVDC Nothing     ctrans = conj . trans@@ -146,11 +149,20 @@  square m = rows m == cols m +-- | determinant of a square matrix, computed from the LU decomposition. det :: Field t => Matrix t -> t-det m | square m = s * (product $ toList $ takeDiag $ u)+det m | square m = s * (product $ toList $ takeDiag $ lu)       | otherwise = error "det of nonsquare matrix"-    where (_,u,_,s) = lu m+    where (lu,perm) = luPacked m+          s = signlp (rows m) perm +-- | LU factorization of a general matrix using lapack's dgetrf or zgetrf.+--+-- If @(l,u,p,s) = lu m@ then @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 :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)+lu = luFact . luPacked+ -- | Inverse of a square matrix using lapacks' dgesv and zgesv. inv :: Field t => Matrix t -> Matrix t inv m | square m = m `linearSolve` ident (rows m)@@ -457,3 +469,35 @@           (.*) = scale           (|+|) = add           (|-|) = sub++------------------------------------------------------------------++signlp r vals = foldl f 1 (zip [0..r-1] vals)+    where f s (a,b) | a /= b    = -s+                    | otherwise =  s++swap (arr,s) (a,b) | a /= b    = (arr // [(a, arr!b),(b,arr!a)],-s)+                   | otherwise = (arr,s)++fixPerm r vals = (fromColumns $ elems res, sign)+    where v = [0..r-1]+          s = toColumns (ident r)+          (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)++triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]+    where el i j = if j-i>=h then v else 1 - v++luFact (lu,perm) | r <= c    = (l ,u ,p, s)+                 | otherwise = (l',u',p, s)+  where+    r = rows lu+    c = cols lu+    tu = triang r c 0 1+    tl = triang r c 0 0+    l = takeColumns r (lu |*| tl) |+| diagRect (constant 1 r) r r+    u = lu |*| tu+    (p,s) = fixPerm r perm+    l' = (lu |*| tl) |+| diagRect (constant 1 c) r c+    u' = takeRows c (lu |*| tu)+    (|+|) = add+    (|*|) = mul

@@ -26,6 +26,8 @@ import Complex import Data.List(transpose,intersperse) import Foreign(Storable)+import Data.Monoid+import Data.Packed.Internal.Vector  ------------------------------------------------------------------ @@ -159,3 +161,12 @@     (**)  = liftMatrix2' (**)     sqrt  = liftMatrix sqrt     pi    = (1><1) [pi]++---------------------------------------------------------------++instance (Storable a) => Monoid (Vector a) where+    mempty = V { dim = 0, fptr = undefined }+    mappend a b = mconcat [a,b]+    mconcat = j . filter ((>0).dim)+        where j [] = mempty+              j l  = join l

@@ -19,6 +19,7 @@     linearSolveR, linearSolveC,     linearSolveLSR, linearSolveLSC,     linearSolveSVDR, linearSolveSVDC,+    luR, luC,     cholS, cholH,     qrR, qrC,     hessR, hessC,@@ -299,7 +300,7 @@         mn = min m n  ------------------------------------------------------------------------------------foreign import ccall safe "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM+foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM  -- | Wrapper for LAPACK's /dgees/, which computes a Schur factorization of a square real matrix.@@ -318,3 +319,21 @@   where n = rows a  -----------------------------------------------------------------------------------+foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM+foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM++-- | Wrapper for LAPACK's /dgetrf/, which computes a LU factorization of a general real matrix.+luR :: Matrix Double -> (Matrix Double, [Int])+luR = luAux dgetrf "luR" . fmat++-- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix.+luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])+luC = luAux zgetrf "luC" . fmat++luAux f st a = unsafePerformIO $ do+    lu <- createMatrix ColumnMajor n m+    piv <- createVector (min n m)+    app3 f mat a vec piv mat lu st+    return (lu, map (pred.round) (toList piv))+  where n = rows a+        m = cols a

@@ -768,3 +768,49 @@     OK     #endif }++//////////////////// LU factorization /////////////////////////++int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) {+    integer m = ar;+    integer n = ac;+    integer mn = MIN(m,n);+    REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);+    DEBUGMSG("lu_l_R");+    integer* auxipiv = (integer*)malloc(mn*sizeof(integer));+    memcpy(rp,ap,m*n*sizeof(double));+    integer res;+    dgetrf_ (&m,&n,rp,&m,auxipiv,&res);+    if(res>0) {+        res = 0; // fixme+    }+    CHECK(res,res);+    int k;+    for (k=0; k<mn; k++) {+        ipivp[k] = auxipiv[k];+    }+    free(auxipiv);+    OK+}++int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {+    integer m = ar;+    integer n = ac;+    integer mn = MIN(m,n);+    REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);+    DEBUGMSG("lu_l_C");+    integer* auxipiv = (integer*)malloc(mn*sizeof(integer));+    memcpy(rp,ap,m*n*sizeof(doublecomplex));+    integer res;+    zgetrf_ (&m,&n,(doublecomplex*)rp,&m,auxipiv,&res);+    if(res>0) {+        res = 0; // fixme+    }+    CHECK(res,res);+    int k;+    for (k=0; k<mn; k++) {+        ipivp[k] = auxipiv[k];+    }+    free(auxipiv);+    OK+}

@@ -14,41 +14,34 @@  int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v)); int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));- int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));  int eig_l_C(KCMAT(a),CMAT(u),CVEC(s),CMAT(v));- int eig_l_R(KDMAT(a),DMAT(u),CVEC(s),DMAT(v));  int eig_l_S(KDMAT(a),DVEC(s),DMAT(v));- int eig_l_H(KCMAT(a),DVEC(s),CMAT(v));  int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x));- int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x));  int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x));- int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x));  int linearSolveSVDR_l(double,KDMAT(a),KDMAT(b),DMAT(x));- int linearSolveSVDC_l(double,KCMAT(a),KCMAT(b),CMAT(x));  int chol_l_H(KCMAT(a),CMAT(r));- int chol_l_S(KDMAT(a),DMAT(r));  int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r));- int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r));  int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r));- int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r));  int schur_l_R(KDMAT(a), DMAT(u), DMAT(s));- int schur_l_C(KCMAT(a), CMAT(u), CMAT(s));++int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r));+int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r));