module Numeric.GSL.Matrix(
eigSg, eigHg,
svdg,
qr, qrPacked, unpackQR,
cholR, cholC,
luSolveR, luSolveC,
luR, luC
) where
import Data.Packed.Internal
import Data.Packed.Matrix(ident)
import Numeric.GSL.Vector
import Foreign
import Complex
eigSg :: Matrix Double -> (Vector Double, Matrix Double)
eigSg = eigSg' . cmat
eigSg' m
| r == 1 = (fromList [cdat m `at` 0], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix RowMajor r r
app3 c_eigS mat m vec l mat v "eigSg"
return (l,v)
where r = rows m
foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double))
eigHg = eigHg' . cmat
eigHg' m
| r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix RowMajor r r
app3 c_eigH mat m vec l mat v "eigHg"
return (l,v)
where r = rows m
foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdg x = if rows x >= cols x
then svd' (cmat x)
else (v, s, u) where (u,s,v) = svd' (cmat (trans x))
svd' x = unsafePerformIO $ do
u <- createMatrix RowMajor r c
s <- createVector c
v <- createMatrix RowMajor c c
app4 c_svd mat x mat u vec s mat v "svdg"
return (u,s,v)
where r = rows x
c = cols x
foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
qr :: Matrix Double -> (Matrix Double, Matrix Double)
qr = qr' . cmat
qr' x = unsafePerformIO $ do
q <- createMatrix RowMajor r r
rot <- createMatrix RowMajor r c
app3 c_qr mat x mat q mat rot "qr"
return (q,rot)
where r = rows x
c = cols x
foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
qrPacked :: Matrix Double -> (Matrix Double, Vector Double)
qrPacked = qrPacked' . cmat
qrPacked' x = unsafePerformIO $ do
qrp <- createMatrix RowMajor r c
tau <- createVector (min r c)
app3 c_qrPacked mat x mat qrp vec tau "qrUnpacked"
return (qrp,tau)
where r = rows x
c = cols x
foreign import ccall "gsl-aux.h QRpacked" c_qrPacked :: TMMV
unpackQR :: (Matrix Double, Vector Double) -> (Matrix Double, Matrix Double)
unpackQR (qrp,tau) = unpackQR' (cmat qrp, tau)
unpackQR' (qrp,tau) = unsafePerformIO $ do
q <- createMatrix RowMajor r r
res <- createMatrix RowMajor r c
app4 c_qrUnpack mat qrp vec tau mat q mat res "qrUnpack"
return (q,res)
where r = rows qrp
c = cols qrp
foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM
type TMMV = Int -> Int -> PD -> TMV
type TMVMM = Int -> Int -> PD -> Int -> PD -> TMM
cholR :: Matrix Double -> Matrix Double
cholR = cholR' . cmat
cholR' x = unsafePerformIO $ do
r <- createMatrix RowMajor n n
app2 c_cholR mat x mat r "cholR"
return r
where n = rows x
foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
cholC = cholC' . cmat
cholC' x = unsafePerformIO $ do
r <- createMatrix RowMajor n n
app2 c_cholC mat x mat r "cholC"
return r
where n = rows x
foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
luSolveR :: Matrix Double -> Matrix Double -> Matrix Double
luSolveR a b = luSolveR' (cmat a) (cmat b)
luSolveR' a b
| n1==n2 && n1==r = unsafePerformIO $ do
s <- createMatrix RowMajor r c
app3 c_luSolveR mat a mat b mat s "luSolveR"
return s
| otherwise = error "luSolveR of nonsquare matrix"
where n1 = rows a
n2 = cols a
r = rows b
c = cols b
foreign import ccall "gsl-aux.h luSolveR" c_luSolveR :: TMMM
luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
luSolveC a b = luSolveC' (cmat a) (cmat b)
luSolveC' a b
| n1==n2 && n1==r = unsafePerformIO $ do
s <- createMatrix RowMajor r c
app3 c_luSolveC mat a mat b mat s "luSolveC"
return s
| otherwise = error "luSolveC of nonsquare matrix"
where n1 = rows a
n2 = cols a
r = rows b
c = cols b
foreign import ccall "gsl-aux.h luSolveC" c_luSolveC :: TCMCMCM
luRaux :: Matrix Double -> Vector Double
luRaux = luRaux' . cmat
luRaux' x = unsafePerformIO $ do
res <- createVector (r*r+r+1)
app2 c_luRaux mat x vec res "luRaux"
return res
where r = rows x
foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV
luCaux :: Matrix (Complex Double) -> Vector (Complex Double)
luCaux = luCaux' . cmat
luCaux' x = unsafePerformIO $ do
res <- createVector (r*r+r+1)
app2 c_luCaux mat x vec res "luCaux"
return res
where r = rows x
foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV
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
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` liftMatrix comp (ident r)
add = liftMatrix2 $ vectorZipC Add
mul = liftMatrix2 $ vectorZipC Mul
triang r c h v = reshape c $ fromList [el i j | i<-[0..r1], j<-[0..c1]]
where el i j = if ji>=h then v else 1 v