module Data.Packed.Internal.Matrix where
import Data.Packed.Internal.Common
import Data.Packed.Internal.Vector
import Foreign hiding (xor)
import Complex
import Foreign.C.String
import Foreign.C.Types
data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
data Matrix t = MC { rows :: !Int
, cols :: !Int
, cdat :: !(Vector t) }
| MF { rows :: !Int
, cols :: !Int
, fdat :: !(Vector t) }
trans :: Matrix t -> Matrix t
trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d }
trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d }
cmat m@MC{} = m
cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c}
fmat m@MF{} = m
fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r}
mat = withMatrix
withMatrix MC {rows = r, cols = c, cdat = d } f =
withForeignPtr (fptr d) $ \p -> do
let m g = do
g (fi r) (fi c) p
f m
withMatrix MF {rows = r, cols = c, fdat = d } f =
withForeignPtr (fptr d) $ \p -> do
let m g = do
g (fi r) (fi c) p
f m
flatten :: Element t => Matrix t -> Vector t
flatten = cdat . cmat
type Mt t s = Int -> Int -> Ptr t -> s
toLists :: (Element t) => Matrix t -> [[t]]
toLists m = partit (cols m) . toList . flatten $ m
fromRows :: Element t => [Vector t] -> Matrix t
fromRows vs = case common dim vs of
Nothing -> error "fromRows applied to [] or to vectors with different sizes"
Just c -> reshape c (join vs)
toRows :: Element t => Matrix t -> [Vector t]
toRows m = toRows' 0 where
v = flatten $ m
r = rows m
c = cols m
toRows' k | k == r*c = []
| otherwise = subVector k c v : toRows' (k+c)
fromColumns :: Element t => [Vector t] -> Matrix t
fromColumns m = trans . fromRows $ m
toColumns :: Element t => Matrix t -> [Vector t]
toColumns m = toRows . trans $ m
(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
infixl 9 @@>
MC {rows = r, cols = c, cdat = v} @@> (i,j)
| safe = if i<0 || i>=r || j<0 || j>=c
then error "matrix indexing out of range"
else v `at` (i*c+j)
| otherwise = v `at` (i*c+j)
MF {rows = r, cols = c, fdat = v} @@> (i,j)
| safe = if i<0 || i>=r || j<0 || j>=c
then error "matrix indexing out of range"
else v `at` (j*r+i)
| otherwise = v `at` (j*r+i)
matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v }
where (d,m) = dim v `divMod` c
r | m==0 = d
| otherwise = error "matrixFromVector"
matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v }
where (d,m) = dim v `divMod` c
r | m==0 = d
| otherwise = error "matrixFromVector"
createMatrix order r c = do
p <- createVector (r*c)
return (matrixFromVector order c p)
reshape :: Element t => Int -> Vector t -> Matrix t
reshape c v = matrixFromVector RowMajor c v
singleton x = reshape 1 (fromList [x])
liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor c (f d)
liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)
liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
liftMatrix2 f m1 m2
| not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
| otherwise = case m1 of
MC {} -> matrixFromVector RowMajor (cols m1) (f (cdat m1) (flatten m2))
MF {} -> matrixFromVector ColumnMajor (cols m1) (f (fdat m1) ((fdat.fmat) m2))
compat :: Matrix a -> Matrix b -> Bool
compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
class (Storable a, Floating a) => Element a where
constantD :: a -> Int -> Vector a
transdata :: Int -> Vector a -> Int -> Vector a
subMatrixD :: (Int,Int)
-> (Int,Int)
-> Matrix a -> Matrix a
diagD :: Vector a -> Matrix a
instance Element Double where
constantD = constantR
transdata = transdataR
subMatrixD = subMatrixR
diagD = diagR
instance Element (Complex Double) where
constantD = constantC
transdata = transdataC
subMatrixD = subMatrixC
diagD = diagC
(>|<) :: (Element a) => Int -> Int -> [a] -> Matrix a
r >|< c = f where
f l | dim v == r*c = matrixFromVector ColumnMajor c v
| otherwise = error $ "inconsistent list size = "
++show (dim v) ++" in ("++show r++"><"++show c++")"
where v = fromList l
transdataR :: Int -> Vector Double -> Int -> Vector Double
transdataR = transdataAux ctransR
transdataC :: Int -> Vector (Complex Double) -> Int -> Vector (Complex Double)
transdataC = transdataAux ctransC
transdataAux fun c1 d c2 =
if noneed
then d
else unsafePerformIO $ do
v <- createVector (dim d)
withForeignPtr (fptr d) $ \pd ->
withForeignPtr (fptr v) $ \pv ->
fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
return v
where r1 = dim d `div` c1
r2 = dim d `div` c2
noneed = r1 == 1 || c1 == 1
foreign import ccall "auxi.h transR" ctransR :: TMM
foreign import ccall "auxi.h transC" ctransC :: TCMCM
subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double
subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do
r <- createMatrix RowMajor rt ct
let x = cmat x'
app2 (c_submatrixR (fi r0) (fi $ r0+rt1) (fi c0) (fi $ c0+ct1)) mat x mat r "subMatrixR"
return r
foreign import ccall "auxi.h submatrixR" c_submatrixR :: CInt -> CInt -> CInt -> CInt -> TMM
subMatrixC :: (Int,Int) -> (Int,Int) -> Matrix (Complex Double) -> Matrix (Complex Double)
subMatrixC (r0,c0) (rt,ct) x =
reshape ct . asComplex . flatten .
subMatrixR (r0,2*c0) (rt,2*ct) .
reshape (2*cols x) . asReal . flatten $ x
subMatrix :: Element a
=> (Int,Int)
-> (Int,Int)
-> Matrix a
-> Matrix a
subMatrix = subMatrixD
diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
m <- createMatrix RowMajor n n
app2 fun vec v mat m msg
return m
diagR :: Vector Double -> Matrix Double
diagR = diagAux c_diagR "diagR"
foreign import ccall "auxi.h diagR" c_diagR :: TVM
diagC :: Vector (Complex Double) -> Matrix (Complex Double)
diagC = diagAux c_diagC "diagC"
foreign import ccall "auxi.h diagC" c_diagC :: TCVCM
diag :: Element a => Vector a -> Matrix a
diag = diagD
constantAux fun x n = unsafePerformIO $ do
v <- createVector n
px <- newArray [x]
app1 (fun px) vec v "constantAux"
free px
return v
constantR :: Double -> Int -> Vector Double
constantR = constantAux cconstantR
foreign import ccall "auxi.h constantR" cconstantR :: Ptr Double -> TV
constantC :: Complex Double -> Int -> Vector (Complex Double)
constantC = constantAux cconstantC
foreign import ccall "auxi.h constantC" cconstantC :: Ptr (Complex Double) -> TCV
constant :: Element a => a -> Int -> Vector a
constant = constantD
conj :: Vector (Complex Double) -> Vector (Complex Double)
conj v = unsafePerformIO $ do
r <- createVector (dim v)
app2 cconjugate vec v vec r "cconjugate"
return r
foreign import ccall "auxi.h conjugate" cconjugate :: TCVCV
toComplex :: (Vector Double, Vector Double) -> Vector (Complex Double)
toComplex (r,i) = asComplex $ flatten $ fromColumns [r,i]
fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double)
fromComplex z = (r,i) where
[r,i] = toColumns $ reshape 2 $ asReal z
comp :: Vector Double -> Vector (Complex Double)
comp v = toComplex (v,constant 0 (dim v))
fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
fromFile filename (r,c) = do
charname <- newCString filename
res <- createMatrix RowMajor r c
app1 (c_gslReadMatrix charname) mat res "gslReadMatrix"
return res
foreign import ccall "auxi.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM