summaryrefslogtreecommitdiff
path: root/lib/LAPACK
diff options
context:
space:
mode:
Diffstat (limited to 'lib/LAPACK')
-rw-r--r--lib/LAPACK/Internal.hs53
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})
174foreign import ccall "lapack-aux.h linearSolveR_l" 174foreign 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.
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
177----------------------------------------------------------------------------- 186-----------------------------------------------------------------------------
178-- zgesv 187-- zgesv
179foreign import ccall "lapack-aux.h linearSolveC_l" 188foreign 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.
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
182----------------------------------------------------------------------------------- 200-----------------------------------------------------------------------------------
183-- dgels 201-- dgels
184foreign import ccall "lapack-aux.h linearSolveLSR_l" 202foreign import ccall "lapack-aux.h linearSolveLSR_l"
@@ -198,12 +216,47 @@ linearSolveLSR_l a@(M {rows = m, cols = n}) b@(M {cols = nrhs}) = unsafePerformI
198foreign import ccall "lapack-aux.h linearSolveLSC_l" 216foreign 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'.
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
201----------------------------------------------------------------------------------- 228-----------------------------------------------------------------------------------
202-- dgelss 229-- dgelss
203foreign import ccall "lapack-aux.h linearSolveSVDR_l" 230foreign 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.
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
206----------------------------------------------------------------------------------- 246-----------------------------------------------------------------------------------
207-- zgelss 247-- zgelss
208foreign import ccall "lapack-aux.h linearSolveSVDC_l" 248foreign 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.
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