diff --git a/examples/experiments/Static.hs b/examples/experiments/Static.hs --- a/examples/experiments/Static.hs +++ b/examples/experiments/Static.hs @@ -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 diff --git a/examples/experiments/listlike.hs b/examples/experiments/listlike.hs --- a/examples/experiments/listlike.hs +++ b/examples/experiments/listlike.hs @@ -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) diff --git a/examples/tests.hs b/examples/tests.hs --- a/examples/tests.hs +++ b/examples/tests.hs @@ -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 diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -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 diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -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 diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -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 diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Instances.hs @@ -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 diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs @@ -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 diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c @@ -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=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