1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
{-# OPTIONS_GHC -fglasgow-exts #-}
-----------------------------------------------------------------------------
-- |
-- Module : LAPACK.Internal
-- Copyright : (c) Alberto Ruiz 2006-7
-- License : GPL-style
--
-- Maintainer : Alberto Ruiz (aruiz at um dot es)
-- Stability : provisional
-- Portability : portable (uses FFI)
--
-- Wrappers for a few LAPACK functions (<http://www.netlib.org/lapack>).
--
-----------------------------------------------------------------------------
module LAPACK.Internal where
import Data.Packed.Internal.Vector
import Data.Packed.Internal.Matrix
import Complex
import Foreign
import Foreign.C.Types
import Foreign.C.String
-----------------------------------------------------------------------------
-- dgesvd
foreign import ccall "lapack-aux.h svd_l_R"
dgesvd :: Double ::> Double ::> (Double :> Double ::> IO Int)
-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
--
-- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
svdR :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
svdR x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdR' x
svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
u <- createMatrix ColumnMajor r r
s <- createVector (min r c)
v <- createMatrix ColumnMajor c c
dgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdR" [fdat x]
return (u,s,trans v)
-----------------------------------------------------------------------------
-- dgesdd
foreign import ccall "lapack-aux.h svd_l_Rdd"
dgesdd :: Double ::> Double ::> (Double :> Double ::> IO Int)
-----------------------------------------------------------------------------
-- zgesvd
foreign import ccall "lapack-aux.h svd_l_C"
zgesvd :: (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
--
-- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@.
svdC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double))
svdC x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdC' x
svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
u <- createMatrix ColumnMajor r r
s <- createVector (min r c)
v <- createMatrix ColumnMajor c c
zgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdC" [fdat x]
return (u,s,trans v)
-----------------------------------------------------------------------------
-- zgeev
foreign import ccall "lapack-aux.h eig_l_C"
zgeev :: (Complex Double) ::> (Complex Double) ::> ((Complex Double) :> (Complex Double) ::> IO Int)
-- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
--
-- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@.
--
-- The eigenvectors are the columns of v.
-- The eigenvalues are not sorted.
eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
eigC (m@M {rows = r})
| r == 1 = (fromList [cdat m `at` 0], singleton 1)
| otherwise = unsafePerformIO $ do
l <- createVector r
v <- createMatrix ColumnMajor r r
dummy <- createMatrix ColumnMajor 1 1
zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m]
return (l,v)
-----------------------------------------------------------------------------
-- dgeev
foreign import ccall "lapack-aux.h eig_l_R"
dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int)
-----------------------------------------------------------------------------
-- dsyev
foreign import ccall "lapack-aux.h eig_l_S"
dsyev :: Double ::> (Double :> Double ::> IO Int)
-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
--
-- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@.
--
-- The eigenvectors are the columns of v.
-- The eigenvalues are sorted in descending order (use eigS' for ascending order).
eigS :: Matrix Double -> (Vector Double, Matrix Double)
eigS m = (s', fliprl v)
where (s,v) = eigS' m
s' = fromList . reverse . toList $ s
eigS' (m@M {rows = r}) = unsafePerformIO $ do
l <- createVector r
v <- createMatrix ColumnMajor r r
dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m]
return (l,v)
-----------------------------------------------------------------------------
-- zheev
foreign import ccall "lapack-aux.h eig_l_H"
zheev :: (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
-----------------------------------------------------------------------------
-- dgesv
foreign import ccall "lapack-aux.h linearSolveR_l"
dgesv :: Double ::> Double ::> Double ::> IO Int
-----------------------------------------------------------------------------
-- zgesv
foreign import ccall "lapack-aux.h linearSolveC_l"
zgesv :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
-----------------------------------------------------------------------------------
-- dgels
foreign import ccall "lapack-aux.h linearSolveLSR_l"
dgels :: Double ::> Double ::> Double ::> IO Int
-- | Wrapper for LAPACK's /dgels/, which obtains the least squared error solution of an overconstrained real linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDR'.
linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSR_l a b
linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
r <- createMatrix ColumnMajor (max m n) nrhs
dgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSR" [fdat a, fdat b]
return r
-----------------------------------------------------------------------------------
-- zgels
foreign import ccall "lapack-aux.h linearSolveLSC_l"
zgels :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
-----------------------------------------------------------------------------------
-- dgelss
foreign import ccall "lapack-aux.h linearSolveSVDR_l"
dgelss :: Double -> Double ::> Double ::> Double ::> IO Int
-----------------------------------------------------------------------------------
-- zgelss
foreign import ccall "lapack-aux.h linearSolveSVDC_l"
zgelss :: Double -> (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
|
