diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Data/Packed.hs | 70 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 12 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/aux.h | 2 | ||||
-rw-r--r-- | lib/Data/Packed/Matrix.hs | 8 | ||||
-rw-r--r-- | lib/GSL/Matrix.hs | 35 | ||||
-rw-r--r-- | lib/GSL/gsl-aux.h | 2 | ||||
-rw-r--r-- | lib/LinearAlgebra.hs | 16 | ||||
-rw-r--r-- | lib/LinearAlgebra/Algorithms.hs | 125 | ||||
-rw-r--r-- | lib/LinearAlgebra/Interface.hs | 7 | ||||
-rw-r--r-- | lib/LinearAlgebra/Linear.hs | 43 |
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 | {- | | ||
4 | Module : Data.Packed | ||
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 : uses ffi | ||
11 | |||
12 | The Vector and Matrix types and some utilities. | ||
13 | |||
14 | -} | ||
15 | ----------------------------------------------------------------------------- | ||
16 | |||
17 | module Data.Packed ( | ||
18 | module Data.Packed.Vector, | ||
19 | module Data.Packed.Matrix, | ||
20 | module Data.Complex, | ||
21 | Container(..) | ||
22 | ) where | ||
23 | |||
24 | import Data.Packed.Vector | ||
25 | import Data.Packed.Matrix | ||
26 | import Data.Complex | ||
27 | import Data.Packed.Internal | ||
28 | |||
29 | -- | conversion utilities | ||
30 | class (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 | |||
38 | instance 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 | |||
46 | instance 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 | |||
54 | instance 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 | |||
64 | instance 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 | |||
24 | import Control.Monad(when) | 24 | import Control.Monad(when) |
25 | import Data.List(transpose,intersperse) | 25 | import Data.List(transpose,intersperse) |
26 | import Data.Maybe(fromJust) | 26 | import Data.Maybe(fromJust) |
27 | import Foreign.C.String | ||
28 | import Foreign.C.Types | ||
27 | 29 | ||
28 | ----------------------------------------------------------------- | 30 | ----------------------------------------------------------------- |
29 | 31 | ||
@@ -371,6 +373,16 @@ fromComplex z = (r,i) where | |||
371 | comp :: Vector Double -> Vector (Complex Double) | 373 | comp :: Vector Double -> Vector (Complex Double) |
372 | comp v = toComplex (v,constant 0 (dim v)) | 374 | comp 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). | ||
377 | fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) | ||
378 | fromFile 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 | ||
384 | foreign 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)); | |||
26 | int diagC(KCVEC(d),CMAT(r)); | 26 | int diagC(KCVEC(d),CMAT(r)); |
27 | 27 | ||
28 | const char * gsl_strerror (const int gsl_errno); | 28 | const char * gsl_strerror (const int gsl_errno); |
29 | |||
30 | int 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 | ||
33 | import Data.Packed.Internal | 34 | import 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. |
210 | readMatrix :: String -> Matrix Double | 211 | readMatrix :: String -> Matrix Double |
211 | readMatrix = fromLists . map (map read). map words . filter (not.null) . lines | 212 | readMatrix = 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. | ||
215 | extractRows :: Field t => [Int] -> Matrix t -> Matrix t | ||
216 | extractRows 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 | ||
25 | import Data.Packed.Internal | 24 | import Data.Packed.Internal |
26 | import Data.Packed.Matrix(fromLists,ident,takeDiag) | 25 | import Data.Packed.Matrix(fromLists,ident,takeDiag) |
27 | import GSL.Vector | 26 | import GSL.Vector |
28 | import Foreign | 27 | import Foreign |
29 | import Foreign.C.Types | ||
30 | import Complex | 28 | import Complex |
31 | import 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 | ||
258 | L \* U obtains a permuted version of the original matrix: | 255 | L \* 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 | ||
301 | extract l is = [l!!i |i<-is] | ||
302 | |||
303 | {- auxiliary function to get triangular matrices | 298 | {- auxiliary function to get triangular matrices |
304 | -} | 299 | -} |
305 | triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]] | 300 | triang 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 | -} | ||
317 | extractRows :: Field t => [Int] -> Matrix t -> Matrix t | ||
318 | extractRows 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). | ||
323 | fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) | ||
324 | fromFile 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 | ||
330 | foreign 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 | ||
53 | int polySolve(KRVEC(a), CVEC(z)); | 53 | int polySolve(KRVEC(a), CVEC(z)); |
54 | 54 | ||
55 | int matrix_fscanf(char*filename, RMAT(a)); | ||
56 | |||
57 | int minimize(double f(int, double*), double tolsize, int maxit, | 55 | int 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 | ----------------------------------------------------------------------------- |
15 | module LinearAlgebra ( | 15 | module 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 | ||
23 | import Data.Packed | ||
24 | import LinearAlgebra.Linear | ||
25 | import LinearAlgebra.Algorithms | ||
24 | import LinearAlgebra.Instances | 26 | import LinearAlgebra.Instances |
25 | import LinearAlgebra.Interface | 27 | import LinearAlgebra.Interface |
26 | import LinearAlgebra.Algorithms | ||
27 | import Data.Packed.Matrix | ||
28 | import Data.Packed.Vector | ||
29 | import 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) | |||
9 | Stability : provisional | 9 | Stability : provisional |
10 | Portability : uses ffi | 10 | Portability : uses ffi |
11 | 11 | ||
12 | A simple interface to the available matrix computations. We have defined some generic functions | 12 | A generic interface for a number of essential functions. Using it some higher level algorithms |
13 | on both real and complex matrices in such a way that higher level algorithms and | 13 | and testing properties can be written for both real and complex matrices. |
14 | testing properties are valid for both base types. | ||
15 | 14 | ||
16 | In any case the specific functions for particular base types can also be explicitly | 15 | In any case, the specific functions for particular base types can also be explicitly |
17 | imported from the LAPACK and GSL.Matrix modules. | 16 | imported from the LAPACK and GSL.Matrix modules. |
18 | 17 | ||
19 | -} | 18 | -} |
20 | ----------------------------------------------------------------------------- | 19 | ----------------------------------------------------------------------------- |
21 | 20 | ||
22 | module LinearAlgebra.Algorithms ( | 21 | module 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 | ||
52 | import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj) | 45 | import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj) |
53 | import Data.Packed.Matrix | 46 | import Data.Packed |
54 | import GSL.Matrix | 47 | import GSL.Matrix |
55 | import GSL.Vector | 48 | import GSL.Vector |
56 | import LAPACK | 49 | import LAPACK |
57 | import Complex | 50 | import Complex |
58 | import LinearAlgebra.Linear | 51 | import 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 |
61 | class (Linear Matrix t) => Optimized t where | 54 | class (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 | ||
71 | instance Optimized Double where | 65 | instance 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 | ||
81 | instance Optimized (Complex Double) where | 75 | instance 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 | ||
86 | eigS :: Matrix Double -> (Vector Double, Matrix Double) | ||
87 | eigS = LAPACK.eigS | ||
88 | |||
89 | -- | eigensystem of a hermitian matrix | ||
90 | eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
91 | eigH = LAPACK.eigH | ||
92 | |||
93 | qr :: Matrix Double -> (Matrix Double, Matrix Double) | ||
94 | qr = GSL.Matrix.qr | ||
95 | |||
91 | square m = rows m == cols m | 96 | square m = rows m == cols m |
92 | 97 | ||
93 | det :: Optimized t => Matrix t -> t | 98 | det :: GenMat t => Matrix t -> t |
94 | det m | square m = s * (product $ toList $ takeDiag $ u) | 99 | det 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 | ||
98 | inv :: Optimized t => Matrix t -> Matrix t | 103 | inv :: GenMat t => Matrix t -> Matrix t |
99 | inv m | square m = m `linearSolve` ident (rows m) | 104 | inv m | square m = m `linearSolve` ident (rows m) |
100 | | otherwise = error "inv of nonsquare matrix" | 105 | | otherwise = error "inv of nonsquare matrix" |
101 | 106 | ||
102 | pinv :: Optimized t => Matrix t -> Matrix t | 107 | pinv :: GenMat t => Matrix t -> Matrix t |
103 | pinv m = linearSolveSVD m (ident (rows m)) | 108 | pinv m = linearSolveSVD m (ident (rows m)) |
104 | 109 | ||
105 | 110 | full :: Field t | |
111 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Matrix Double, Matrix t) | ||
106 | full svd m = (u, d ,v) where | 112 | full 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 | ||
118 | economy :: Field t | ||
119 | => (Matrix t -> (Matrix t, Vector Double, Matrix t)) -> Matrix t -> (Matrix t, Vector Double, Matrix t) | ||
112 | economy svd m = (u', subVector 0 d s, v') where | 120 | economy 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 | -} | ||
134 | eps :: Double | 135 | eps :: Double |
135 | eps = 2.22044604925031e-16 | 136 | eps = 2.22044604925031e-16 |
136 | 137 | ||
137 | {- | The imaginary unit | 138 | -- | The imaginary unit: @i == 0.0 :+ 1.0@ |
138 | |||
139 | @> 'ident' 3 \<\> i | ||
140 | 1.i 0. 0. | ||
141 | 0. 1.i 0. | ||
142 | 0. 0. 1.i@ | ||
143 | |||
144 | -} | ||
145 | i :: Complex Double | 139 | i :: Complex Double |
146 | i = 0:+1 | 140 | i = 0:+1 |
147 | 141 | ||
148 | 142 | ||
149 | -- | matrix product | 143 | -- | matrix product |
150 | mXm :: (Num t, Optimized t) => Matrix t -> Matrix t -> Matrix t | 144 | mXm :: (Num t, GenMat t) => Matrix t -> Matrix t -> Matrix t |
151 | mXm = multiply | 145 | mXm = multiply |
152 | 146 | ||
153 | -- | matrix - vector product | 147 | -- | matrix - vector product |
154 | mXv :: (Num t, Optimized t) => Matrix t -> Vector t -> Vector t | 148 | mXv :: (Num t, GenMat t) => Matrix t -> Vector t -> Vector t |
155 | mXv m v = flatten $ m `mXm` (asColumn v) | 149 | mXv m v = flatten $ m `mXm` (asColumn v) |
156 | 150 | ||
157 | -- | vector - matrix product | 151 | -- | vector - matrix product |
158 | vXm :: (Num t, Optimized t) => Vector t -> Matrix t -> Vector t | 152 | vXm :: (Num t, GenMat t) => Vector t -> Matrix t -> Vector t |
159 | vXm v m = flatten $ (asRow v) `mXm` m | 153 | vXm v m = flatten $ (asRow v) `mXm` m |
160 | 154 | ||
161 | 155 | ||
@@ -167,15 +161,6 @@ norm2 = toScalarR Norm2 | |||
167 | norm1 :: Vector Double -> Double | 161 | norm1 :: Vector Double -> Double |
168 | norm1 = toScalarR AbsSum | 162 | norm1 = toScalarR AbsSum |
169 | 163 | ||
170 | vectorMax :: Vector Double -> Double | ||
171 | vectorMax = toScalarR Max | ||
172 | vectorMin :: Vector Double -> Double | ||
173 | vectorMin = toScalarR Min | ||
174 | vectorMaxIndex :: Vector Double -> Int | ||
175 | vectorMaxIndex = round . toScalarR MaxIdx | ||
176 | vectorMinIndex :: Vector Double -> Int | ||
177 | vectorMinIndex = round . toScalarR MinIdx | ||
178 | |||
179 | data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int | 164 | data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int |
180 | 165 | ||
181 | pnormRV PNorm2 = norm2 | 166 | pnormRV 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 | ||
205 | class Normed t where | 190 | class 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. |
225 | nullspacePrec :: Optimized t | 210 | nullspacePrec :: 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@). |
238 | nullVector :: Optimized t => Matrix t -> Vector t | 223 | nullVector :: GenMat t => Matrix t -> Vector t |
239 | nullVector = last . nullspacePrec 1 | 224 | nullVector = 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 |
244 | multiplicative factor of the default tolerance used by GNU-Octave (see 'pinv'). | 229 | multiplicative 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'). | |||
258 | 0. 0. 1.@ | 243 | 0. 0. 1.@ |
259 | 244 | ||
260 | -} | 245 | -} |
261 | pinvTol :: Double -> Matrix Double -> Matrix Double | 246 | --pinvTol :: Double -> Matrix Double -> Matrix Double |
262 | pinvTol t m = v' `mXm` diag s' `mXm` trans u' where | 247 | pinvTol 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 | ||
17 | module LinearAlgebra.Interface( | 17 | module LinearAlgebra.Interface( |
18 | (<>),(<.>), | 18 | (<>),(<.>), |
19 | (<\>), | ||
19 | (.*),(*/), | 20 | (.*),(*/), |
20 | (<|>),(<->), | 21 | (<|>),(<->), |
21 | ) where | 22 | ) where |
@@ -23,6 +24,7 @@ module LinearAlgebra.Interface( | |||
23 | import LinearAlgebra.Linear | 24 | import LinearAlgebra.Linear |
24 | import Data.Packed.Vector | 25 | import Data.Packed.Vector |
25 | import Data.Packed.Matrix | 26 | import Data.Packed.Matrix |
27 | import LinearAlgebra.Algorithms | ||
26 | 28 | ||
27 | class Mul a b c | a b -> c where | 29 | class Mul a b c | a b -> c where |
28 | infixl 7 <> | 30 | infixl 7 <> |
@@ -59,6 +61,11 @@ a .* x = scale a x | |||
59 | infixl 7 */ | 61 | infixl 7 */ |
60 | v */ x = scale (recip x) v | 62 | v */ 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 | ||
66 | infixl 7 <\> | ||
67 | m <\> v = flatten (linearSolveSVD m (reshape 1 v)) | ||
68 | |||
62 | ------------------------------------------------ | 69 | ------------------------------------------------ |
63 | 70 | ||
64 | class Joinable a b where | 71 | class 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) | |||
9 | Stability : provisional | 9 | Stability : provisional |
10 | Portability : uses ffi | 10 | Portability : uses ffi |
11 | 11 | ||
12 | Basic optimized operations on vectors and matrices. | ||
12 | 13 | ||
13 | -} | 14 | -} |
14 | ----------------------------------------------------------------------------- | 15 | ----------------------------------------------------------------------------- |
15 | -- #hide | ||
16 | 16 | ||
17 | module LinearAlgebra.Linear ( | 17 | module LinearAlgebra.Linear ( |
18 | Linear(..), | 18 | Linear(..), |
@@ -21,25 +21,22 @@ module LinearAlgebra.Linear ( | |||
21 | 21 | ||
22 | 22 | ||
23 | import Data.Packed.Internal | 23 | import Data.Packed.Internal |
24 | import Data.Packed.Matrix | 24 | import Data.Packed |
25 | import GSL.Vector | 25 | import GSL.Vector |
26 | import Complex | 26 | import Complex |
27 | 27 | ||
28 | -- | basic optimized operations | 28 | -- | basic optimized operations |
29 | class (Field e) => Linear c e where | 29 | class (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 | ||
44 | instance Linear Vector Double where | 41 | instance 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 | ||
59 | instance Linear Vector (Complex Double) where | 50 | instance 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 | ||
74 | instance Linear Matrix Double where | 59 | instance 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 | ||
91 | instance Linear Matrix (Complex Double) where | 68 | instance 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 | ||