summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/LAPACK.hs261
-rw-r--r--lib/LAPACK/Internal.hs262
2 files changed, 259 insertions, 264 deletions
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-09-22 10:27:21 +0000