summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Data/Packed.hs70
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs12
-rw-r--r--lib/Data/Packed/Internal/aux.h2
-rw-r--r--lib/Data/Packed/Matrix.hs8
-rw-r--r--lib/GSL/Matrix.hs35
-rw-r--r--lib/GSL/gsl-aux.h2
-rw-r--r--lib/LinearAlgebra.hs16
-rw-r--r--lib/LinearAlgebra/Algorithms.hs125
-rw-r--r--lib/LinearAlgebra/Interface.hs7
-rw-r--r--lib/LinearAlgebra/Linear.hs43
10 files changed, 169 insertions, 151 deletions
diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs
new file mode 100644
index 0000000..668d2f7
--- /dev/null
+++ b/lib/Data/Packed.hs
@@ -0,0 +1,70 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Data.Packed
5Copyright : (c) Alberto Ruiz 2006-7
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses ffi
11
12The Vector and Matrix types and some utilities.
13
14-}
15-----------------------------------------------------------------------------
16
17module Data.Packed (
18 module Data.Packed.Vector,
19 module Data.Packed.Matrix,
20 module Data.Complex,
21 Container(..)
22) where
23
24import Data.Packed.Vector
25import Data.Packed.Matrix
26import Data.Complex
27import Data.Packed.Internal
28
29-- | conversion utilities
30class (Field e) => Container c e where
31 toComplex :: RealFloat e => (c e, c e) -> c (Complex e)
32 fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
33 comp :: RealFloat e => c e -> c (Complex e)
34 conj :: RealFloat e => c (Complex e) -> c (Complex e)
35 real :: c Double -> c e
36 complex :: c e -> c (Complex Double)
37
38instance Container Vector Double where
39 toComplex = Data.Packed.Internal.toComplex
40 fromComplex = Data.Packed.Internal.fromComplex
41 comp = Data.Packed.Internal.comp
42 conj = Data.Packed.Internal.conj
43 real = id
44 complex = Data.Packed.comp
45
46instance Container Vector (Complex Double) where
47 toComplex = undefined -- can't match
48 fromComplex = undefined
49 comp = undefined
50 conj = undefined
51 real = Data.Packed.comp
52 complex = id
53
54instance Container Matrix Double where
55 toComplex = uncurry $ liftMatrix2 $ curry Data.Packed.toComplex
56 fromComplex z = (reshape c r, reshape c i)
57 where (r,i) = Data.Packed.fromComplex (cdat z)
58 c = cols z
59 comp = liftMatrix Data.Packed.Internal.comp
60 conj = liftMatrix Data.Packed.Internal.conj
61 real = id
62 complex = Data.Packed.comp
63
64instance Container Matrix (Complex Double) where
65 toComplex = undefined
66 fromComplex = undefined
67 comp = undefined
68 conj = undefined
69 real = Data.Packed.comp
70 complex = id
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 63ebddf..e76500b 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -24,6 +24,8 @@ import Complex
24import Control.Monad(when) 24import Control.Monad(when)
25import Data.List(transpose,intersperse) 25import Data.List(transpose,intersperse)
26import Data.Maybe(fromJust) 26import Data.Maybe(fromJust)
27import Foreign.C.String
28import Foreign.C.Types
27 29
28----------------------------------------------------------------- 30-----------------------------------------------------------------
29 31
@@ -371,6 +373,16 @@ fromComplex z = (r,i) where
371comp :: Vector Double -> Vector (Complex Double) 373comp :: Vector Double -> Vector (Complex Double)
372comp v = toComplex (v,constant 0 (dim v)) 374comp v = toComplex (v,constant 0 (dim v))
373 375
376-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance).
377fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
378fromFile filename (r,c) = do
379 charname <- newCString filename
380 res <- createMatrix RowMajor r c
381 c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
382 --free charname -- TO DO: free the auxiliary CString
383 return res
384foreign import ccall "aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
385
374------------------------------------------------------------------------- 386-------------------------------------------------------------------------
375 387
376-- Generic definitions 388-- Generic definitions
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/aux.h
index 83111e5..73334e3 100644
--- a/lib/Data/Packed/Internal/aux.h
+++ b/lib/Data/Packed/Internal/aux.h
@@ -26,3 +26,5 @@ int diagR(KRVEC(d),RMAT(r));
26int diagC(KCVEC(d),CMAT(r)); 26int diagC(KCVEC(d),CMAT(r));
27 27
28const char * gsl_strerror (const int gsl_errno); 28const char * gsl_strerror (const int gsl_errno);
29
30int matrix_fscanf(char*filename, RMAT(a));
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 404fde7..2b93348 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -25,9 +25,10 @@ module Data.Packed.Matrix (
25 fromBlocks, 25 fromBlocks,
26 flipud, fliprl, 26 flipud, fliprl,
27 subMatrix, takeRows, dropRows, takeColumns, dropColumns, 27 subMatrix, takeRows, dropRows, takeColumns, dropColumns,
28 extractRows,
28 ident, diag, diagRect, takeDiag, 29 ident, diag, diagRect, takeDiag,
29 liftMatrix, liftMatrix2, 30 liftMatrix, liftMatrix2,
30 format, dispR, readMatrix, fromArray2D 31 format, dispR, readMatrix, fromFile, fromArray2D
31) where 32) where
32 33
33import Data.Packed.Internal 34import Data.Packed.Internal
@@ -209,3 +210,8 @@ dispC d m = disp m (shfc d)
209-- | creates a matrix from a table of numbers. 210-- | creates a matrix from a table of numbers.
210readMatrix :: String -> Matrix Double 211readMatrix :: String -> Matrix Double
211readMatrix = fromLists . map (map read). map words . filter (not.null) . lines 212readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
213
214-- | rearranges the rows of a matrix according to the order given in a list of integers.
215extractRows :: Field t => [Int] -> Matrix t -> Matrix t
216extractRows l m = fromRows $ extract (toRows $ m) l
217 where extract l is = [l!!i |i<-is]
diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs
index c1bce37..3a54226 100644
--- a/lib/GSL/Matrix.hs
+++ b/lib/GSL/Matrix.hs
@@ -18,17 +18,14 @@ module GSL.Matrix(
18 qr, 18 qr,
19 cholR, -- cholC, 19 cholR, -- cholC,
20 luSolveR, luSolveC, 20 luSolveR, luSolveC,
21 luR, luC, 21 luR, luC
22 fromFile, extractRows
23) where 22) where
24 23
25import Data.Packed.Internal 24import Data.Packed.Internal
26import Data.Packed.Matrix(fromLists,ident,takeDiag) 25import Data.Packed.Matrix(fromLists,ident,takeDiag)
27import GSL.Vector 26import GSL.Vector
28import Foreign 27import Foreign
29import Foreign.C.Types
30import Complex 28import Complex
31import Foreign.C.String
32 29
33{- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/. 30{- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/.
34 31
@@ -257,7 +254,7 @@ p is a permutation:
257 254
258L \* U obtains a permuted version of the original matrix: 255L \* U obtains a permuted version of the original matrix:
259 256
260@\> 'extractRows' p m 257@\> extractRows p m
261 2.+3.i -7. 0. 258 2.+3.i -7. 0.
262 1. 2. -3. 259 1. 2. -3.
263 1. -1.i 2.i 260 1. -1.i 2.i
@@ -298,35 +295,7 @@ luC m = (l,u,p, fromIntegral s') where
298 add = liftMatrix2 $ vectorZipC Add 295 add = liftMatrix2 $ vectorZipC Add
299 mul = liftMatrix2 $ vectorZipC Mul 296 mul = liftMatrix2 $ vectorZipC Mul
300 297
301extract l is = [l!!i |i<-is]
302
303{- auxiliary function to get triangular matrices 298{- auxiliary function to get triangular matrices
304-} 299-}
305triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]] 300triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]
306 where el i j = if j-i>=h then v else 1 - v 301 where el i j = if j-i>=h then v else 1 - v
307
308{- | rearranges the rows of a matrix according to the order given in a list of integers.
309
310> > extractRows [3,3,0,1] (ident 4)
311> 0. 0. 0. 1.
312> 0. 0. 0. 1.
313> 1. 0. 0. 0.
314> 0. 1. 0. 0.
315
316-}
317extractRows :: Field t => [Int] -> Matrix t -> Matrix t
318extractRows l m = fromRows $ extract (toRows $ m) l
319
320--------------------------------------------------------------
321
322-- | loads a matrix efficiently from formatted ASCII text file (the number of rows and columns must be known in advance).
323fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
324fromFile filename (r,c) = do
325 charname <- newCString filename
326 res <- createMatrix RowMajor r c
327 c_gslReadMatrix charname // mat dat res // check "gslReadMatrix" []
328 --free charname -- TO DO: free the auxiliary CString
329 return res
330foreign import ccall "gsl-aux.h matrix_fscanf" c_gslReadMatrix:: Ptr CChar -> TM
331
332---------------------------------------------------------------------------
diff --git a/lib/GSL/gsl-aux.h b/lib/GSL/gsl-aux.h
index dd1a247..3ccac25 100644
--- a/lib/GSL/gsl-aux.h
+++ b/lib/GSL/gsl-aux.h
@@ -52,8 +52,6 @@ int integrate_qags(double f(double,void*), double a, double b, double prec, int
52 52
53int polySolve(KRVEC(a), CVEC(z)); 53int polySolve(KRVEC(a), CVEC(z));
54 54
55int matrix_fscanf(char*filename, RMAT(a));
56
57int minimize(double f(int, double*), double tolsize, int maxit, 55int minimize(double f(int, double*), double tolsize, int maxit,
58 KRVEC(xi), KRVEC(sz), RMAT(sol)); 56 KRVEC(xi), KRVEC(sz), RMAT(sol));
59 57
diff --git a/lib/LinearAlgebra.hs b/lib/LinearAlgebra.hs
index d8f44b4..d6ef647 100644
--- a/lib/LinearAlgebra.hs
+++ b/lib/LinearAlgebra.hs
@@ -13,17 +13,15 @@ Basic matrix computations implemented by BLAS, LAPACK and GSL.
13-} 13-}
14----------------------------------------------------------------------------- 14-----------------------------------------------------------------------------
15module LinearAlgebra ( 15module LinearAlgebra (
16 module Data.Packed.Vector, 16 module Data.Packed,
17 module Data.Packed.Matrix, 17 module LinearAlgebra.Linear,
18 module LinearAlgebra.Instances,
19 module LinearAlgebra.Interface,
20 module LinearAlgebra.Algorithms, 18 module LinearAlgebra.Algorithms,
21 module Complex 19 module LinearAlgebra.Instances,
20 module LinearAlgebra.Interface
22) where 21) where
23 22
23import Data.Packed
24import LinearAlgebra.Linear
25import LinearAlgebra.Algorithms
24import LinearAlgebra.Instances 26import LinearAlgebra.Instances
25import LinearAlgebra.Interface 27import LinearAlgebra.Interface
26import LinearAlgebra.Algorithms
27import Data.Packed.Matrix
28import Data.Packed.Vector
29import Complex \ No newline at end of file
diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index 162426f..a67f822 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -9,106 +9,114 @@ Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional 9Stability : provisional
10Portability : uses ffi 10Portability : uses ffi
11 11
12A simple interface to the available matrix computations. We have defined some generic functions 12A generic interface for a number of essential functions. Using it some higher level algorithms
13on both real and complex matrices in such a way that higher level algorithms and 13and testing properties can be written for both real and complex matrices.
14testing properties are valid for both base types.
15 14
16In any case the specific functions for particular base types can also be explicitly 15In any case, the specific functions for particular base types can also be explicitly
17imported from the LAPACK and GSL.Matrix modules. 16imported from the LAPACK and GSL.Matrix modules.
18 17
19-} 18-}
20----------------------------------------------------------------------------- 19-----------------------------------------------------------------------------
21 20
22module LinearAlgebra.Algorithms ( 21module LinearAlgebra.Algorithms (
23-- * Basic Linear Algebra 22-- * Linear Systems
24 scale, 23 linearSolve,
25 add, 24 inv, pinv,
26 multiply, dot, outer, 25 pinvTol, det,
27 linearSolve, linearSolveSVD,
28-- * Matrix factorizations 26-- * Matrix factorizations
29 svd, lu, eig, eigSH, 27-- ** Singular value decomposition
30 qr, chol, 28 svd,
31-- * Utilities 29 full, economy,
32 Normed(..), NormType(..), 30-- ** Eigensystems
33 det,inv,pinv,full,economy, 31 eig, LinearAlgebra.Algorithms.eigS, LinearAlgebra.Algorithms.eigH,
34 pinvTol, 32-- ** Other
35-- pinvTolg, 33 LinearAlgebra.Algorithms.qr, chol,
34-- * Nullspace
36 nullspacePrec, 35 nullspacePrec,
37 nullVector, 36 nullVector,
38-- * Conversions
39 toComplex, fromComplex,
40 comp, real, complex,
41 conj, ctrans,
42-- * Misc 37-- * Misc
43 eps, i, 38 eps, i,
44 scaleRecip, 39 ctrans,
45 addConstant, 40 Normed(..), NormType(..),
46 sub, 41 GenMat(linearSolveSVD,lu,eigSH)
47 mul,
48 divide,
49) where 42) where
50 43
51 44
52import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj) 45import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
53import Data.Packed.Matrix 46import Data.Packed
54import GSL.Matrix 47import GSL.Matrix
55import GSL.Vector 48import GSL.Vector
56import LAPACK 49import LAPACK
57import Complex 50import Complex
58import LinearAlgebra.Linear 51import LinearAlgebra.Linear
59 52
60-- | Base types for which some optimized matrix computations are available 53-- | matrix computations available for both real and complex matrices
61class (Linear Matrix t) => Optimized t where 54class (Linear Matrix t) => GenMat t where
62 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t) 55 svd :: Matrix t -> (Matrix t, Vector Double, Matrix t)
63 lu :: Matrix t -> (Matrix t, Matrix t, [Int], t) 56 lu :: Matrix t -> (Matrix t, Matrix t, [Int], t)
64 linearSolve :: Matrix t -> Matrix t -> Matrix t 57 linearSolve :: Matrix t -> Matrix t -> Matrix t
65 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t 58 linearSolveSVD :: Matrix t -> Matrix t -> Matrix t
66 ctrans :: Matrix t -> Matrix t
67 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 59 eig :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
68 eigSH :: Matrix t -> (Vector Double, Matrix t) 60 eigSH :: Matrix t -> (Vector Double, Matrix t)
69 chol :: Matrix t -> Matrix t 61 chol :: Matrix t -> Matrix t
62 -- | conjugate transpose
63 ctrans :: Matrix t -> Matrix t
70 64
71instance Optimized Double where 65instance GenMat Double where
72 svd = svdR 66 svd = svdR
73 lu = luR 67 lu = luR
74 linearSolve = linearSolveR 68 linearSolve = linearSolveR
75 linearSolveSVD = linearSolveSVDR Nothing 69 linearSolveSVD = linearSolveSVDR Nothing
76 ctrans = trans 70 ctrans = trans
77 eig = eigR 71 eig = eigR
78 eigSH = eigS 72 eigSH = LAPACK.eigS
79 chol = cholR 73 chol = cholR
80 74
81instance Optimized (Complex Double) where 75instance GenMat (Complex Double) where
82 svd = svdC 76 svd = svdC
83 lu = luC 77 lu = luC
84 linearSolve = linearSolveC 78 linearSolve = linearSolveC
85 linearSolveSVD = linearSolveSVDC Nothing 79 linearSolveSVD = linearSolveSVDC Nothing
86 ctrans = conjTrans 80 ctrans = conjTrans
87 eig = eigC 81 eig = eigC
88 eigSH = eigH 82 eigSH = LAPACK.eigH
89 chol = error "cholC not yet implemented" -- waiting for GSL-1.10 83 chol = error "cholC not yet implemented" -- waiting for GSL-1.10
90 84
85-- | eigensystem of a symmetric matrix
86eigS :: Matrix Double -> (Vector Double, Matrix Double)
87eigS = LAPACK.eigS
88
89-- | eigensystem of a hermitian matrix
90eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
91eigH = LAPACK.eigH
92
93qr :: Matrix Double -> (Matrix Double, Matrix Double)
94qr = GSL.Matrix.qr
95
91square m = rows m == cols m 96square m = rows m == cols m
92 97
93det :: Optimized t => Matrix t -> t 98det :: GenMat t => Matrix t -> t
94det m | square m = s * (product $ toList $ takeDiag $ u) 99det m | square m = s * (product $ toList $ takeDiag $ u)
95 | otherwise = error "det of nonsquare matrix" 100 | otherwise = error "det of nonsquare matrix"
96 where (_,u,_,s) = lu m 101 where (_,u,_,s) = lu m
97 102
98inv :: Optimized t => Matrix t -> Matrix t 103inv :: GenMat t => Matrix t -> Matrix t
99inv m | square m = m `linearSolve` ident (rows m) 104inv m | square m = m `linearSolve` ident (rows m)
100 | otherwise = error "inv of nonsquare matrix" 105 | otherwise = error "inv of nonsquare matrix"
101 106
102pinv :: Optimized t => Matrix t -> Matrix t 107pinv :: GenMat t => Matrix t -> Matrix t
103pinv m = linearSolveSVD m (ident (rows m)) 108pinv m = linearSolveSVD m (ident (rows m))
104 109
105 110full :: Field t
111 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t)
106full svd m = (u, d ,v) where 112full svd m = (u, d ,v) where
107 (u,s,v) = svd m 113 (u,s,v) = svd m
108 d = diagRect s r c 114 d = diagRect s r c
109 r = rows m 115 r = rows m
110 c = cols m 116 c = cols m
111 117
118economy :: Field t
119 => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t)
112economy svd m = (u', subVector 0 d s, v') where 120economy svd m = (u', subVector 0 d s, v') where
113 (u,s,v) = svd m 121 (u,s,v) = svd m
114 sl@(g:_) = toList (complex s) 122 sl@(g:_) = toList (complex s)
@@ -123,39 +131,25 @@ economy svd m = (u', subVector 0 d s, v') where
123 v' = takeColumns d v 131 v' = takeColumns d v
124 132
125 133
126{- | Machine precision of a Double. 134-- | The machine precision of a Double: @eps == 2.22044604925031e-16@ (the value used by GNU-Octave).
127
128>> eps
129> 2.22044604925031e-16
130
131(The value used by GNU-Octave)
132
133-}
134eps :: Double 135eps :: Double
135eps = 2.22044604925031e-16 136eps = 2.22044604925031e-16
136 137
137{- | The imaginary unit 138-- | The imaginary unit: @i == 0.0 :+ 1.0@
138
139@> 'ident' 3 \<\> i
1401.i 0. 0.
141 0. 1.i 0.
142 0. 0. 1.i@
143
144-}
145i :: Complex Double 139i :: Complex Double
146i = 0:+1 140i = 0:+1
147 141
148 142
149-- | matrix product 143-- | matrix product
150mXm :: (Num t, Optimized t) => Matrix t -> Matrix t -> Matrix t 144mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t
151mXm = multiply 145mXm = multiply
152 146
153-- | matrix - vector product 147-- | matrix - vector product
154mXv :: (Num t, Optimized t) => Matrix t -> Vector t -> Vector t 148mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t
155mXv m v = flatten $ m `mXm` (asColumn v) 149mXv m v = flatten $ m `mXm` (asColumn v)
156 150
157-- | vector - matrix product 151-- | vector - matrix product
158vXm :: (Num t, Optimized t) => Vector t -> Matrix t -> Vector t 152vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t
159vXm v m = flatten $ (asRow v) `mXm` m 153vXm v m = flatten $ (asRow v) `mXm` m
160 154
161 155
@@ -167,15 +161,6 @@ norm2 = toScalarR Norm2
167norm1 :: Vector Double -> Double 161norm1 :: Vector Double -> Double
168norm1 = toScalarR AbsSum 162norm1 = toScalarR AbsSum
169 163
170vectorMax :: Vector Double -> Double
171vectorMax = toScalarR Max
172vectorMin :: Vector Double -> Double
173vectorMin = toScalarR Min
174vectorMaxIndex :: Vector Double -> Int
175vectorMaxIndex = round . toScalarR MaxIdx
176vectorMinIndex :: Vector Double -> Int
177vectorMinIndex = round . toScalarR MinIdx
178
179data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int 164data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int
180 165
181pnormRV PNorm2 = norm2 166pnormRV PNorm2 = norm2
@@ -199,7 +184,7 @@ pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` const
199--pnormCM _ _ = error "p norm not yet defined" 184--pnormCM _ _ = error "p norm not yet defined"
200 185
201-- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'. 186-- -- | computes the p-norm of a matrix or vector (with the same definitions as GNU-octave). pnorm 0 denotes \\inf-norm. See also 'norm'.
202--pnorm :: (Container t, Optimized a) => Int -> t a -> Double 187--pnorm :: (Container t, GenMat a) => Int -> t a -> Double
203--pnorm = pnormG 188--pnorm = pnormG
204 189
205class Normed t where 190class Normed t where
@@ -222,7 +207,7 @@ instance Normed (Matrix (Complex Double)) where
222----------------------------------------------------------------------- 207-----------------------------------------------------------------------
223 208
224-- | The nullspace of a matrix from its SVD decomposition. 209-- | The nullspace of a matrix from its SVD decomposition.
225nullspacePrec :: Optimized t 210nullspacePrec :: GenMat t
226 => Double -- ^ relative tolerance in 'eps' units 211 => Double -- ^ relative tolerance in 'eps' units
227 -> Matrix t -- ^ input matrix 212 -> Matrix t -- ^ input matrix
228 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace 213 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
@@ -235,12 +220,12 @@ nullspacePrec t m = ns where
235 ns = drop rank $ toRows $ ctrans v 220 ns = drop rank $ toRows $ ctrans v
236 221
237-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@). 222-- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance (shortcut for @last . nullspacePrec 1@).
238nullVector :: Optimized t => Matrix t -> Vector t 223nullVector :: GenMat t => Matrix t -> Vector t
239nullVector = last . nullspacePrec 1 224nullVector = last . nullspacePrec 1
240 225
241------------------------------------------------------------------------ 226------------------------------------------------------------------------
242 227
243{- | Pseudoinverse of a real matrix with the desired tolerance, expressed as a 228{- Pseudoinverse of a real matrix with the desired tolerance, expressed as a
244multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv'). 229multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
245 230
246@\> let m = 'fromLists' [[1,0, 0] 231@\> let m = 'fromLists' [[1,0, 0]
@@ -258,7 +243,7 @@ multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv').
2580. 0. 1.@ 2430. 0. 1.@
259 244
260-} 245-}
261pinvTol :: Double -> Matrix Double -> Matrix Double 246--pinvTol :: Double -> Matrix Double -> Matrix Double
262pinvTol t m = v' `mXm` diag s' `mXm` trans u' where 247pinvTol t m = v' `mXm` diag s' `mXm` trans u' where
263 (u,s,v) = svdR m 248 (u,s,v) = svdR m
264 sl@(g:_) = toList s 249 sl@(g:_) = toList s
diff --git a/lib/LinearAlgebra/Interface.hs b/lib/LinearAlgebra/Interface.hs
index 3392d54..0c65a8b 100644
--- a/lib/LinearAlgebra/Interface.hs
+++ b/lib/LinearAlgebra/Interface.hs
@@ -16,6 +16,7 @@ Operators for frequent operations.
16 16
17module LinearAlgebra.Interface( 17module LinearAlgebra.Interface(
18 (<>),(<.>), 18 (<>),(<.>),
19 (<\>),
19 (.*),(*/), 20 (.*),(*/),
20 (<|>),(<->), 21 (<|>),(<->),
21) where 22) where
@@ -23,6 +24,7 @@ module LinearAlgebra.Interface(
23import LinearAlgebra.Linear 24import LinearAlgebra.Linear
24import Data.Packed.Vector 25import Data.Packed.Vector
25import Data.Packed.Matrix 26import Data.Packed.Matrix
27import LinearAlgebra.Algorithms
26 28
27class Mul a b c | a b -> c where 29class Mul a b c | a b -> c where
28 infixl 7 <> 30 infixl 7 <>
@@ -59,6 +61,11 @@ a .* x = scale a x
59infixl 7 */ 61infixl 7 */
60v */ x = scale (recip x) v 62v */ x = scale (recip x) v
61 63
64-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
65(<\>) :: (GenMat a) => Matrix a -> Vector a -> Vector a
66infixl 7 <\>
67m <\> v = flatten (linearSolveSVD m (reshape 1 v))
68
62------------------------------------------------ 69------------------------------------------------
63 70
64class Joinable a b where 71class Joinable a b where
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index 85daa4a..3e6c55d 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -9,10 +9,10 @@ Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional 9Stability : provisional
10Portability : uses ffi 10Portability : uses ffi
11 11
12Basic optimized operations on vectors and matrices.
12 13
13-} 14-}
14----------------------------------------------------------------------------- 15-----------------------------------------------------------------------------
15-- #hide
16 16
17module LinearAlgebra.Linear ( 17module LinearAlgebra.Linear (
18 Linear(..), 18 Linear(..),
@@ -21,25 +21,22 @@ module LinearAlgebra.Linear (
21 21
22 22
23import Data.Packed.Internal 23import Data.Packed.Internal
24import Data.Packed.Matrix 24import Data.Packed
25import GSL.Vector 25import GSL.Vector
26import Complex 26import Complex
27 27
28-- | basic optimized operations 28-- | basic optimized operations
29class (Field e) => Linear c e where 29class (Container c e) => Linear c e where
30 scale :: e -> c e -> c e 30 scale :: e -> c e -> c e
31 scaleRecip :: e -> c e -> c e
32 addConstant :: e -> c e -> c e 31 addConstant :: e -> c e -> c e
33 add :: c e -> c e -> c e 32 add :: c e -> c e -> c e
34 sub :: c e -> c e -> c e 33 sub :: c e -> c e -> c e
34 -- | element by element multiplication
35 mul :: c e -> c e -> c e 35 mul :: c e -> c e -> c e
36 -- | element by element division
36 divide :: c e -> c e -> c e 37 divide :: c e -> c e -> c e
37 toComplex :: RealFloat e => (c e, c e) -> c (Complex e) 38 -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
38 fromComplex :: RealFloat e => c (Complex e) -> (c e, c e) 39 scaleRecip :: e -> c e -> c e
39 comp :: RealFloat e => c e -> c (Complex e)
40 conj :: RealFloat e => c (Complex e) -> c (Complex e)
41 real :: c Double -> c e
42 complex :: c e -> c (Complex Double)
43 40
44instance Linear Vector Double where 41instance Linear Vector Double where
45 scale = vectorMapValR Scale 42 scale = vectorMapValR Scale
@@ -49,12 +46,6 @@ instance Linear Vector Double where
49 sub = vectorZipR Sub 46 sub = vectorZipR Sub
50 mul = vectorZipR Mul 47 mul = vectorZipR Mul
51 divide = vectorZipR Div 48 divide = vectorZipR Div
52 toComplex = Data.Packed.Internal.toComplex
53 fromComplex = Data.Packed.Internal.fromComplex
54 comp = Data.Packed.Internal.comp
55 conj = Data.Packed.Internal.conj
56 real = id
57 complex = LinearAlgebra.Linear.comp
58 49
59instance Linear Vector (Complex Double) where 50instance Linear Vector (Complex Double) where
60 scale = vectorMapValC Scale 51 scale = vectorMapValC Scale
@@ -64,12 +55,6 @@ instance Linear Vector (Complex Double) where
64 sub = vectorZipC Sub 55 sub = vectorZipC Sub
65 mul = vectorZipC Mul 56 mul = vectorZipC Mul
66 divide = vectorZipC Div 57 divide = vectorZipC Div
67 toComplex = undefined -- can't match
68 fromComplex = undefined
69 comp = undefined
70 conj = undefined
71 real = LinearAlgebra.Linear.comp
72 complex = id
73 58
74instance Linear Matrix Double where 59instance Linear Matrix Double where
75 scale x = liftMatrix (scale x) 60 scale x = liftMatrix (scale x)
@@ -79,14 +64,6 @@ instance Linear Matrix Double where
79 sub = liftMatrix2 sub 64 sub = liftMatrix2 sub
80 mul = liftMatrix2 mul 65 mul = liftMatrix2 mul
81 divide = liftMatrix2 divide 66 divide = liftMatrix2 divide
82 toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex
83 fromComplex z = (reshape c r, reshape c i)
84 where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z)
85 c = cols z
86 comp = liftMatrix Data.Packed.Internal.comp
87 conj = liftMatrix Data.Packed.Internal.conj
88 real = id
89 complex = LinearAlgebra.Linear.comp
90 67
91instance Linear Matrix (Complex Double) where 68instance Linear Matrix (Complex Double) where
92 scale x = liftMatrix (scale x) 69 scale x = liftMatrix (scale x)
@@ -96,12 +73,6 @@ instance Linear Matrix (Complex Double) where
96 sub = liftMatrix2 sub 73 sub = liftMatrix2 sub
97 mul = liftMatrix2 mul 74 mul = liftMatrix2 mul
98 divide = liftMatrix2 divide 75 divide = liftMatrix2 divide
99 toComplex = undefined
100 fromComplex = undefined
101 comp = undefined
102 conj = undefined
103 real = LinearAlgebra.Linear.comp
104 complex = id
105 76
106-------------------------------------------------- 77--------------------------------------------------
107 78