diff options
Diffstat (limited to 'lib/LAPACK/Internal.hs')
-rw-r--r-- | lib/LAPACK/Internal.hs | 53 |
1 files changed, 53 insertions, 0 deletions
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs index ba50e6b..ec46b66 100644 --- a/lib/LAPACK/Internal.hs +++ b/lib/LAPACK/Internal.hs | |||
@@ -174,11 +174,29 @@ eigH' (m@M {rows = r}) | |||
174 | foreign import ccall "lapack-aux.h linearSolveR_l" | 174 | foreign import ccall "lapack-aux.h linearSolveR_l" |
175 | dgesv :: Double ::> Double ::> Double ::> IO Int | 175 | dgesv :: Double ::> Double ::> Double ::> IO Int |
176 | 176 | ||
177 | -- | Wrapper for LAPACK's /dgesv/, which solves a general real linear system (for several right-hand sides) internally using the lu decomposition. | ||
178 | linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
179 | linearSolveR 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 | |||
177 | ----------------------------------------------------------------------------- | 186 | ----------------------------------------------------------------------------- |
178 | -- zgesv | 187 | -- zgesv |
179 | foreign import ccall "lapack-aux.h linearSolveC_l" | 188 | foreign import ccall "lapack-aux.h linearSolveC_l" |
180 | zgesv :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int | 189 | zgesv :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int |
181 | 190 | ||
191 | -- | Wrapper for LAPACK's /zgesv/, which solves a general complex linear system (for several right-hand sides) internally using the lu decomposition. | ||
192 | linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
193 | linearSolveC 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 | |||
182 | ----------------------------------------------------------------------------------- | 200 | ----------------------------------------------------------------------------------- |
183 | -- dgels | 201 | -- dgels |
184 | foreign import ccall "lapack-aux.h linearSolveLSR_l" | 202 | foreign import ccall "lapack-aux.h linearSolveLSR_l" |
@@ -198,12 +216,47 @@ linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformI | |||
198 | foreign import ccall "lapack-aux.h linearSolveLSC_l" | 216 | foreign import ccall "lapack-aux.h linearSolveLSC_l" |
199 | zgels :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int | 217 | zgels :: (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int |
200 | 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'. | ||
220 | linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
221 | linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ linearSolveLSC_l a b | ||
222 | |||
223 | linearSolveLSC_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 | |||
201 | ----------------------------------------------------------------------------------- | 228 | ----------------------------------------------------------------------------------- |
202 | -- dgelss | 229 | -- dgelss |
203 | foreign import ccall "lapack-aux.h linearSolveSVDR_l" | 230 | foreign import ccall "lapack-aux.h linearSolveSVDR_l" |
204 | dgelss :: Double -> Double ::> Double ::> Double ::> IO Int | 231 | dgelss :: Double -> Double ::> Double ::> Double ::> IO Int |
205 | 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. | ||
234 | linearSolveSVDR :: Maybe Double -- ^ rcond | ||
235 | -> Matrix Double -- ^ coefficient matrix | ||
236 | -> Matrix Double -- ^ right hand sides (as columns) | ||
237 | -> Matrix Double -- ^ solution vectors (as columns) | ||
238 | linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDR_l rcond a b | ||
239 | linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) a b | ||
240 | |||
241 | linearSolveSVDR_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 | |||
206 | ----------------------------------------------------------------------------------- | 246 | ----------------------------------------------------------------------------------- |
207 | -- zgelss | 247 | -- zgelss |
208 | foreign import ccall "lapack-aux.h linearSolveSVDC_l" | 248 | foreign import ccall "lapack-aux.h linearSolveSVDC_l" |
209 | zgelss :: Double -> (Complex Double) ::> (Complex Double) ::> (Complex Double) ::> IO Int | 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. | ||
252 | linearSolveSVDC :: 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) | ||
256 | linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ linearSolveSVDC_l rcond a b | ||
257 | linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) a b | ||
258 | |||
259 | linearSolveSVDC_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 | ||