From d14515a4a50d5b5335f9c1525432b68ab67fa7c8 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 14 Sep 2007 18:23:20 +0000 Subject: more refactoring --- lib/Data/Packed/Internal/Matrix.hs | 7 +++++- lib/Data/Packed/Matrix.hs | 4 ++-- lib/GSL/Matrix.hs | 2 +- lib/GSLHaskell.hs | 44 +++++++++----------------------------- lib/LinearAlgebra.hs | 15 ++++++++++++- lib/LinearAlgebra/Linear.hs | 40 ++++++++++++++++++++++++++++++---- 6 files changed, 69 insertions(+), 43 deletions(-) (limited to 'lib') diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 89a162a..58b325c 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -348,7 +348,7 @@ foreign import ccall safe "aux.h constantC" @> constant 2 7 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ -} -constant :: Field a => a -> Int -> Vector a +constant :: Double -> Int -> Vector Double constant = constantD -------------------------------------------------------------------------- @@ -361,6 +361,11 @@ conj v = asComplex $ cdat $ reshape 2 (asReal v) `multiply` diag (fromList [1,-1 toComplex :: (Vector Double, Vector Double) -> Vector (Complex Double) toComplex (r,i) = asComplex $ cdat $ fromColumns [r,i] +-- | the inverse of 'toComplex' +fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double) +fromComplex z = (r,i) where + [r,i] = toColumns $ reshape 2 $ asReal z + -- | converts a real vector into a complex representation (with zero imaginary parts) comp :: Vector Double -> Vector (Complex Double) comp v = toComplex (v,constant 0 (dim v)) diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs index 01e8133..fc08ce4 100644 --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs @@ -75,12 +75,12 @@ diagRect s r c | r == c = diag s | r < c = trans $ diagRect s c r | r > c = joinVert [diag s , zeros (r-c,c)] - where zeros (r,c) = reshape c $ constant 0 (r*c) + where zeros (r,c) = reshape c $ constantD 0 (r*c) takeDiag :: (Field t) => Matrix t -> Vector t takeDiag m = fromList [cdat m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] -ident :: (Num t, Field t) => Int -> Matrix t +ident :: Int -> Matrix Double ident n = diag (constant 1 n) ------------------------------------------------------------ diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs index 15710df..ec8ceea 100644 --- a/lib/GSL/Matrix.hs +++ b/lib/GSL/Matrix.hs @@ -289,7 +289,7 @@ luC m = (l,u,p, fromIntegral s') where 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 + l = (triang r r 0 0 `mul` lu) `add` liftMatrix comp (ident r) add = liftMatrix2 $ vectorZipC Add mul = liftMatrix2 $ vectorZipC Mul diff --git a/lib/GSLHaskell.hs b/lib/GSLHaskell.hs index e65ff28..3158458 100644 --- a/lib/GSLHaskell.hs +++ b/lib/GSLHaskell.hs @@ -18,6 +18,7 @@ module GSLHaskell( module Data.Packed.Vector, module Data.Packed.Matrix, module LinearAlgebra.Algorithms, + module LinearAlgebra.Linear, module LAPACK, module GSL.Integration, module GSL.Differentiation, @@ -28,18 +29,15 @@ module GSLHaskell( module GSL.Special, module Graphics.Plot, module Complex, - Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->), GSLHaskell.constant, - fromArray2D, fromComplex, toComplex, GSLHaskell.pnorm, scale, outer + Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->), + fromArray2D, GSLHaskell.pnorm, ) where - -import LAPACK import GSL.Integration import GSL.Differentiation import GSL.Fourier import GSL.Polynomials import GSL.Minimization -import GSL.Matrix import Graphics.Plot import Complex import GSL.Special(setErrorHandlerOff, @@ -48,14 +46,16 @@ import GSL.Special(setErrorHandlerOff, bessel_J0_e, exp_e10_e, gamma) -import Data.Packed.Internal hiding (dsp) -import Data.Packed.Vector hiding (constant) +--import Data.Packed.Internal hiding (dsp,comp) +import Data.Packed.Vector import Data.Packed.Matrix import Data.Packed.Matrix hiding ((><)) import GSL.Vector -import LinearAlgebra.Linear import qualified LinearAlgebra.Algorithms +import LAPACK +import GSL.Matrix import LinearAlgebra.Algorithms hiding (pnorm) +import LinearAlgebra.Linear import Complex import Numeric(showGFloat) import Data.List(transpose,intersperse) @@ -69,7 +69,7 @@ adaptScalar f1 f2 f3 x y | otherwise = f2 x y liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t -liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (cdat m1) (cdat m2)) +liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2)) | otherwise = error "nonconformant matrices in liftMatrix2'" compat' :: Matrix a -> Matrix b -> Bool @@ -89,7 +89,7 @@ instance (Linear Vector a) => Num (Vector a) where fromInteger = fromList . return . fromInteger instance (Eq a, Field a) => Eq (Matrix a) where - a == b = rows a == rows b && cols a == cols b && cdat a == cdat b && fdat a == fdat b + a == b = cols a == cols b && flatten a == flatten b instance (Field a, Linear Vector a) => Num (Matrix a) where (+) = liftMatrix2' (+) @@ -377,9 +377,6 @@ size = dim gmap :: (Storable a, Storable b) => (a->b) -> Vector a -> Vector b gmap f v = liftVector f v -constant :: Double -> Int -> Vector Double -constant = constantR - -- shows a Double with n digits after the decimal point shf :: (RealFloat a) => Int -> a -> String shf dec n | abs n < 1e-10 = "0." @@ -473,27 +470,6 @@ fromArray2D m = (r> Vector (Complex Double) -toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] - --- | extracts the real and imaginary parts of a complex vector -fromComplexV :: Vector (Complex Double) -> (Vector Double, Vector Double) -fromComplexV m = (a,b) where [a,b] = toColumns $ reshape 2 $ asReal m - --- | creates a complex matrix from matrices with real and imaginary parts -toComplexM :: (Matrix Double, Matrix Double) -> Matrix (Complex Double) -toComplexM (r,i) = reshape (cols r) $ asComplex $ flatten $ fromColumns [flatten r, flatten i] - --- | extracts the real and imaginary parts of a complex matrix -fromComplexM :: Matrix (Complex Double) -> (Matrix Double, Matrix Double) -fromComplexM m = (reshape c a, reshape c b) - where c = cols m - [a,b] = toColumns $ reshape 2 $ asReal $ flatten m - -fromComplex :: Matrix (Complex Double) -> (Matrix Double, Matrix Double) -fromComplex = fromComplexM - pnorm :: (Normed t1, Num t) => t -> t1 -> Double pnorm 0 = LinearAlgebra.Algorithms.pnorm Infinity diff --git a/lib/LinearAlgebra.hs b/lib/LinearAlgebra.hs index b0c8b9d..3b56fc4 100644 --- a/lib/LinearAlgebra.hs +++ b/lib/LinearAlgebra.hs @@ -13,6 +13,19 @@ Some linear algebra algorithms, implemented by means of BLAS, LAPACK or GSL. -} ----------------------------------------------------------------------------- module LinearAlgebra ( - + module Data.Packed.Vector, + module Data.Packed.Matrix, + module LinearAlgebra.Linear, + module LAPACK, + module GSL.Matrix, + module LinearAlgebra.Algorithms, + module Complex ) where +import LinearAlgebra.Linear +import LinearAlgebra.Algorithms +import LAPACK +import GSL.Matrix +import Data.Packed.Matrix +import Data.Packed.Vector +import Complex \ No newline at end of file diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs index 53e011f..a2071ed 100644 --- a/lib/LinearAlgebra/Linear.hs +++ b/lib/LinearAlgebra/Linear.hs @@ -15,8 +15,6 @@ Portability : uses ffi module LinearAlgebra.Linear ( Linear(..), - toComplex, comp, - conj, multiply, dot, outer ) where @@ -33,6 +31,10 @@ class (Field e) => Linear c e where add :: c e -> c e -> c e sub :: c e -> c e -> c e mul :: c e -> c e -> c e + toComplex :: RealFloat e => (c e, c e) -> c (Complex e) + fromComplex :: RealFloat e => c (Complex e) -> (c e, c e) + comp :: RealFloat e => c e -> c (Complex e) + conj :: RealFloat e => c (Complex e) -> c (Complex e) instance Linear Vector Double where scale = vectorMapValR Scale @@ -40,6 +42,10 @@ instance Linear Vector Double where add = vectorZipR Add sub = vectorZipR Sub mul = vectorZipR Mul + toComplex = Data.Packed.Internal.toComplex + fromComplex = Data.Packed.Internal.fromComplex + comp = Data.Packed.Internal.comp + conj = Data.Packed.Internal.conj instance Linear Vector (Complex Double) where scale = vectorMapValC Scale @@ -47,6 +53,34 @@ instance Linear Vector (Complex Double) where add = vectorZipC Add sub = vectorZipC Sub mul = vectorZipC Mul + toComplex = undefined -- can't match + fromComplex = undefined + comp = undefined + conj = undefined + +instance Linear Matrix Double where + scale x = liftMatrix (scale x) + addConstant x = liftMatrix (addConstant x) + add = liftMatrix2 add + sub = liftMatrix2 sub + mul = liftMatrix2 mul + toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex + fromComplex z = (reshape c r, reshape c i) + where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z) + c = cols z + comp = liftMatrix Data.Packed.Internal.comp + conj = liftMatrix Data.Packed.Internal.conj + +instance Linear Matrix (Complex Double) where + scale x = liftMatrix (scale x) + addConstant x = liftMatrix (addConstant x) + add = liftMatrix2 add + sub = liftMatrix2 sub + mul = liftMatrix2 mul + toComplex = undefined + fromComplex = undefined + comp = undefined + conj = undefined -------------------------------------------------- @@ -58,8 +92,6 @@ dot u v = dat (multiply r c) `at` 0 c = asColumn v - - {- | Outer product of two vectors. @\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] -- cgit v1.2.3