From 1871acb835b4fc164bcff3f6e7467884b87fbd0f Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 25 Jun 2007 07:32:56 +0000 Subject: l.a. algorithms, etc. --- lib/GSL/Matrix.hs | 311 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ lib/GSL/Vector.hs | 7 +- 2 files changed, 317 insertions(+), 1 deletion(-) create mode 100644 lib/GSL/Matrix.hs (limited to 'lib/GSL') diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs new file mode 100644 index 0000000..919c2d9 --- /dev/null +++ b/lib/GSL/Matrix.hs @@ -0,0 +1,311 @@ +----------------------------------------------------------------------------- +-- | +-- Module : GSL.Matrix +-- Copyright : (c) Alberto Ruiz 2007 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable (uses FFI) +-- +-- A few linear algebra computations based on the GSL (). +-- +----------------------------------------------------------------------------- + +module GSL.Matrix( + eigSg, eigHg, + svdg, + qr, + chol, + luSolveR, luSolveC, + luR, luC, + fromFile +) where + +import Data.Packed.Internal +import Data.Packed.Matrix(fromLists,ident,takeDiag) +import GSL.Vector +import Foreign +import Foreign.C.Types +import Complex +import Foreign.C.String + +{- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/. + +> > let (l,v) = eigS $ 'fromLists' [[1,2],[2,1]] +> > l +> 3.000 -1.000 +> +> > v +> 0.707 -0.707 +> 0.707 0.707 +> +> > v <> diag l <> trans v +> 1.000 2.000 +> 2.000 1.000 + +-} +eigSg :: Matrix Double -> (Vector Double, Matrix Double) +eigSg (m@M {rows = r}) + | r == 1 = (fromList [cdat m `at` 0], singleton 1) + | otherwise = unsafePerformIO $ do + l <- createVector r + v <- createMatrix RowMajor r r + c_eigS // mat cdat m // vec l // mat dat v // check "eigSg" [cdat m] + return (l,v) +foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM + +------------------------------------------------------------------ + + + +{- | eigendecomposition of a complex hermitian matrix using /gsl_eigen_hermv/ + +> > let (l,v) = eigH $ 'fromLists' [[1,2+i],[2-i,3]] +> +> > l +> 4.449 -0.449 +> +> > v +> -0.544 0.839 +> (-0.751,0.375) (-0.487,0.243) +> +> > v <> diag l <> (conjTrans) v +> 1.000 (2.000,1.000) +> (2.000,-1.000) 3.000 + +-} +eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double)) +eigHg (m@M {rows = r}) + | r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1) + | otherwise = unsafePerformIO $ do + l <- createVector r + v <- createMatrix RowMajor r r + c_eigH // mat cdat m // vec l // mat dat v // check "eigHg" [cdat m] + return (l,v) +foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM + + +{- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/: + +@\> let (u,s,v) = svdg $ 'fromLists' [[1,2,3],[-4,1,7]] +\ +\> u +0.310 -0.951 +0.951 0.310 +\ +\> s +8.497 2.792 +\ +\> v +-0.411 -0.785 + 0.185 -0.570 + 0.893 -0.243 +\ +\> u \<\> 'diag' s \<\> 'trans' v + 1. 2. 3. +-4. 1. 7.@ + +-} +svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +svdg x@M {rows = r, cols = c} = if r>=c + then svd' x + else (v, s, u) where (u,s,v) = svd' (trans x) + +svd' x@M {rows = r, cols = c} = unsafePerformIO $ do + u <- createMatrix RowMajor r c + s <- createVector c + v <- createMatrix RowMajor c c + c_svd // mat cdat x // mat dat u // vec s // mat dat v // check "svdg" [cdat x] + return (u,s,v) +foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM + +{- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/. + +@\> let (q,r) = qr $ 'fromLists' [[1,3,5,7],[2,0,-2,4]] +\ +\> q +-0.447 -0.894 +-0.894 0.447 +\ +\> r +-2.236 -1.342 -0.447 -6.708 + 0. -2.683 -5.367 -4.472 +\ +\> q \<\> r +1.000 3.000 5.000 7.000 +2.000 0. -2.000 4.000@ + +-} +qr :: Matrix Double -> (Matrix Double, Matrix Double) +qr x@M {rows = r, cols = c} = unsafePerformIO $ do + q <- createMatrix RowMajor r r + rot <- createMatrix RowMajor r c + c_qr // mat cdat x // mat dat q // mat dat rot // check "qr" [cdat x] + return (q,rot) +foreign import ccall "gsl-aux.h QR" c_qr :: TMMM + +{- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/. + +@\> let c = chol $ 'fromLists' [[5,4],[4,5]] +\ +\> c +2.236 0. +1.789 1.342 +\ +\> c \<\> 'trans' c +5.000 4.000 +4.000 5.000@ + +-} +chol :: Matrix Double -> Matrix Double +--chol x@(M r _ p) = createM [p] "chol" r r $ m c_chol x +chol x@M {rows = r} = unsafePerformIO $ do + res <- createMatrix RowMajor r r + c_chol // mat cdat x // mat dat res // check "chol" [cdat x] + return res +foreign import ccall "gsl-aux.h chol" c_chol :: TMM + +-------------------------------------------------------- + +{- -| efficient multiplication by the inverse of a matrix (for real matrices) +-} +luSolveR :: Matrix Double -> Matrix Double -> Matrix Double +luSolveR a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c}) + | n1==n2 && n1==r = unsafePerformIO $ do + s <- createMatrix RowMajor r c + c_luSolveR // mat cdat a // mat cdat b // mat dat s // check "luSolveR" [cdat a, cdat b] + return s + | otherwise = error "luSolveR of nonsquare matrix" + +foreign import ccall "gsl-aux.h luSolveR" c_luSolveR :: TMMM + +{- -| efficient multiplication by the inverse of a matrix (for complex matrices). +-} +luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +luSolveC a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c}) + | n1==n2 && n1==r = unsafePerformIO $ do + s <- createMatrix RowMajor r c + c_luSolveC // mat cdat a // mat cdat b // mat dat s // check "luSolveC" [cdat a, cdat b] + return s + | otherwise = error "luSolveC of nonsquare matrix" + +foreign import ccall "gsl-aux.h luSolveC" c_luSolveC :: TCMCMCM + +{- | lu decomposition of real matrix (packed as a vector including l, u, the permutation and sign) +-} +luRaux :: Matrix Double -> Vector Double +luRaux x@M {rows = r, cols = c} = unsafePerformIO $ do + res <- createVector (r*r+r+1) + c_luRaux // mat cdat x // vec res // check "luRaux" [cdat x] + return res +foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV + +{- | lu decomposition of complex matrix (packed as a vector including l, u, the permutation and sign) +-} +luCaux :: Matrix (Complex Double) -> Vector (Complex Double) +luCaux x@M {rows = r, cols = c} = unsafePerformIO $ do + res <- createVector (r*r+r+1) + c_luCaux // mat cdat x // vec res // check "luCaux" [cdat x] + return res +foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV + +{- | The LU decomposition of a square matrix. Is based on /gsl_linalg_LU_decomp/ and /gsl_linalg_complex_LU_decomp/ as described in . + +@\> let m = 'fromLists' [[1,2,-3],[2+3*i,-7,0],[1,-i,2*i]] +\> let (l,u,p,s) = luR m@ + +L is the lower triangular: + +@\> l + 1. 0. 0. +0.154-0.231i 1. 0. +0.154-0.231i 0.624-0.522i 1.@ + +U is the upper triangular: + +@\> u +2.+3.i -7. 0. + 0. 3.077-1.615i -3. + 0. 0. 1.873+0.433i@ + +p is a permutation: + +@\> p +[1,0,2]@ + +L \* U obtains a permuted version of the original matrix: + +@\> 'extractRows' p m + 2.+3.i -7. 0. + 1. 2. -3. + 1. -1.i 2.i +\ +\> l \<\> u + 2.+3.i -7. 0. + 1. 2. -3. + 1. -1.i 2.i@ + +s is the sign of the permutation, required to obtain sign of the determinant: + +@\> s * product ('toList' $ 'takeDiag' u) +(-18.0) :+ (-16.000000000000004) +\> 'LinearAlgebra.Algorithms.det' m +(-18.0) :+ (-16.000000000000004)@ + + -} +luR :: Matrix Double -> (Matrix Double, Matrix Double, [Int], Double) +luR m = (l,u,p, fromIntegral s') where + r = rows m + v = luRaux m + lu = reshape r $ subVector 0 (r*r) v + s':p = map round . toList . subVector (r*r) (r+1) $ v + u = triang r r 0 1`mul` lu + l = (triang r r 0 0 `mul` lu) `add` ident r + add = liftMatrix2 $ vectorZipR Add + mul = liftMatrix2 $ vectorZipR Mul + +-- | Complex version of 'luR'. +luC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double), [Int], Complex Double) +luC m = (l,u,p, fromIntegral s') where + r = rows m + v = luCaux m + lu = reshape r $ subVector 0 (r*r) v + s':p = map (round.realPart) . toList . subVector (r*r) (r+1) $ v + u = triang r r 0 1 `mul` lu + l = (triang r r 0 0 `mul` lu) `add` ident r + add = liftMatrix2 $ vectorZipC Add + mul = liftMatrix2 $ vectorZipC Mul + +extract l is = [l!!i |i<-is] + +{- auxiliary function to get triangular matrices +-} +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 + +{- | rearranges the rows of a matrix according to the order given in a list of integers. + +> > extractRows [3,3,0,1] (ident 4) +> 0. 0. 0. 1. +> 0. 0. 0. 1. +> 1. 0. 0. 0. +> 0. 1. 0. 0. + +-} +extractRows :: Field t => [Int] -> Matrix t -> Matrix t +extractRows l m = fromRows $ extract (toRows $ m) l + +-------------------------------------------------------------- + +-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance). +fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) +fromFile filename (r,c) = do + charname <- newCString filename + res <- createMatrix RowMajor r c + c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" [] + --free charname -- TO DO: free the auxiliary CString + return res +foreign import ccall "gsl-aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM + +--------------------------------------------------------------------------- diff --git a/lib/GSL/Vector.hs b/lib/GSL/Vector.hs index a772b34..a074254 100644 --- a/lib/GSL/Vector.hs +++ b/lib/GSL/Vector.hs @@ -18,7 +18,7 @@ module GSL.Vector ( FunCodeV(..), vectorMapR, vectorMapC, FunCodeSV(..), vectorMapValR, vectorMapValC, FunCodeVV(..), vectorZipR, vectorZipC, - scale, addConstant, add, mul, + scale, addConstant, add, sub, mul, ) where import Data.Packed.Internal @@ -84,6 +84,11 @@ add u v | isReal baseOf v = scast $ vectorZipR Add (scast u) (scast v) | isComp baseOf v = scast $ vectorZipC Add (scast u) (scast v) | otherwise = fromList $ zipWith (+) (toList u) (toList v) +sub :: (Num a, Field a) => Vector a -> Vector a -> Vector a +sub u v | isReal baseOf v = scast $ vectorZipR Sub (scast u) (scast v) + | isComp baseOf v = scast $ vectorZipC Sub (scast u) (scast v) + | otherwise = fromList $ zipWith (-) (toList u) (toList v) + mul :: (Num a, Field a) => Vector a -> Vector a -> Vector a mul u v | isReal baseOf v = scast $ vectorZipR Mul (scast u) (scast v) | isComp baseOf v = scast $ vectorZipC Mul (scast u) (scast v) -- cgit v1.2.3