summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--examples/tests.hs10
-rw-r--r--lib/LAPACK.hs261
-rw-r--r--lib/LAPACK/Internal.hs262
3 files changed, 265 insertions, 268 deletions
diff --git a/examples/tests.hs b/examples/tests.hs
index 047794c..b075704 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -159,10 +159,10 @@ addV v1 v2 = fromList $ zipWith (+) (toList v1) (toList v2)
159 159
160type BaseType = Double 160type BaseType = Double
161 161
162svdTestR prod m = u <> s <> trans v |~| m 162svdTestR fun prod m = u <> s <> trans v |~| m
163 && u <> trans u |~| ident (rows m) 163 && u <> trans u |~| ident (rows m)
164 && v <> trans v |~| ident (cols m) 164 && v <> trans v |~| ident (cols m)
165 where (u,s,v) = svdR m 165 where (u,s,v) = fun m
166 (<>) = prod 166 (<>) = prod
167 167
168 168
@@ -243,8 +243,10 @@ main = do
243 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| trans (mulF (trans m2) (trans m1 :: Matrix BaseType)) 243 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| trans (mulF (trans m2) (trans m1 :: Matrix BaseType))
244 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| multiplyG m1 (m2 :: Matrix BaseType) 244 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 |=| multiplyG m1 (m2 :: Matrix BaseType)
245 putStrLn "--------- SVD ---------" 245 putStrLn "--------- SVD ---------"
246 quickCheck (svdTestR mulC) 246 quickCheck (svdTestR svdR mulC)
247 quickCheck (svdTestR mulF) 247 quickCheck (svdTestR svdR mulF)
248 quickCheck (svdTestR svdRdd mulC)
249 quickCheck (svdTestR svdRdd mulF)
248 quickCheck (svdTestC mulC) 250 quickCheck (svdTestC mulC)
249 quickCheck (svdTestC mulF) 251 quickCheck (svdTestC mulF)
250 putStrLn "--------- EIG ---------" 252 putStrLn "--------- EIG ---------"
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index 0019fbe..e6249a1 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -1,3 +1,4 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
1----------------------------------------------------------------------------- 2-----------------------------------------------------------------------------
2-- | 3-- |
3-- Module : LAPACK 4-- Module : LAPACK
@@ -13,11 +14,267 @@
13----------------------------------------------------------------------------- 14-----------------------------------------------------------------------------
14 15
15module LAPACK ( 16module LAPACK (
16 svdR, svdR', svdC, svdC', 17 svdR, svdR', svdRdd, svdRdd', svdC, svdC',
17 eigC, eigR, eigS, eigH, 18 eigC, eigR, eigS, eigH,
18 linearSolveR, linearSolveC, 19 linearSolveR, linearSolveC,
19 linearSolveLSR, linearSolveLSC, 20 linearSolveLSR, linearSolveLSC,
20 linearSolveSVDR, linearSolveSVDC, 21 linearSolveSVDR, linearSolveSVDC,
21) where 22) where
22 23
23import LAPACK.Internal 24import Data.Packed.Internal.Vector
25import Data.Packed.Internal.Matrix
26import Complex
27import Foreign
28
29-----------------------------------------------------------------------------
30-- dgesvd
31foreign import ccall "LAPACK/lapack-aux.h svd_l_R"
32 dgesvd :: Double ::> Double ::> (Double :> Double ::> IO Int)
33
34-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
35--
36-- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
37svdR :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
38svdR x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdR' x
39
40svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
41 u <- createMatrix ColumnMajor r r
42 s <- createVector (min r c)
43 v <- createMatrix ColumnMajor c c
44 dgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdR" [fdat x]
45 return (u,s,trans v)
46
47-----------------------------------------------------------------------------
48-- dgesdd
49foreign import ccall "LAPACK/lapack-aux.h svd_l_Rdd"
50 dgesdd :: Double ::> Double ::> (Double :> Double ::> IO Int)
51
52-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
53--
54-- @(u,s,v)=svdRdd m@ so that @m=u \<\> s \<\> 'trans' v@.
55svdRdd :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
56svdRdd x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdRdd' x
57
58svdRdd' x@M {rows = r, cols = c} = unsafePerformIO $ do
59 u <- createMatrix ColumnMajor r r
60 s <- createVector (min r c)
61 v <- createMatrix ColumnMajor c c
62 dgesdd // mat fdat x // mat dat u // vec s // mat dat v // check "svdRdd" [fdat x]
63 return (u,s,trans v)
64
65-----------------------------------------------------------------------------
66-- zgesvd
67foreign import ccall "LAPACK/lapack-aux.h svd_l_C"
68 zgesvd :: (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
69
70-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
71--
72-- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@.
73svdC :: Matrix (Complex Double)
74 -> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double))
75svdC x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdC' x
76
77svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
78 u <- createMatrix ColumnMajor r r
79 s <- createVector (min r c)
80 v <- createMatrix ColumnMajor c c
81 zgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdC" [fdat x]
82 return (u,s,trans v)
83
84-----------------------------------------------------------------------------
85-- zgeev
86foreign import ccall "LAPACK/lapack-aux.h eig_l_C"
87 zgeev :: (Complex Double) ::> (Complex Double) ::> ((Complex Double) :> (Complex Double) ::> IO Int)
88
89-- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
90--
91-- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@.
92--
93-- The eigenvectors are the columns of v.
94-- The eigenvalues are not sorted.
95eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
96eigC (m@M {rows = r})
97 | r == 1 = (fromList [cdat m `at` 0], singleton 1)
98 | otherwise = unsafePerformIO $ do
99 l <- createVector r
100 v <- createMatrix ColumnMajor r r
101 dummy <- createMatrix ColumnMajor 1 1
102 zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m]
103 return (l,v)
104
105-----------------------------------------------------------------------------
106-- dgeev
107foreign import ccall "LAPACK/lapack-aux.h eig_l_R"
108 dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int)
109
110-- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
111--
112-- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@.
113--
114-- The eigenvectors are the columns of v.
115-- The eigenvalues are not sorted.
116eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
117eigR (m@M {rows = r}) = (s', v'')
118 where (s,v) = eigRaux m
119 s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
120 v' = toRows $ trans v
121 v'' = fromColumns $ fixeig (toList s') v'
122
123eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
124eigRaux (m@M {rows = r})
125 | r == 1 = (fromList [(cdat m `at` 0):+0], singleton 1)
126 | otherwise = unsafePerformIO $ do
127 l <- createVector r
128 v <- createMatrix ColumnMajor r r
129 dummy <- createMatrix ColumnMajor 1 1
130 dgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigR" [fdat m]
131 return (l,v)
132
133fixeig [] _ = []
134fixeig [r] [v] = [comp v]
135fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
136 | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs
137 | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs)
138
139scale r v = fromList [r] `outer` v
140
141-----------------------------------------------------------------------------
142-- dsyev
143foreign import ccall "LAPACK/lapack-aux.h eig_l_S"
144 dsyev :: Double ::> (Double :> Double ::> IO Int)
145
146-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
147--
148-- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@.
149--
150-- The eigenvectors are the columns of v.
151-- The eigenvalues are sorted in descending order (use eigS' for ascending order).
152eigS :: Matrix Double -> (Vector Double, Matrix Double)
153eigS m = (s', fliprl v)
154 where (s,v) = eigS' m
155 s' = fromList . reverse . toList $ s
156
157eigS' (m@M {rows = r})
158 | r == 1 = (fromList [cdat m `at` 0], singleton 1)
159 | otherwise = unsafePerformIO $ do
160 l <- createVector r
161 v <- createMatrix ColumnMajor r r
162 dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m]
163 return (l,v)
164
165-----------------------------------------------------------------------------
166-- zheev
167foreign import ccall "LAPACK/lapack-aux.h eig_l_H"
168 zheev :: (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
169
170-- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix:
171--
172-- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@.
173--
174-- The eigenvectors are the columns of v.
175-- The eigenvalues are sorted in descending order.
176eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
177eigH m = (s', fliprl v)
178 where (s,v) = eigH' m
179 s' = fromList . reverse . toList $ s
180
181eigH' (m@M {rows = r})
182 | r == 1 = (fromList [realPart (cdat m `at` 0)], singleton 1)
183 | otherwise = unsafePerformIO $ do
184 l <- createVector r
185 v <- createMatrix ColumnMajor r r
186 zheev // mat fdat m // vec l // mat dat v // check "eigH" [fdat m]
187 return (l,v)
188
189-----------------------------------------------------------------------------
190-- dgesv
191foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l"
192 dgesv :: Double ::> Double ::> Double ::> IO Int
193
194-- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition.
195linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
196linearSolveR a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
197 | n1==n2 && n1==r = unsafePerformIO $ do
198 s <- createMatrix ColumnMajor r c
199 dgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveR" [fdat a, fdat b]
200 return s
201 | otherwise = error "linearSolveR of nonsquare matrix"
202
203-----------------------------------------------------------------------------
204-- zgesv
205foreign import ccall "LAPACK/lapack-aux.h linearSolveC_l"
206 zgesv :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
207
208-- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition.
209linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
210linearSolveC a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
211 | n1==n2 && n1==r = unsafePerformIO $ do
212 s <- createMatrix ColumnMajor r c
213 zgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveC" [fdat a, fdat b]
214 return s
215 | otherwise = error "linearSolveC of nonsquare matrix"
216
217-----------------------------------------------------------------------------------
218-- dgels
219foreign import ccall "LAPACK/lapack-aux.h linearSolveLSR_l"
220 dgels :: Double ::> Double ::> Double ::> IO Int
221
222-- | 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'.
223linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
224linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSR_l a b
225
226linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
227 r <- createMatrix ColumnMajor (max m n) nrhs
228 dgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSR" [fdat a, fdat b]
229 return r
230
231-----------------------------------------------------------------------------------
232-- zgels
233foreign import ccall "LAPACK/lapack-aux.h linearSolveLSC_l"
234 zgels :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
235
236-- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'.
237linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
238linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSC_l a b
239
240linearSolveLSC_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
241 r <- createMatrix ColumnMajor (max m n) nrhs
242 zgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSC" [fdat a, fdat b]
243 return r
244
245-----------------------------------------------------------------------------------
246-- dgelss
247foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDR_l"
248 dgelss :: Double -> Double ::> Double ::> Double ::> IO Int
249
250-- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
251linearSolveSVDR :: Maybe Double -- ^ rcond
252 -> Matrix Double -- ^ coefficient matrix
253 -> Matrix Double -- ^ right hand sides (as columns)
254 -> Matrix Double -- ^ solution vectors (as columns)
255linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDR_l rcond a b
256linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) a b
257
258linearSolveSVDR_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
259 r <- createMatrix ColumnMajor (max m n) nrhs
260 dgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDR" [fdat a, fdat b]
261 return r
262
263-----------------------------------------------------------------------------------
264-- zgelss
265foreign import ccall "LAPACK/lapack-aux.h linearSolveSVDC_l"
266 zgelss :: Double -> (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
267
268-- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
269linearSolveSVDC :: Maybe Double -- ^ rcond
270 -> Matrix (Complex Double) -- ^ coefficient matrix
271 -> Matrix (Complex Double) -- ^ right hand sides (as columns)
272 -> Matrix (Complex Double) -- ^ solution vectors (as columns)
273linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDC_l rcond a b
274linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) a b
275
276linearSolveSVDC_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
277 r <- createMatrix ColumnMajor (max m n) nrhs
278 zgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDC" [fdat a, fdat b]
279 return r
280
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs
deleted file mode 100644
index ec46b66..0000000
--- a/lib/LAPACK/Internal.hs
+++ /dev/null
@@ -1,262 +0,0 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : LAPACK.Internal
5-- Copyright : (c) Alberto Ruiz 2006-7
6-- License : GPL-style
7--
8-- Maintainer : Alberto Ruiz (aruiz at um dot es)
9-- Stability : provisional
10-- Portability : portable (uses FFI)
11--
12-- Wrappers for a few LAPACK functions (<http://www.netlib.org/lapack>).
13--
14-----------------------------------------------------------------------------
15
16module LAPACK.Internal where
17
18import Data.Packed.Internal.Vector
19import Data.Packed.Internal.Matrix
20import Complex
21import Foreign
22import Foreign.C.Types
23import Foreign.C.String
24
25-----------------------------------------------------------------------------
26-- dgesvd
27foreign import ccall "lapack-aux.h svd_l_R"
28 dgesvd :: Double ::> Double ::> (Double :> Double ::> IO Int)
29
30-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
31--
32-- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
33svdR :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
34svdR x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdR' x
35
36svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
37 u <- createMatrix ColumnMajor r r
38 s <- createVector (min r c)
39 v <- createMatrix ColumnMajor c c
40 dgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdR" [fdat x]
41 return (u,s,trans v)
42
43-----------------------------------------------------------------------------
44-- dgesdd
45foreign import ccall "lapack-aux.h svd_l_Rdd"
46 dgesdd :: Double ::> Double ::> (Double :> Double ::> IO Int)
47
48-----------------------------------------------------------------------------
49-- zgesvd
50foreign import ccall "lapack-aux.h svd_l_C"
51 zgesvd :: (Complex Double) ::> (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
52
53-- | Wrapper for LAPACK's /zgesvd/, which computes the full svd decomposition of a complex matrix.
54--
55-- @(u,s,v)=svdC m@ so that @m=u \<\> s \<\> 'trans' v@.
56svdC :: Matrix (Complex Double)
57 -> (Matrix (Complex Double), Matrix Double, Matrix (Complex Double))
58svdC x@M {rows = r, cols = c} = (u, diagRect s r c, v) where (u,s,v) = svdC' x
59
60svdC' x@M {rows = r, cols = c} = unsafePerformIO $ do
61 u <- createMatrix ColumnMajor r r
62 s <- createVector (min r c)
63 v <- createMatrix ColumnMajor c c
64 zgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdC" [fdat x]
65 return (u,s,trans v)
66
67-----------------------------------------------------------------------------
68-- zgeev
69foreign import ccall "lapack-aux.h eig_l_C"
70 zgeev :: (Complex Double) ::> (Complex Double) ::> ((Complex Double) :> (Complex Double) ::> IO Int)
71
72-- | Wrapper for LAPACK's /zgeev/, which computes the eigenvalues and right eigenvectors of a general complex matrix:
73--
74-- if @(l,v)=eigC m@ then @m \<\> v = v \<\> diag l@.
75--
76-- The eigenvectors are the columns of v.
77-- The eigenvalues are not sorted.
78eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
79eigC (m@M {rows = r})
80 | r == 1 = (fromList [cdat m `at` 0], singleton 1)
81 | otherwise = unsafePerformIO $ do
82 l <- createVector r
83 v <- createMatrix ColumnMajor r r
84 dummy <- createMatrix ColumnMajor 1 1
85 zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m]
86 return (l,v)
87
88-----------------------------------------------------------------------------
89-- dgeev
90foreign import ccall "lapack-aux.h eig_l_R"
91 dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int)
92
93-- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix:
94--
95-- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@.
96--
97-- The eigenvectors are the columns of v.
98-- The eigenvalues are not sorted.
99eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
100eigR (m@M {rows = r}) = (s', v'')
101 where (s,v) = eigRaux m
102 s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
103 v' = toRows $ trans v
104 v'' = fromColumns $ fixeig (toList s') v'
105
106eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
107eigRaux (m@M {rows = r})
108 | r == 1 = (fromList [(cdat m `at` 0):+0], singleton 1)
109 | otherwise = unsafePerformIO $ do
110 l <- createVector r
111 v <- createMatrix ColumnMajor r r
112 dummy <- createMatrix ColumnMajor 1 1
113 dgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigR" [fdat m]
114 return (l,v)
115
116fixeig [] _ = []
117fixeig [r] [v] = [comp v]
118fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
119 | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs
120 | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs)
121
122scale r v = fromList [r] `outer` v
123
124-----------------------------------------------------------------------------
125-- dsyev
126foreign import ccall "lapack-aux.h eig_l_S"
127 dsyev :: Double ::> (Double :> Double ::> IO Int)
128
129-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
130--
131-- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@.
132--
133-- The eigenvectors are the columns of v.
134-- The eigenvalues are sorted in descending order (use eigS' for ascending order).
135eigS :: Matrix Double -> (Vector Double, Matrix Double)
136eigS m = (s', fliprl v)
137 where (s,v) = eigS' m
138 s' = fromList . reverse . toList $ s
139
140eigS' (m@M {rows = r})
141 | r == 1 = (fromList [cdat m `at` 0], singleton 1)
142 | otherwise = unsafePerformIO $ do
143 l <- createVector r
144 v <- createMatrix ColumnMajor r r
145 dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m]
146 return (l,v)
147
148-----------------------------------------------------------------------------
149-- zheev
150foreign import ccall "lapack-aux.h eig_l_H"
151 zheev :: (Complex Double) ::> (Double :> (Complex Double) ::> IO Int)
152
153-- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix:
154--
155-- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@.
156--
157-- The eigenvectors are the columns of v.
158-- The eigenvalues are sorted in descending order.
159eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
160eigH m = (s', fliprl v)
161 where (s,v) = eigH' m
162 s' = fromList . reverse . toList $ s
163
164eigH' (m@M {rows = r})
165 | r == 1 = (fromList [realPart (cdat m `at` 0)], singleton 1)
166 | otherwise = unsafePerformIO $ do
167 l <- createVector r
168 v <- createMatrix ColumnMajor r r
169 zheev // mat fdat m // vec l // mat dat v // check "eigH" [fdat m]
170 return (l,v)
171
172-----------------------------------------------------------------------------
173-- dgesv
174foreign import ccall "lapack-aux.h linearSolveR_l"
175 dgesv :: Double ::> Double ::> Double ::> IO Int
176
177-- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition.
178linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
179linearSolveR a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
180 | n1==n2 && n1==r = unsafePerformIO $ do
181 s <- createMatrix ColumnMajor r c
182 dgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveR" [fdat a, fdat b]
183 return s
184 | otherwise = error "linearSolveR of nonsquare matrix"
185
186-----------------------------------------------------------------------------
187-- zgesv
188foreign import ccall "lapack-aux.h linearSolveC_l"
189 zgesv :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
190
191-- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition.
192linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
193linearSolveC a@(M {rows = n1, cols = n2}) b@(M {rows = r, cols = c})
194 | n1==n2 && n1==r = unsafePerformIO $ do
195 s <- createMatrix ColumnMajor r c
196 zgesv // mat fdat a // mat fdat b // mat dat s // check "linearSolveC" [fdat a, fdat b]
197 return s
198 | otherwise = error "linearSolveC of nonsquare matrix"
199
200-----------------------------------------------------------------------------------
201-- dgels
202foreign import ccall "lapack-aux.h linearSolveLSR_l"
203 dgels :: Double ::> Double ::> Double ::> IO Int
204
205-- | 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'.
206linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
207linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSR_l a b
208
209linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
210 r <- createMatrix ColumnMajor (max m n) nrhs
211 dgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSR" [fdat a, fdat b]
212 return r
213
214-----------------------------------------------------------------------------------
215-- zgels
216foreign import ccall "lapack-aux.h linearSolveLSC_l"
217 zgels :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
218
219-- | Wrapper for LAPACK's /zgels/, which obtains the least squared error solution of an overconstrained complex linear system or the minimum norm solution of an underdetermined system, for several right-hand sides. For rank deficient systems use 'linearSolveSVDC'.
220linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
221linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSC_l a b
222
223linearSolveLSC_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
224 r <- createMatrix ColumnMajor (max m n) nrhs
225 zgels // mat fdat a // mat fdat b // mat dat r // check "linearSolveLSC" [fdat a, fdat b]
226 return r
227
228-----------------------------------------------------------------------------------
229-- dgelss
230foreign import ccall "lapack-aux.h linearSolveSVDR_l"
231 dgelss :: Double -> Double ::> Double ::> Double ::> IO Int
232
233-- | Wrapper for LAPACK's /dgelss/, which obtains the minimum norm solution to a real linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
234linearSolveSVDR :: Maybe Double -- ^ rcond
235 -> Matrix Double -- ^ coefficient matrix
236 -> Matrix Double -- ^ right hand sides (as columns)
237 -> Matrix Double -- ^ solution vectors (as columns)
238linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDR_l rcond a b
239linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) a b
240
241linearSolveSVDR_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
242 r <- createMatrix ColumnMajor (max m n) nrhs
243 dgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDR" [fdat a, fdat b]
244 return r
245
246-----------------------------------------------------------------------------------
247-- zgelss
248foreign import ccall "lapack-aux.h linearSolveSVDC_l"
249 zgelss :: Double -> (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int
250
251-- | Wrapper for LAPACK's /zgelss/, which obtains the minimum norm solution to a complex linear least squares problem Ax=B using the svd, for several right-hand sides. Admits rank deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
252linearSolveSVDC :: Maybe Double -- ^ rcond
253 -> Matrix (Complex Double) -- ^ coefficient matrix
254 -> Matrix (Complex Double) -- ^ right hand sides (as columns)
255 -> Matrix (Complex Double) -- ^ solution vectors (as columns)
256linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDC_l rcond a b
257linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) a b
258
259linearSolveSVDC_l rcond a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformIO $ do
260 r <- createMatrix ColumnMajor (max m n) nrhs
261 zgelss rcond // mat fdat a // mat fdat b // mat dat r // check "linearSolveSVDC" [fdat a, fdat b]
262 return r
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 09:24:04 +0000