module Numeric.LinearAlgebra.LAPACK (
svdR, svdRdd, svdC,
eigC, eigR, eigS, eigH, eigS', eigH',
linearSolveR, linearSolveC,
linearSolveLSR, linearSolveLSC,
linearSolveSVDR, linearSolveSVDC,
luR, luC,
cholS, cholH,
qrR, qrC,
hessR, hessC,
schurR, schurC
) where
import Data.Packed.Internal
import Data.Packed.Internal.Vector
import Data.Packed.Internal.Matrix
import Data.Packed.Vector
import Data.Packed.Matrix
import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale))
import Complex
import Foreign
foreign import ccall "LAPACK/lapack-aux.h svd_l_R" dgesvd :: TMMVM
foreign import ccall "LAPACK/lapack-aux.h svd_l_C" zgesvd :: TCMCMVCM
foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd" dgesdd :: TMMVM
svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdR = svdAux dgesvd "svdR" . fmat
svdRdd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdRdd = svdAux dgesdd "svdRdd" . fmat
svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdC = svdAux zgesvd "svdC" . fmat
svdAux f st x = unsafePerformIO $ do
u <- createMatrix ColumnMajor r r
s <- createVector (min r c)
v <- createMatrix ColumnMajor c c
app4 f mat x mat u vec s mat v st
return (u,s,trans v)
where r = rows x
c = cols x
eigAux f st m
| r == 1 = (fromList [flatten m `at` 0], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix ColumnMajor r r
dummy <- createMatrix ColumnMajor 1 1
app4 f mat m mat dummy vec l mat v st
return (l,v)
where r = rows m
foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
eigC = eigAux zgeev "eigC" . fmat
eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
eigR m = (s', v'')
where (s,v) = eigRaux (fmat m)
s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
v' = toRows $ trans v
v'' = fromColumns $ fixeig (toList s') v'
r = rows m
eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
eigRaux m
| r == 1 = (fromList [(flatten m `at` 0):+0], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix ColumnMajor r r
dummy <- createMatrix ColumnMajor 1 1
app4 dgeev mat m mat dummy vec l mat v "eigR"
return (l,v)
where r = rows m
fixeig [] _ = []
fixeig [_] [v] = [comp v]
fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
| r1 == r2 && i1 == (i2) = toComplex (v1,v2) : toComplex (v1,scale (1) v2) : fixeig r vs
| otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs)
where scale = vectorMapValR Scale
fixeig _ _ = error "fixeig with impossible inputs"
eigS :: Matrix Double -> (Vector Double, Matrix Double)
eigS m = (s', fliprl v)
where (s,v) = eigS' (fmat m)
s' = fromList . reverse . toList $ s
eigS' m
| r == 1 = (fromList [flatten m `at` 0], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix ColumnMajor r r
app3 dsyev mat m vec l mat v "eigS"
return (l,v)
where r = rows m
eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH m = (s', fliprl v)
where (s,v) = eigH' (fmat m)
s' = fromList . reverse . toList $ s
eigH' m
| r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix ColumnMajor r r
app3 zheev mat m vec l mat v "eigH"
return (l,v)
where r = rows m
foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l" zgesv :: TCMCMCM
linearSolveSQAux f st a b
| n1==n2 && n1==r = unsafePerformIO $ do
s <- createMatrix ColumnMajor r c
app3 f mat a mat b mat s st
return s
| otherwise = error $ st ++ " of nonsquare matrix"
where n1 = rows a
n2 = cols a
r = rows b
c = cols b
linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b)
linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b)
foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l" dgels :: TMMM
foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l" zgels :: TCMCMCM
foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l" dgelss :: Double -> TMMM
foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
linearSolveAux f st a b = unsafePerformIO $ do
r <- createMatrix ColumnMajor (max m n) nrhs
app3 f mat a mat b mat r st
return r
where m = rows a
n = cols a
nrhs = cols b
linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $
linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b)
linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $
linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b)
linearSolveSVDR :: Maybe Double
-> Matrix Double
-> Matrix Double
-> Matrix Double
linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b)
linearSolveSVDR Nothing a b = linearSolveSVDR (Just (1)) (fmat a) (fmat b)
linearSolveSVDC :: Maybe Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b)
linearSolveSVDC Nothing a b = linearSolveSVDC (Just (1)) (fmat a) (fmat b)
foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM
foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
cholH = cholAux zpotrf "cholH" . fmat
cholS :: Matrix Double -> Matrix Double
cholS = cholAux dpotrf "cholS" . fmat
cholAux f st a = unsafePerformIO $ do
r <- createMatrix ColumnMajor n n
app2 f mat a mat r st
return r
where n = rows a
foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM
foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM
qrR :: Matrix Double -> (Matrix Double, Vector Double)
qrR = qrAux dgeqr2 "qrR" . fmat
qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
qrC = qrAux zgeqr2 "qrC" . fmat
qrAux f st a = unsafePerformIO $ do
r <- createMatrix ColumnMajor m n
tau <- createVector mn
app3 f mat a vec tau mat r st
return (r,tau)
where m = rows a
n = cols a
mn = min m n
foreign import ccall "LAPACK/lapack-aux.h hess_l_R" dgehrd :: TMVM
foreign import ccall "LAPACK/lapack-aux.h hess_l_C" zgehrd :: TCMCVCM
hessR :: Matrix Double -> (Matrix Double, Vector Double)
hessR = hessAux dgehrd "hessR" . fmat
hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
hessC = hessAux zgehrd "hessC" . fmat
hessAux f st a = unsafePerformIO $ do
r <- createMatrix ColumnMajor m n
tau <- createVector (mn1)
app3 f mat a vec tau mat r st
return (r,tau)
where m = rows a
n = cols a
mn = min m n
foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM
schurR :: Matrix Double -> (Matrix Double, Matrix Double)
schurR = schurAux dgees "schurR" . fmat
schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
schurC = schurAux zgees "schurC" . fmat
schurAux f st a = unsafePerformIO $ do
u <- createMatrix ColumnMajor n n
s <- createMatrix ColumnMajor n n
app3 f mat a mat u mat s st
return (u,s)
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
luR :: Matrix Double -> (Matrix Double, [Int])
luR = luAux dgetrf "luR" . fmat
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