summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric/LinearAlgebra
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Numeric/LinearAlgebra')
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs746
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs555
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c1489
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h60
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs295
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs114
6 files changed, 3259 insertions, 0 deletions
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs
new file mode 100644
index 0000000..8c4b610
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -0,0 +1,746 @@
1{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
2{-# LANGUAGE CPP #-}
3{-# LANGUAGE MultiParamTypeClasses #-}
4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE TypeFamilies #-}
6
7-----------------------------------------------------------------------------
8{- |
9Module : Numeric.LinearAlgebra.Algorithms
10Copyright : (c) Alberto Ruiz 2006-9
11License : GPL-style
12
13Maintainer : Alberto Ruiz (aruiz at um dot es)
14Stability : provisional
15Portability : uses ffi
16
17High level generic interface to common matrix computations.
18
19Specific functions for particular base types can also be explicitly
20imported from "Numeric.LinearAlgebra.LAPACK".
21
22-}
23-----------------------------------------------------------------------------
24{-# OPTIONS_HADDOCK hide #-}
25
26
27module Numeric.LinearAlgebra.Algorithms (
28-- * Supported types
29 Field(),
30-- * Linear Systems
31 linearSolve,
32 luSolve,
33 cholSolve,
34 linearSolveLS,
35 linearSolveSVD,
36 inv, pinv, pinvTol,
37 det, invlndet,
38 rank, rcond,
39-- * Matrix factorizations
40-- ** Singular value decomposition
41 svd,
42 fullSVD,
43 thinSVD,
44 compactSVD,
45 singularValues,
46 leftSV, rightSV,
47-- ** Eigensystems
48 eig, eigSH, eigSH',
49 eigenvalues, eigenvaluesSH, eigenvaluesSH',
50 geigSH',
51-- ** QR
52 qr, rq, qrRaw, qrgr,
53-- ** Cholesky
54 chol, cholSH, mbCholSH,
55-- ** Hessenberg
56 hess,
57-- ** Schur
58 schur,
59-- ** LU
60 lu, luPacked,
61-- * Matrix functions
62 expm,
63 sqrtm,
64 matFunc,
65-- * Nullspace
66 nullspacePrec,
67 nullVector,
68 nullspaceSVD,
69 orth,
70-- * Norms
71 Normed(..), NormType(..),
72 relativeError,
73-- * Misc
74 eps, peps, i,
75-- * Util
76 haussholder,
77 unpackQR, unpackHess,
78 ranksv
79) where
80
81
82import Data.Packed.Internal hiding ((//))
83import Data.Packed.Matrix
84import Numeric.LinearAlgebra.LAPACK as LAPACK
85import Data.List(foldl1')
86import Data.Array
87import Numeric.ContainerBoot
88
89
90{- | Generic linear algebra functions for double precision real and complex matrices.
91
92(Single precision data can be converted using 'single' and 'double').
93
94-}
95class (Product t,
96 Convert t,
97 Container Vector t,
98 Container Matrix t,
99 Normed Matrix t,
100 Normed Vector t,
101 Floating t,
102 RealOf t ~ Double) => Field t where
103 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
104 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
105 sv' :: Matrix t -> Vector Double
106 luPacked' :: Matrix t -> (Matrix t, [Int])
107 luSolve' :: (Matrix t, [Int]) -> Matrix t -> Matrix t
108 linearSolve' :: Matrix t -> Matrix t -> Matrix t
109 cholSolve' :: Matrix t -> Matrix t -> Matrix t
110 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
111 linearSolveLS' :: Matrix t -> Matrix t -> Matrix t
112 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
113 eigSH'' :: Matrix t -> (Vector Double, Matrix t)
114 eigOnly :: Matrix t -> Vector (Complex Double)
115 eigOnlySH :: Matrix t -> Vector Double
116 cholSH' :: Matrix t -> Matrix t
117 mbCholSH' :: Matrix t -> Maybe (Matrix t)
118 qr' :: Matrix t -> (Matrix t, Vector t)
119 qrgr' :: Int -> (Matrix t, Vector t) -> Matrix t
120 hess' :: Matrix t -> (Matrix t, Matrix t)
121 schur' :: Matrix t -> (Matrix t, Matrix t)
122
123
124instance Field Double where
125 svd' = svdRd
126 thinSVD' = thinSVDRd
127 sv' = svR
128 luPacked' = luR
129 luSolve' (l_u,perm) = lusR l_u perm
130 linearSolve' = linearSolveR -- (luSolve . luPacked) ??
131 cholSolve' = cholSolveR
132 linearSolveLS' = linearSolveLSR
133 linearSolveSVD' = linearSolveSVDR Nothing
134 eig' = eigR
135 eigSH'' = eigS
136 eigOnly = eigOnlyR
137 eigOnlySH = eigOnlyS
138 cholSH' = cholS
139 mbCholSH' = mbCholS
140 qr' = qrR
141 qrgr' = qrgrR
142 hess' = unpackHess hessR
143 schur' = schurR
144
145instance Field (Complex Double) where
146#ifdef NOZGESDD
147 svd' = svdC
148 thinSVD' = thinSVDC
149#else
150 svd' = svdCd
151 thinSVD' = thinSVDCd
152#endif
153 sv' = svC
154 luPacked' = luC
155 luSolve' (l_u,perm) = lusC l_u perm
156 linearSolve' = linearSolveC
157 cholSolve' = cholSolveC
158 linearSolveLS' = linearSolveLSC
159 linearSolveSVD' = linearSolveSVDC Nothing
160 eig' = eigC
161 eigOnly = eigOnlyC
162 eigSH'' = eigH
163 eigOnlySH = eigOnlyH
164 cholSH' = cholH
165 mbCholSH' = mbCholH
166 qr' = qrC
167 qrgr' = qrgrC
168 hess' = unpackHess hessC
169 schur' = schurC
170
171--------------------------------------------------------------
172
173square m = rows m == cols m
174
175vertical m = rows m >= cols m
176
177exactHermitian m = m `equal` ctrans m
178
179--------------------------------------------------------------
180
181-- | Full singular value decomposition.
182svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
183svd = {-# SCC "svd" #-} svd'
184
185-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
186--
187-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.
188thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
189thinSVD = {-# SCC "thinSVD" #-} thinSVD'
190
191-- | Singular values only.
192singularValues :: Field t => Matrix t -> Vector Double
193singularValues = {-# SCC "singularValues" #-} sv'
194
195-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
196--
197-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
198fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
199fullSVD m = (u,d,v) where
200 (u,s,v) = svd m
201 d = diagRect 0 s r c
202 r = rows m
203 c = cols m
204
205-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
206compactSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
207compactSVD m = (u', subVector 0 d s, v') where
208 (u,s,v) = thinSVD m
209 d = rankSVD (1*eps) m s `max` 1
210 u' = takeColumns d u
211 v' = takeColumns d v
212
213
214-- | Singular values and all right singular vectors.
215rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
216rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
217 | otherwise = let (_,s,v) = svd m in (s,v)
218
219-- | Singular values and all left singular vectors.
220leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
221leftSV m | vertical m = let (u,s,_) = svd m in (u,s)
222 | otherwise = let (u,s,_) = thinSVD m in (u,s)
223
224
225--------------------------------------------------------------
226
227-- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
228luPacked :: Field t => Matrix t -> (Matrix t, [Int])
229luPacked = {-# SCC "luPacked" #-} luPacked'
230
231-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
232luSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
233luSolve = {-# SCC "luSolve" #-} luSolve'
234
235-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
236-- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
237linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
238linearSolve = {-# SCC "linearSolve" #-} linearSolve'
239
240-- | Solve a symmetric or Hermitian positive definite linear system using a precomputed Cholesky decomposition obtained by 'chol'.
241cholSolve :: Field t => Matrix t -> Matrix t -> Matrix t
242cholSolve = {-# SCC "cholSolve" #-} cholSolve'
243
244-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value.
245linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
246linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'
247
248
249-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
250linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
251linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
252
253--------------------------------------------------------------
254
255-- | Eigenvalues and eigenvectors of a general square matrix.
256--
257-- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
258eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
259eig = {-# SCC "eig" #-} eig'
260
261-- | Eigenvalues of a general square matrix.
262eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
263eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
264
265-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
266eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
267eigSH' = {-# SCC "eigSH'" #-} eigSH''
268
269-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
270eigenvaluesSH' :: Field t => Matrix t -> Vector Double
271eigenvaluesSH' = {-# SCC "eigenvaluesSH'" #-} eigOnlySH
272
273-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
274--
275-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
276eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
277eigSH m | exactHermitian m = eigSH' m
278 | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
279
280-- | Eigenvalues of a complex hermitian or real symmetric matrix.
281eigenvaluesSH :: Field t => Matrix t -> Vector Double
282eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m
283 | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
284
285--------------------------------------------------------------
286
287-- | QR factorization.
288--
289-- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
290qr :: Field t => Matrix t -> (Matrix t, Matrix t)
291qr = {-# SCC "qr" #-} unpackQR . qr'
292
293qrRaw m = qr' m
294
295{- | generate a matrix with k orthogonal columns from the output of qrRaw
296-}
297qrgr n (a,t)
298 | dim t > min (cols a) (rows a) || n < 0 || n > dim t = error "qrgr expects k <= min(rows,cols)"
299 | otherwise = qrgr' n (a,t)
300
301-- | RQ factorization.
302--
303-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
304rq :: Field t => Matrix t -> (Matrix t, Matrix t)
305rq m = {-# SCC "rq" #-} (r,q) where
306 (q',r') = qr $ trans $ rev1 m
307 r = rev2 (trans r')
308 q = rev2 (trans q')
309 rev1 = flipud . fliprl
310 rev2 = fliprl . flipud
311
312-- | Hessenberg factorization.
313--
314-- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
315-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
316hess :: Field t => Matrix t -> (Matrix t, Matrix t)
317hess = hess'
318
319-- | Schur factorization.
320--
321-- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
322-- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
323-- upper triangular in 2x2 blocks.
324--
325-- \"Anything that the Jordan decomposition can do, the Schur decomposition
326-- can do better!\" (Van Loan)
327schur :: Field t => Matrix t -> (Matrix t, Matrix t)
328schur = schur'
329
330
331-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
332mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
333mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH'
334
335-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
336cholSH :: Field t => Matrix t -> Matrix t
337cholSH = {-# SCC "cholSH" #-} cholSH'
338
339-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
340--
341-- If @c = chol m@ then @c@ is upper triangular and @m == ctrans c \<> c@.
342chol :: Field t => Matrix t -> Matrix t
343chol m | exactHermitian m = cholSH m
344 | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
345
346
347-- | Joint computation of inverse and logarithm of determinant of a square matrix.
348invlndet :: Field t
349 => Matrix t
350 -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
351invlndet m | square m = (im,(ladm,sdm))
352 | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
353 where
354 lp@(lup,perm) = luPacked m
355 s = signlp (rows m) perm
356 dg = toList $ takeDiag $ lup
357 ladm = sum $ map (log.abs) dg
358 sdm = s* product (map signum dg)
359 im = luSolve lp (ident (rows m))
360
361
362-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
363det :: Field t => Matrix t -> t
364det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
365 | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"
366 where (lup,perm) = luPacked m
367 s = signlp (rows m) perm
368
369-- | Explicit LU factorization of a general matrix.
370--
371-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
372-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
373lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
374lu = luFact . luPacked
375
376-- | Inverse of a square matrix. See also 'invlndet'.
377inv :: Field t => Matrix t -> Matrix t
378inv m | square m = m `linearSolve` ident (rows m)
379 | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"
380
381
382-- | Pseudoinverse of a general matrix with default tolerance ('pinvTol' 1, similar to GNU-Octave).
383pinv :: Field t => Matrix t -> Matrix t
384pinv = pinvTol 1
385
386{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.
387
388@
389m = (3><3) [ 1, 0, 0
390 , 0, 1, 0
391 , 0, 0, 1e-10] :: Matrix Double
392@
393
394>>> pinv m
3951. 0. 0.
3960. 1. 0.
3970. 0. 10000000000.
398
399>>> pinvTol 1E8 m
4001. 0. 0.
4010. 1. 0.
4020. 0. 1.
403
404-}
405
406pinvTol :: Field t => Double -> Matrix t -> Matrix t
407pinvTol t m = conj v' `mXm` diag s' `mXm` ctrans u' where
408 (u,s,v) = thinSVD m
409 sl@(g:_) = toList s
410 s' = real . fromList . map rec $ sl
411 rec x = if x <= g*tol then x else 1/x
412 tol = (fromIntegral (max r c) * g * t * eps)
413 r = rows m
414 c = cols m
415 d = dim s
416 u' = takeColumns d u
417 v' = takeColumns d v
418
419
420-- | Numeric rank of a matrix from the SVD decomposition.
421rankSVD :: Element t
422 => Double -- ^ numeric zero (e.g. 1*'eps')
423 -> Matrix t -- ^ input matrix m
424 -> Vector Double -- ^ 'sv' of m
425 -> Int -- ^ rank of m
426rankSVD teps m s = ranksv teps (max (rows m) (cols m)) (toList s)
427
428-- | Numeric rank of a matrix from its singular values.
429ranksv :: Double -- ^ numeric zero (e.g. 1*'eps')
430 -> Int -- ^ maximum dimension of the matrix
431 -> [Double] -- ^ singular values
432 -> Int -- ^ rank of m
433ranksv teps maxdim s = k where
434 g = maximum s
435 tol = fromIntegral maxdim * g * teps
436 s' = filter (>tol) s
437 k = if g > teps then length s' else 0
438
439-- | The machine precision of a Double: @eps = 2.22044604925031e-16@ (the value used by GNU-Octave).
440eps :: Double
441eps = 2.22044604925031e-16
442
443
444-- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1
445peps :: RealFloat x => x
446peps = x where x = 2.0 ** fromIntegral (1 - floatDigits x)
447
448
449-- | The imaginary unit: @i = 0.0 :+ 1.0@
450i :: Complex Double
451i = 0:+1
452
453-----------------------------------------------------------------------
454
455-- | The nullspace of a matrix from its precomputed SVD decomposition.
456nullspaceSVD :: Field t
457 => Either Double Int -- ^ Left \"numeric\" zero (eg. 1*'eps'),
458 -- or Right \"theoretical\" matrix rank.
459 -> Matrix t -- ^ input matrix m
460 -> (Vector Double, Matrix t) -- ^ 'rightSV' of m
461 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
462nullspaceSVD hint a (s,v) = vs where
463 tol = case hint of
464 Left t -> t
465 _ -> eps
466 k = case hint of
467 Right t -> t
468 _ -> rankSVD tol a s
469 vs = drop k $ toRows $ ctrans v
470
471
472-- | The nullspace of a matrix. See also 'nullspaceSVD'.
473nullspacePrec :: Field t
474 => Double -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
475 -> Matrix t -- ^ input matrix
476 -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
477nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
478
479-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
480nullVector :: Field t => Matrix t -> Vector t
481nullVector = last . nullspacePrec 1
482
483orth :: Field t => Matrix t -> [Vector t]
484-- ^ Return an orthonormal basis of the range space of a matrix
485orth m = take r $ toColumns u
486 where
487 (u,s,_) = compactSVD m
488 r = ranksv eps (max (rows m) (cols m)) (toList s)
489
490------------------------------------------------------------------------
491
492-- many thanks, quickcheck!
493
494haussholder :: (Field a) => a -> Vector a -> Matrix a
495haussholder tau v = ident (dim v) `sub` (tau `scale` (w `mXm` ctrans w))
496 where w = asColumn v
497
498
499zh k v = fromList $ replicate (k-1) 0 ++ (1:drop k xs)
500 where xs = toList v
501
502zt 0 v = v
503zt k v = vjoin [subVector 0 (dim v - k) v, konst' 0 k]
504
505
506unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t)
507unpackQR (pq, tau) = {-# SCC "unpackQR" #-} (q,r)
508 where cs = toColumns pq
509 m = rows pq
510 n = cols pq
511 mn = min m n
512 r = fromColumns $ zipWith zt ([m-1, m-2 .. 1] ++ repeat 0) cs
513 vs = zipWith zh [1..mn] cs
514 hs = zipWith haussholder (toList tau) vs
515 q = foldl1' mXm hs
516
517unpackHess :: (Field t) => (Matrix t -> (Matrix t,Vector t)) -> Matrix t -> (Matrix t, Matrix t)
518unpackHess hf m
519 | rows m == 1 = ((1><1)[1],m)
520 | otherwise = (uH . hf) m
521
522uH (pq, tau) = (p,h)
523 where cs = toColumns pq
524 m = rows pq
525 n = cols pq
526 mn = min m n
527 h = fromColumns $ zipWith zt ([m-2, m-3 .. 1] ++ repeat 0) cs
528 vs = zipWith zh [2..mn] cs
529 hs = zipWith haussholder (toList tau) vs
530 p = foldl1' mXm hs
531
532--------------------------------------------------------------------------
533
534-- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
535rcond :: Field t => Matrix t -> Double
536rcond m = last s / head s
537 where s = toList (singularValues m)
538
539-- | Number of linearly independent rows or columns.
540rank :: Field t => Matrix t -> Int
541rank m = rankSVD eps m (singularValues m)
542
543{-
544expm' m = case diagonalize (complex m) of
545 Just (l,v) -> v `mXm` diag (exp l) `mXm` inv v
546 Nothing -> error "Sorry, expm not yet implemented for non-diagonalizable matrices"
547 where exp = vectorMapC Exp
548-}
549
550diagonalize m = if rank v == n
551 then Just (l,v)
552 else Nothing
553 where n = rows m
554 (l,v) = if exactHermitian m
555 then let (l',v') = eigSH m in (real l', v')
556 else eig m
557
558-- | Generic matrix functions for diagonalizable matrices. For instance:
559--
560-- @logm = matFunc log@
561--
562matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
563matFunc f m = case diagonalize m of
564 Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v
565 Nothing -> error "Sorry, matFunc requires a diagonalizable matrix"
566
567--------------------------------------------------------------
568
569golubeps :: Integer -> Integer -> Double
570golubeps p q = a * fromIntegral b / fromIntegral c where
571 a = 2^^(3-p-q)
572 b = fact p * fact q
573 c = fact (p+q) * fact (p+q+1)
574 fact n = product [1..n]
575
576epslist :: [(Int,Double)]
577epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
578
579geps delta = head [ k | (k,g) <- epslist, g<delta]
580
581
582{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
583 based on a scaled Pade approximation.
584-}
585expm :: Field t => Matrix t -> Matrix t
586expm = expGolub
587
588expGolub :: Field t => Matrix t -> Matrix t
589expGolub m = iterate msq f !! j
590 where j = max 0 $ floor $ logBase 2 $ pnorm Infinity m
591 a = m */ fromIntegral ((2::Int)^j)
592 q = geps eps -- 7 steps
593 eye = ident (rows m)
594 work (k,c,x,n,d) = (k',c',x',n',d')
595 where k' = k+1
596 c' = c * fromIntegral (q-k+1) / fromIntegral ((2*q-k+1)*k)
597 x' = a <> x
598 n' = n |+| (c' .* x')
599 d' = d |+| (((-1)^k * c') .* x')
600 (_,_,_,nf,df) = iterate work (1,1,eye,eye,eye) !! q
601 f = linearSolve df nf
602 msq x = x <> x
603
604 (<>) = multiply
605 v */ x = scale (recip x) v
606 (.*) = scale
607 (|+|) = add
608
609--------------------------------------------------------------
610
611{- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia.
612It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@.
613
614@m = (2><2) [4,9
615 ,0,4] :: Matrix Double@
616
617>>> sqrtm m
618(2><2)
619 [ 2.0, 2.25
620 , 0.0, 2.0 ]
621
622-}
623sqrtm :: Field t => Matrix t -> Matrix t
624sqrtm = sqrtmInv
625
626sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
627 where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a
628 | otherwise = fixedPoint (b:rest)
629 fixedPoint _ = error "fixedpoint with impossible inputs"
630 f (y,z) = (0.5 .* (y |+| inv z),
631 0.5 .* (inv y |+| z))
632 (.*) = scale
633 (|+|) = add
634 (|-|) = sub
635
636------------------------------------------------------------------
637
638signlp r vals = foldl f 1 (zip [0..r-1] vals)
639 where f s (a,b) | a /= b = -s
640 | otherwise = s
641
642swap (arr,s) (a,b) | a /= b = (arr // [(a, arr!b),(b,arr!a)],-s)
643 | otherwise = (arr,s)
644
645fixPerm r vals = (fromColumns $ elems res, sign)
646 where v = [0..r-1]
647 s = toColumns (ident r)
648 (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
649
650triang r c h v = (r><c) [el s t | s<-[0..r-1], t<-[0..c-1]]
651 where el p q = if q-p>=h then v else 1 - v
652
653luFact (l_u,perm) | r <= c = (l ,u ,p, s)
654 | otherwise = (l',u',p, s)
655 where
656 r = rows l_u
657 c = cols l_u
658 tu = triang r c 0 1
659 tl = triang r c 0 0
660 l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst' 1 r) r r
661 u = l_u |*| tu
662 (p,s) = fixPerm r perm
663 l' = (l_u |*| tl) |+| diagRect 0 (konst' 1 c) r c
664 u' = takeRows c (l_u |*| tu)
665 (|+|) = add
666 (|*|) = mul
667
668---------------------------------------------------------------------------
669
670data NormType = Infinity | PNorm1 | PNorm2 | Frobenius
671
672class (RealFloat (RealOf t)) => Normed c t where
673 pnorm :: NormType -> c t -> RealOf t
674
675instance Normed Vector Double where
676 pnorm PNorm1 = norm1
677 pnorm PNorm2 = norm2
678 pnorm Infinity = normInf
679 pnorm Frobenius = norm2
680
681instance Normed Vector (Complex Double) where
682 pnorm PNorm1 = norm1
683 pnorm PNorm2 = norm2
684 pnorm Infinity = normInf
685 pnorm Frobenius = pnorm PNorm2
686
687instance Normed Vector Float where
688 pnorm PNorm1 = norm1
689 pnorm PNorm2 = norm2
690 pnorm Infinity = normInf
691 pnorm Frobenius = pnorm PNorm2
692
693instance Normed Vector (Complex Float) where
694 pnorm PNorm1 = norm1
695 pnorm PNorm2 = norm2
696 pnorm Infinity = normInf
697 pnorm Frobenius = pnorm PNorm2
698
699
700instance Normed Matrix Double where
701 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
702 pnorm PNorm2 = (@>0) . singularValues
703 pnorm Infinity = pnorm PNorm1 . trans
704 pnorm Frobenius = pnorm PNorm2 . flatten
705
706instance Normed Matrix (Complex Double) where
707 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
708 pnorm PNorm2 = (@>0) . singularValues
709 pnorm Infinity = pnorm PNorm1 . trans
710 pnorm Frobenius = pnorm PNorm2 . flatten
711
712instance Normed Matrix Float where
713 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
714 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
715 pnorm Infinity = pnorm PNorm1 . trans
716 pnorm Frobenius = pnorm PNorm2 . flatten
717
718instance Normed Matrix (Complex Float) where
719 pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns
720 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
721 pnorm Infinity = pnorm PNorm1 . trans
722 pnorm Frobenius = pnorm PNorm2 . flatten
723
724-- | Approximate number of common digits in the maximum element.
725relativeError :: (Normed c t, Container c t) => c t -> c t -> Int
726relativeError x y = dig (norm (x `sub` y) / norm x)
727 where norm = pnorm Infinity
728 dig r = round $ -logBase 10 (realToFrac r :: Double)
729
730----------------------------------------------------------------------
731
732-- | Generalized symmetric positive definite eigensystem Av = lBv,
733-- for A and B symmetric, B positive definite (conditions not checked).
734geigSH' :: Field t
735 => Matrix t -- ^ A
736 -> Matrix t -- ^ B
737 -> (Vector Double, Matrix t)
738geigSH' a b = (l,v')
739 where
740 u = cholSH b
741 iu = inv u
742 c = ctrans iu <> a <> iu
743 (l,v) = eigSH' c
744 v' = iu <> v
745 (<>) = mXm
746
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs
new file mode 100644
index 0000000..11394a6
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs
@@ -0,0 +1,555 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.LinearAlgebra.LAPACK
4-- Copyright : (c) Alberto Ruiz 2006-7
5-- License : GPL-style
6--
7-- Maintainer : Alberto Ruiz (aruiz at um dot es)
8-- Stability : provisional
9-- Portability : portable (uses FFI)
10--
11-- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).
12--
13-----------------------------------------------------------------------------
14{-# OPTIONS_HADDOCK hide #-}
15
16module Numeric.LinearAlgebra.LAPACK (
17 -- * Matrix product
18 multiplyR, multiplyC, multiplyF, multiplyQ,
19 -- * Linear systems
20 linearSolveR, linearSolveC,
21 lusR, lusC,
22 cholSolveR, cholSolveC,
23 linearSolveLSR, linearSolveLSC,
24 linearSolveSVDR, linearSolveSVDC,
25 -- * SVD
26 svR, svRd, svC, svCd,
27 svdR, svdRd, svdC, svdCd,
28 thinSVDR, thinSVDRd, thinSVDC, thinSVDCd,
29 rightSVR, rightSVC, leftSVR, leftSVC,
30 -- * Eigensystems
31 eigR, eigC, eigS, eigS', eigH, eigH',
32 eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH,
33 -- * LU
34 luR, luC,
35 -- * Cholesky
36 cholS, cholH, mbCholS, mbCholH,
37 -- * QR
38 qrR, qrC, qrgrR, qrgrC,
39 -- * Hessenberg
40 hessR, hessC,
41 -- * Schur
42 schurR, schurC
43) where
44
45import Data.Packed.Internal
46import Data.Packed.Matrix
47import Numeric.Conversion
48import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale))
49
50import Foreign.Ptr(nullPtr)
51import Foreign.C.Types
52import Control.Monad(when)
53import System.IO.Unsafe(unsafePerformIO)
54
55-----------------------------------------------------------------------------------
56
57foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM
58foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM
59foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM
60foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM
61
62isT Matrix{order = ColumnMajor} = 0
63isT Matrix{order = RowMajor} = 1
64
65tt x@Matrix{order = ColumnMajor} = x
66tt x@Matrix{order = RowMajor} = trans x
67
68multiplyAux f st a b = unsafePerformIO $ do
69 when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++
70 show (rows a,cols a) ++ " x " ++ show (rows b, cols b)
71 s <- createMatrix ColumnMajor (rows a) (cols b)
72 app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st
73 return s
74
75-- | Matrix product based on BLAS's /dgemm/.
76multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
77multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b
78
79-- | Matrix product based on BLAS's /zgemm/.
80multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
81multiplyC a b = multiplyAux zgemmc "zgemmc" a b
82
83-- | Matrix product based on BLAS's /sgemm/.
84multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
85multiplyF a b = multiplyAux sgemmc "sgemmc" a b
86
87-- | Matrix product based on BLAS's /cgemm/.
88multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)
89multiplyQ a b = multiplyAux cgemmc "cgemmc" a b
90
91-----------------------------------------------------------------------------
92foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM
93foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM
94foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM
95foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM
96
97-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
98svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
99svdR = svdAux dgesvd "svdR" . fmat
100
101-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
102svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
103svdRd = svdAux dgesdd "svdRdd" . fmat
104
105-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
106svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
107svdC = svdAux zgesvd "svdC" . fmat
108
109-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
110svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
111svdCd = svdAux zgesdd "svdCdd" . fmat
112
113svdAux f st x = unsafePerformIO $ do
114 u <- createMatrix ColumnMajor r r
115 s <- createVector (min r c)
116 v <- createMatrix ColumnMajor c c
117 app4 f mat x mat u vec s mat v st
118 return (u,s,trans v)
119 where r = rows x
120 c = cols x
121
122
123-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
124thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
125thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat
126
127-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
128thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
129thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat
130
131-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
132thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
133thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat
134
135-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
136thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
137thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat
138
139thinSVDAux f st x = unsafePerformIO $ do
140 u <- createMatrix ColumnMajor r q
141 s <- createVector q
142 v <- createMatrix ColumnMajor q c
143 app4 f mat x mat u vec s mat v st
144 return (u,s,trans v)
145 where r = rows x
146 c = cols x
147 q = min r c
148
149
150-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
151svR :: Matrix Double -> Vector Double
152svR = svAux dgesvd "svR" . fmat
153
154-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
155svC :: Matrix (Complex Double) -> Vector Double
156svC = svAux zgesvd "svC" . fmat
157
158-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
159svRd :: Matrix Double -> Vector Double
160svRd = svAux dgesdd "svRd" . fmat
161
162-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
163svCd :: Matrix (Complex Double) -> Vector Double
164svCd = svAux zgesdd "svCd" . fmat
165
166svAux f st x = unsafePerformIO $ do
167 s <- createVector q
168 app2 g mat x vec s st
169 return s
170 where r = rows x
171 c = cols x
172 q = min r c
173 g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr
174
175
176-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
177rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
178rightSVR = rightSVAux dgesvd "rightSVR" . fmat
179
180-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
181rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
182rightSVC = rightSVAux zgesvd "rightSVC" . fmat
183
184rightSVAux f st x = unsafePerformIO $ do
185 s <- createVector q
186 v <- createMatrix ColumnMajor c c
187 app3 g mat x vec s mat v st
188 return (s,trans v)
189 where r = rows x
190 c = cols x
191 q = min r c
192 g ra ca pa = f ra ca pa 0 0 nullPtr
193
194
195-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'.
196leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
197leftSVR = leftSVAux dgesvd "leftSVR" . fmat
198
199-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
200leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
201leftSVC = leftSVAux zgesvd "leftSVC" . fmat
202
203leftSVAux f st x = unsafePerformIO $ do
204 u <- createMatrix ColumnMajor r r
205 s <- createVector q
206 app3 g mat x mat u vec s st
207 return (u,s)
208 where r = rows x
209 c = cols x
210 q = min r c
211 g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr
212
213-----------------------------------------------------------------------------
214
215foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM
216foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM
217foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM
218foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM
219
220eigAux f st m = unsafePerformIO $ do
221 l <- createVector r
222 v <- createMatrix ColumnMajor r r
223 app3 g mat m vec l mat v st
224 return (l,v)
225 where r = rows m
226 g ra ca pa = f ra ca pa 0 0 nullPtr
227
228
229-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
230-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
231eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
232eigC = eigAux zgeev "eigC" . fmat
233
234eigOnlyAux f st m = unsafePerformIO $ do
235 l <- createVector r
236 app2 g mat m vec l st
237 return l
238 where r = rows m
239 g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr
240
241-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
242-- The eigenvalues are not sorted.
243eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
244eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat
245
246-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.
247-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
248eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
249eigR m = (s', v'')
250 where (s,v) = eigRaux (fmat m)
251 s' = fixeig1 s
252 v' = toRows $ trans v
253 v'' = fromColumns $ fixeig (toList s') v'
254
255eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
256eigRaux m = unsafePerformIO $ do
257 l <- createVector r
258 v <- createMatrix ColumnMajor r r
259 app3 g mat m vec l mat v "eigR"
260 return (l,v)
261 where r = rows m
262 g ra ca pa = dgeev ra ca pa 0 0 nullPtr
263
264fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s))
265 where r = dim s
266
267fixeig [] _ = []
268fixeig [_] [v] = [comp' v]
269fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
270 | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs
271 | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs)
272 where scale = vectorMapValR Scale
273fixeig _ _ = error "fixeig with impossible inputs"
274
275
276-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
277-- The eigenvalues are not sorted.
278eigOnlyR :: Matrix Double -> Vector (Complex Double)
279eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
280
281
282-----------------------------------------------------------------------------
283
284eigSHAux f st m = unsafePerformIO $ do
285 l <- createVector r
286 v <- createMatrix ColumnMajor r r
287 app3 f mat m vec l mat v st
288 return (l,v)
289 where r = rows m
290
291-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.
292-- The eigenvectors are the columns of v.
293-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order).
294eigS :: Matrix Double -> (Vector Double, Matrix Double)
295eigS m = (s', fliprl v)
296 where (s,v) = eigS' (fmat m)
297 s' = fromList . reverse . toList $ s
298
299-- | 'eigS' in ascending order
300eigS' :: Matrix Double -> (Vector Double, Matrix Double)
301eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
302
303-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.
304-- The eigenvectors are the columns of v.
305-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order).
306eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
307eigH m = (s', fliprl v)
308 where (s,v) = eigH' (fmat m)
309 s' = fromList . reverse . toList $ s
310
311-- | 'eigH' in ascending order
312eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
313eigH' = eigSHAux (zheev 1) "eigH'" . fmat
314
315
316-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.
317-- The eigenvalues are sorted in descending order.
318eigOnlyS :: Matrix Double -> Vector Double
319eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat
320
321-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.
322-- The eigenvalues are sorted in descending order.
323eigOnlyH :: Matrix (Complex Double) -> Vector Double
324eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat
325
326vrev = flatten . flipud . reshape 1
327
328-----------------------------------------------------------------------------
329foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM
330foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM
331foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM
332foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM
333
334linearSolveSQAux f st a b
335 | n1==n2 && n1==r = unsafePerformIO $ do
336 s <- createMatrix ColumnMajor r c
337 app3 f mat a mat b mat s st
338 return s
339 | otherwise = error $ st ++ " of nonsquare matrix"
340 where n1 = rows a
341 n2 = cols a
342 r = rows b
343 c = cols b
344
345-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.
346linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
347linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b)
348
349-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'.
350linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
351linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b)
352
353
354-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.
355cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
356cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b)
357
358-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'.
359cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
360cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b)
361
362-----------------------------------------------------------------------------------
363foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM
364foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM
365foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM
366foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM
367
368linearSolveAux f st a b = unsafePerformIO $ do
369 r <- createMatrix ColumnMajor (max m n) nrhs
370 app3 f mat a mat b mat r st
371 return r
372 where m = rows a
373 n = cols a
374 nrhs = cols b
375
376-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'.
377linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
378linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $
379 linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b)
380
381-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'.
382linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
383linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $
384 linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b)
385
386-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. 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.
387linearSolveSVDR :: Maybe Double -- ^ rcond
388 -> Matrix Double -- ^ coefficient matrix
389 -> Matrix Double -- ^ right hand sides (as columns)
390 -> Matrix Double -- ^ solution vectors (as columns)
391linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
392 linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b)
393linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b)
394
395-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. 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.
396linearSolveSVDC :: Maybe Double -- ^ rcond
397 -> Matrix (Complex Double) -- ^ coefficient matrix
398 -> Matrix (Complex Double) -- ^ right hand sides (as columns)
399 -> Matrix (Complex Double) -- ^ solution vectors (as columns)
400linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
401 linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b)
402linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
403
404-----------------------------------------------------------------------------------
405foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM
406foreign import ccall unsafe "chol_l_S" dpotrf :: TMM
407
408cholAux f st a = do
409 r <- createMatrix ColumnMajor n n
410 app2 f mat a mat r st
411 return r
412 where n = rows a
413
414-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
415cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
416cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat
417
418-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
419cholS :: Matrix Double -> Matrix Double
420cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat
421
422-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version).
423mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
424mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat
425
426-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version).
427mbCholS :: Matrix Double -> Maybe (Matrix Double)
428mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat
429
430-----------------------------------------------------------------------------------
431foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM
432foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM
433
434-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
435qrR :: Matrix Double -> (Matrix Double, Vector Double)
436qrR = qrAux dgeqr2 "qrR" . fmat
437
438-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.
439qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
440qrC = qrAux zgeqr2 "qrC" . fmat
441
442qrAux f st a = unsafePerformIO $ do
443 r <- createMatrix ColumnMajor m n
444 tau <- createVector mn
445 app3 f mat a vec tau mat r st
446 return (r,tau)
447 where
448 m = rows a
449 n = cols a
450 mn = min m n
451
452foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM
453foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM
454
455-- | build rotation from reflectors
456qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
457qrgrR = qrgrAux dorgqr "qrgrR"
458-- | build rotation from reflectors
459qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
460qrgrC = qrgrAux zungqr "qrgrC"
461
462qrgrAux f st n (a, tau) = unsafePerformIO $ do
463 res <- createMatrix ColumnMajor (rows a) n
464 app3 f mat (fmat a) vec (subVector 0 n tau') mat res st
465 return res
466 where
467 tau' = vjoin [tau, constantD 0 n]
468
469-----------------------------------------------------------------------------------
470foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM
471foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM
472
473-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
474hessR :: Matrix Double -> (Matrix Double, Vector Double)
475hessR = hessAux dgehrd "hessR" . fmat
476
477-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.
478hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
479hessC = hessAux zgehrd "hessC" . fmat
480
481hessAux f st a = unsafePerformIO $ do
482 r <- createMatrix ColumnMajor m n
483 tau <- createVector (mn-1)
484 app3 f mat a vec tau mat r st
485 return (r,tau)
486 where m = rows a
487 n = cols a
488 mn = min m n
489
490-----------------------------------------------------------------------------------
491foreign import ccall unsafe "schur_l_R" dgees :: TMMM
492foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM
493
494-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
495schurR :: Matrix Double -> (Matrix Double, Matrix Double)
496schurR = schurAux dgees "schurR" . fmat
497
498-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.
499schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
500schurC = schurAux zgees "schurC" . fmat
501
502schurAux f st a = unsafePerformIO $ do
503 u <- createMatrix ColumnMajor n n
504 s <- createMatrix ColumnMajor n n
505 app3 f mat a mat u mat s st
506 return (u,s)
507 where n = rows a
508
509-----------------------------------------------------------------------------------
510foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM
511foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM
512
513-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
514luR :: Matrix Double -> (Matrix Double, [Int])
515luR = luAux dgetrf "luR" . fmat
516
517-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
518luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
519luC = luAux zgetrf "luC" . fmat
520
521luAux f st a = unsafePerformIO $ do
522 lu <- createMatrix ColumnMajor n m
523 piv <- createVector (min n m)
524 app3 f mat a vec piv mat lu st
525 return (lu, map (pred.round) (toList piv))
526 where n = rows a
527 m = cols a
528
529-----------------------------------------------------------------------------------
530type TW a = CInt -> PD -> a
531type TQ a = CInt -> CInt -> PC -> a
532
533foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM
534foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt))))
535
536-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/.
537lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
538lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b)
539
540-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.
541lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
542lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b)
543
544lusAux f st a piv b
545 | n1==n2 && n2==n =unsafePerformIO $ do
546 x <- createMatrix ColumnMajor n m
547 app4 f mat a vec piv' mat b mat x st
548 return x
549 | otherwise = error $ st ++ " on LU factorization of nonsquare matrix"
550 where n1 = rows a
551 n2 = cols a
552 n = rows b
553 m = cols b
554 piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double
555
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
new file mode 100644
index 0000000..e5e45ef
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -0,0 +1,1489 @@
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#include <math.h>
5#include <time.h>
6#include "lapack-aux.h"
7
8#define MACRO(B) do {B} while (0)
9#define ERROR(CODE) MACRO(return CODE;)
10#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
11
12#define MIN(A,B) ((A)<(B)?(A):(B))
13#define MAX(A,B) ((A)>(B)?(A):(B))
14
15// #define DBGL
16
17#ifdef DBGL
18#define DEBUGMSG(M) printf("\nLAPACK "M"\n");
19#else
20#define DEBUGMSG(M)
21#endif
22
23#define OK return 0;
24
25// #ifdef DBGL
26// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL);
27// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
28// #else
29// #define DEBUGMSG(M)
30// #define OK return 0;
31// #endif
32
33#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \
34 for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");}
35
36#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
37
38#define BAD_SIZE 2000
39#define BAD_CODE 2001
40#define MEM 2002
41#define BAD_FILE 2003
42#define SINGULAR 2004
43#define NOCONVER 2005
44#define NODEFPOS 2006
45#define NOSPRTD 2007
46
47//---------------------------------------
48void asm_finit() {
49#ifdef i386
50
51// asm("finit");
52
53 static unsigned char buf[108];
54 asm("FSAVE %0":"=m" (buf));
55
56 #if FPUDEBUG
57 if(buf[8]!=255 || buf[9]!=255) { // print warning in red
58 printf("%c[;31mWarning: FPU TAG = %x %x\%c[0m\n",0x1B,buf[8],buf[9],0x1B);
59 }
60 #endif
61
62 #if NANDEBUG
63 asm("FRSTOR %0":"=m" (buf));
64 #endif
65
66#endif
67}
68
69//---------------------------------------
70
71#if NANDEBUG
72
73#define CHECKNANR(M,msg) \
74{ int k; \
75for(k=0; k<(M##r * M##c); k++) { \
76 if(M##p[k] != M##p[k]) { \
77 printf(msg); \
78 TRACEMAT(M) \
79 /*exit(1);*/ \
80 } \
81} \
82}
83
84#define CHECKNANC(M,msg) \
85{ int k; \
86for(k=0; k<(M##r * M##c); k++) { \
87 if( M##p[k].r != M##p[k].r \
88 || M##p[k].i != M##p[k].i) { \
89 printf(msg); \
90 /*exit(1);*/ \
91 } \
92} \
93}
94
95#else
96#define CHECKNANC(M,msg)
97#define CHECKNANR(M,msg)
98#endif
99
100//---------------------------------------
101
102//////////////////// real svd ////////////////////////////////////
103
104/* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
105 doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
106 ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
107 integer *info);
108
109int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
110 integer m = ar;
111 integer n = ac;
112 integer q = MIN(m,n);
113 REQUIRES(sn==q,BAD_SIZE);
114 REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
115 char* jobu = "A";
116 if (up==NULL) {
117 jobu = "N";
118 } else {
119 if (uc==q) {
120 jobu = "S";
121 }
122 }
123 REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
124 char* jobvt = "A";
125 integer ldvt = n;
126 if (vp==NULL) {
127 jobvt = "N";
128 } else {
129 if (vr==q) {
130 jobvt = "S";
131 ldvt = q;
132 }
133 }
134 DEBUGMSG("svd_l_R");
135 double *B = (double*)malloc(m*n*sizeof(double));
136 CHECK(!B,MEM);
137 memcpy(B,ap,m*n*sizeof(double));
138 integer lwork = -1;
139 integer res;
140 // ask for optimal lwork
141 double ans;
142 dgesvd_ (jobu,jobvt,
143 &m,&n,B,&m,
144 sp,
145 up,&m,
146 vp,&ldvt,
147 &ans, &lwork,
148 &res);
149 lwork = ceil(ans);
150 double * work = (double*)malloc(lwork*sizeof(double));
151 CHECK(!work,MEM);
152 dgesvd_ (jobu,jobvt,
153 &m,&n,B,&m,
154 sp,
155 up,&m,
156 vp,&ldvt,
157 work, &lwork,
158 &res);
159 CHECK(res,res);
160 free(work);
161 free(B);
162 OK
163}
164
165// (alternative version)
166
167/* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
168 a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
169 doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
170 integer *iwork, integer *info);
171
172int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) {
173 integer m = ar;
174 integer n = ac;
175 integer q = MIN(m,n);
176 REQUIRES(sn==q,BAD_SIZE);
177 REQUIRES((up == NULL && vp == NULL)
178 || (ur==m && vc==n
179 && ((uc == q && vr == q)
180 || (uc == m && vc==n))),BAD_SIZE);
181 char* jobz = "A";
182 integer ldvt = n;
183 if (up==NULL) {
184 jobz = "N";
185 } else {
186 if (uc==q && vr == q) {
187 jobz = "S";
188 ldvt = q;
189 }
190 }
191 DEBUGMSG("svd_l_Rdd");
192 double *B = (double*)malloc(m*n*sizeof(double));
193 CHECK(!B,MEM);
194 memcpy(B,ap,m*n*sizeof(double));
195 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
196 CHECK(!iwk,MEM);
197 integer lwk = -1;
198 integer res;
199 // ask for optimal lwk
200 double ans;
201 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
202 lwk = ans;
203 double * workv = (double*)malloc(lwk*sizeof(double));
204 CHECK(!workv,MEM);
205 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
206 CHECK(res,res);
207 free(iwk);
208 free(workv);
209 free(B);
210 OK
211}
212
213//////////////////// complex svd ////////////////////////////////////
214
215// not in clapack.h
216
217int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
218 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
219 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
220 integer *lwork, doublereal *rwork, integer *info);
221
222int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
223 integer m = ar;
224 integer n = ac;
225 integer q = MIN(m,n);
226 REQUIRES(sn==q,BAD_SIZE);
227 REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE);
228 char* jobu = "A";
229 if (up==NULL) {
230 jobu = "N";
231 } else {
232 if (uc==q) {
233 jobu = "S";
234 }
235 }
236 REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE);
237 char* jobvt = "A";
238 integer ldvt = n;
239 if (vp==NULL) {
240 jobvt = "N";
241 } else {
242 if (vr==q) {
243 jobvt = "S";
244 ldvt = q;
245 }
246 }DEBUGMSG("svd_l_C");
247 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
248 CHECK(!B,MEM);
249 memcpy(B,ap,m*n*sizeof(doublecomplex));
250
251 double *rwork = (double*) malloc(5*q*sizeof(double));
252 CHECK(!rwork,MEM);
253 integer lwork = -1;
254 integer res;
255 // ask for optimal lwork
256 doublecomplex ans;
257 zgesvd_ (jobu,jobvt,
258 &m,&n,B,&m,
259 sp,
260 up,&m,
261 vp,&ldvt,
262 &ans, &lwork,
263 rwork,
264 &res);
265 lwork = ceil(ans.r);
266 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
267 CHECK(!work,MEM);
268 zgesvd_ (jobu,jobvt,
269 &m,&n,B,&m,
270 sp,
271 up,&m,
272 vp,&ldvt,
273 work, &lwork,
274 rwork,
275 &res);
276 CHECK(res,res);
277 free(work);
278 free(rwork);
279 free(B);
280 OK
281}
282
283int zgesdd_ (char *jobz, integer *m, integer *n,
284 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
285 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
286 integer *lwork, doublereal *rwork, integer* iwork, integer *info);
287
288int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
289 //printf("entro\n");
290 integer m = ar;
291 integer n = ac;
292 integer q = MIN(m,n);
293 REQUIRES(sn==q,BAD_SIZE);
294 REQUIRES((up == NULL && vp == NULL)
295 || (ur==m && vc==n
296 && ((uc == q && vr == q)
297 || (uc == m && vc==n))),BAD_SIZE);
298 char* jobz = "A";
299 integer ldvt = n;
300 if (up==NULL) {
301 jobz = "N";
302 } else {
303 if (uc==q && vr == q) {
304 jobz = "S";
305 ldvt = q;
306 }
307 }
308 DEBUGMSG("svd_l_Cdd");
309 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
310 CHECK(!B,MEM);
311 memcpy(B,ap,m*n*sizeof(doublecomplex));
312 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
313 CHECK(!iwk,MEM);
314 int lrwk;
315 if (0 && *jobz == 'N') {
316 lrwk = 5*q; // does not work, crash at free below
317 } else {
318 lrwk = 5*q*q + 7*q;
319 }
320 double *rwk = (double*)malloc(lrwk*sizeof(double));;
321 CHECK(!rwk,MEM);
322 //printf("%s %ld %d\n",jobz,q,lrwk);
323 integer lwk = -1;
324 integer res;
325 // ask for optimal lwk
326 doublecomplex ans;
327 zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
328 lwk = ans.r;
329 //printf("lwk = %ld\n",lwk);
330 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
331 CHECK(!workv,MEM);
332 zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
333 //printf("res = %ld\n",res);
334 CHECK(res,res);
335 free(workv); // printf("freed workv\n");
336 free(rwk); // printf("freed rwk\n");
337 free(iwk); // printf("freed iwk\n");
338 free(B); // printf("freed B, salgo\n");
339 OK
340}
341
342//////////////////// general complex eigensystem ////////////
343
344/* Subroutine */ int zgeev_(char *jobvl, char *jobvr, integer *n,
345 doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl,
346 integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work,
347 integer *lwork, doublereal *rwork, integer *info);
348
349int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) {
350 integer n = ar;
351 REQUIRES(ac==n && sn==n, BAD_SIZE);
352 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
353 char jobvl = up==NULL?'N':'V';
354 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
355 char jobvr = vp==NULL?'N':'V';
356 DEBUGMSG("eig_l_C");
357 doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
358 CHECK(!B,MEM);
359 memcpy(B,ap,n*n*sizeof(doublecomplex));
360 double *rwork = (double*) malloc(2*n*sizeof(double));
361 CHECK(!rwork,MEM);
362 integer lwork = -1;
363 integer res;
364 // ask for optimal lwork
365 doublecomplex ans;
366 //printf("ask zgeev\n");
367 zgeev_ (&jobvl,&jobvr,
368 &n,B,&n,
369 sp,
370 up,&n,
371 vp,&n,
372 &ans, &lwork,
373 rwork,
374 &res);
375 lwork = ceil(ans.r);
376 //printf("ans = %d\n",lwork);
377 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
378 CHECK(!work,MEM);
379 //printf("zgeev\n");
380 zgeev_ (&jobvl,&jobvr,
381 &n,B,&n,
382 sp,
383 up,&n,
384 vp,&n,
385 work, &lwork,
386 rwork,
387 &res);
388 CHECK(res,res);
389 free(work);
390 free(rwork);
391 free(B);
392 OK
393}
394
395
396
397//////////////////// general real eigensystem ////////////
398
399/* Subroutine */ int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
400 a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
401 integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
402 integer *lwork, integer *info);
403
404int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
405 integer n = ar;
406 REQUIRES(ac==n && sn==n, BAD_SIZE);
407 REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
408 char jobvl = up==NULL?'N':'V';
409 REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
410 char jobvr = vp==NULL?'N':'V';
411 DEBUGMSG("eig_l_R");
412 double *B = (double*)malloc(n*n*sizeof(double));
413 CHECK(!B,MEM);
414 memcpy(B,ap,n*n*sizeof(double));
415 integer lwork = -1;
416 integer res;
417 // ask for optimal lwork
418 double ans;
419 //printf("ask dgeev\n");
420 dgeev_ (&jobvl,&jobvr,
421 &n,B,&n,
422 (double*)sp, (double*)sp+n,
423 up,&n,
424 vp,&n,
425 &ans, &lwork,
426 &res);
427 lwork = ceil(ans);
428 //printf("ans = %d\n",lwork);
429 double * work = (double*)malloc(lwork*sizeof(double));
430 CHECK(!work,MEM);
431 //printf("dgeev\n");
432 dgeev_ (&jobvl,&jobvr,
433 &n,B,&n,
434 (double*)sp, (double*)sp+n,
435 up,&n,
436 vp,&n,
437 work, &lwork,
438 &res);
439 CHECK(res,res);
440 free(work);
441 free(B);
442 OK
443}
444
445
446//////////////////// symmetric real eigensystem ////////////
447
448/* Subroutine */ int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
449 integer *lda, doublereal *w, doublereal *work, integer *lwork,
450 integer *info);
451
452int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
453 integer n = ar;
454 REQUIRES(ac==n && sn==n, BAD_SIZE);
455 REQUIRES(vr==n && vc==n, BAD_SIZE);
456 char jobz = wantV?'V':'N';
457 DEBUGMSG("eig_l_S");
458 memcpy(vp,ap,n*n*sizeof(double));
459 integer lwork = -1;
460 char uplo = 'U';
461 integer res;
462 // ask for optimal lwork
463 double ans;
464 //printf("ask dsyev\n");
465 dsyev_ (&jobz,&uplo,
466 &n,vp,&n,
467 sp,
468 &ans, &lwork,
469 &res);
470 lwork = ceil(ans);
471 //printf("ans = %d\n",lwork);
472 double * work = (double*)malloc(lwork*sizeof(double));
473 CHECK(!work,MEM);
474 dsyev_ (&jobz,&uplo,
475 &n,vp,&n,
476 sp,
477 work, &lwork,
478 &res);
479 CHECK(res,res);
480 free(work);
481 OK
482}
483
484//////////////////// hermitian complex eigensystem ////////////
485
486/* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
487 *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork,
488 doublereal *rwork, integer *info);
489
490int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) {
491 integer n = ar;
492 REQUIRES(ac==n && sn==n, BAD_SIZE);
493 REQUIRES(vr==n && vc==n, BAD_SIZE);
494 char jobz = wantV?'V':'N';
495 DEBUGMSG("eig_l_H");
496 memcpy(vp,ap,2*n*n*sizeof(double));
497 double *rwork = (double*) malloc((3*n-2)*sizeof(double));
498 CHECK(!rwork,MEM);
499 integer lwork = -1;
500 char uplo = 'U';
501 integer res;
502 // ask for optimal lwork
503 doublecomplex ans;
504 //printf("ask zheev\n");
505 zheev_ (&jobz,&uplo,
506 &n,vp,&n,
507 sp,
508 &ans, &lwork,
509 rwork,
510 &res);
511 lwork = ceil(ans.r);
512 //printf("ans = %d\n",lwork);
513 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
514 CHECK(!work,MEM);
515 zheev_ (&jobz,&uplo,
516 &n,vp,&n,
517 sp,
518 work, &lwork,
519 rwork,
520 &res);
521 CHECK(res,res);
522 free(work);
523 free(rwork);
524 OK
525}
526
527//////////////////// general real linear system ////////////
528
529/* Subroutine */ int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
530 *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);
531
532int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
533 integer n = ar;
534 integer nhrs = bc;
535 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
536 DEBUGMSG("linearSolveR_l");
537 double*AC = (double*)malloc(n*n*sizeof(double));
538 memcpy(AC,ap,n*n*sizeof(double));
539 memcpy(xp,bp,n*nhrs*sizeof(double));
540 integer * ipiv = (integer*)malloc(n*sizeof(integer));
541 integer res;
542 dgesv_ (&n,&nhrs,
543 AC, &n,
544 ipiv,
545 xp, &n,
546 &res);
547 if(res>0) {
548 return SINGULAR;
549 }
550 CHECK(res,res);
551 free(ipiv);
552 free(AC);
553 OK
554}
555
556//////////////////// general complex linear system ////////////
557
558/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a,
559 integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *
560 info);
561
562int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
563 integer n = ar;
564 integer nhrs = bc;
565 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
566 DEBUGMSG("linearSolveC_l");
567 doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
568 memcpy(AC,ap,n*n*sizeof(doublecomplex));
569 memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
570 integer * ipiv = (integer*)malloc(n*sizeof(integer));
571 integer res;
572 zgesv_ (&n,&nhrs,
573 AC, &n,
574 ipiv,
575 xp, &n,
576 &res);
577 if(res>0) {
578 return SINGULAR;
579 }
580 CHECK(res,res);
581 free(ipiv);
582 free(AC);
583 OK
584}
585
586//////// symmetric positive definite real linear system using Cholesky ////////////
587
588/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs,
589 doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
590 info);
591
592int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
593 integer n = ar;
594 integer nhrs = bc;
595 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
596 DEBUGMSG("cholSolveR_l");
597 memcpy(xp,bp,n*nhrs*sizeof(double));
598 integer res;
599 dpotrs_ ("U",
600 &n,&nhrs,
601 (double*)ap, &n,
602 xp, &n,
603 &res);
604 CHECK(res,res);
605 OK
606}
607
608//////// Hermitian positive definite real linear system using Cholesky ////////////
609
610/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs,
611 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
612 integer *info);
613
614int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
615 integer n = ar;
616 integer nhrs = bc;
617 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
618 DEBUGMSG("cholSolveC_l");
619 memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
620 integer res;
621 zpotrs_ ("U",
622 &n,&nhrs,
623 (doublecomplex*)ap, &n,
624 xp, &n,
625 &res);
626 CHECK(res,res);
627 OK
628}
629
630//////////////////// least squares real linear system ////////////
631
632/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer *
633 nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
634 doublereal *work, integer *lwork, integer *info);
635
636int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) {
637 integer m = ar;
638 integer n = ac;
639 integer nrhs = bc;
640 integer ldb = xr;
641 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
642 DEBUGMSG("linearSolveLSR_l");
643 double*AC = (double*)malloc(m*n*sizeof(double));
644 memcpy(AC,ap,m*n*sizeof(double));
645 if (m>=n) {
646 memcpy(xp,bp,m*nrhs*sizeof(double));
647 } else {
648 int k;
649 for(k = 0; k<nrhs; k++) {
650 memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
651 }
652 }
653 integer res;
654 integer lwork = -1;
655 double ans;
656 dgels_ ("N",&m,&n,&nrhs,
657 AC,&m,
658 xp,&ldb,
659 &ans,&lwork,
660 &res);
661 lwork = ceil(ans);
662 //printf("ans = %d\n",lwork);
663 double * work = (double*)malloc(lwork*sizeof(double));
664 dgels_ ("N",&m,&n,&nrhs,
665 AC,&m,
666 xp,&ldb,
667 work,&lwork,
668 &res);
669 if(res>0) {
670 return SINGULAR;
671 }
672 CHECK(res,res);
673 free(work);
674 free(AC);
675 OK
676}
677
678//////////////////// least squares complex linear system ////////////
679
680/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer *
681 nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
682 doublecomplex *work, integer *lwork, integer *info);
683
684int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) {
685 integer m = ar;
686 integer n = ac;
687 integer nrhs = bc;
688 integer ldb = xr;
689 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
690 DEBUGMSG("linearSolveLSC_l");
691 doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
692 memcpy(AC,ap,m*n*sizeof(doublecomplex));
693 if (m>=n) {
694 memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
695 } else {
696 int k;
697 for(k = 0; k<nrhs; k++) {
698 memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
699 }
700 }
701 integer res;
702 integer lwork = -1;
703 doublecomplex ans;
704 zgels_ ("N",&m,&n,&nrhs,
705 AC,&m,
706 xp,&ldb,
707 &ans,&lwork,
708 &res);
709 lwork = ceil(ans.r);
710 //printf("ans = %d\n",lwork);
711 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
712 zgels_ ("N",&m,&n,&nrhs,
713 AC,&m,
714 xp,&ldb,
715 work,&lwork,
716 &res);
717 if(res>0) {
718 return SINGULAR;
719 }
720 CHECK(res,res);
721 free(work);
722 free(AC);
723 OK
724}
725
726//////////////////// least squares real linear system using SVD ////////////
727
728/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs,
729 doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
730 s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
731 integer *info);
732
733int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) {
734 integer m = ar;
735 integer n = ac;
736 integer nrhs = bc;
737 integer ldb = xr;
738 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
739 DEBUGMSG("linearSolveSVDR_l");
740 double*AC = (double*)malloc(m*n*sizeof(double));
741 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
742 memcpy(AC,ap,m*n*sizeof(double));
743 if (m>=n) {
744 memcpy(xp,bp,m*nrhs*sizeof(double));
745 } else {
746 int k;
747 for(k = 0; k<nrhs; k++) {
748 memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
749 }
750 }
751 integer res;
752 integer lwork = -1;
753 integer rank;
754 double ans;
755 dgelss_ (&m,&n,&nrhs,
756 AC,&m,
757 xp,&ldb,
758 S,
759 &rcond,&rank,
760 &ans,&lwork,
761 &res);
762 lwork = ceil(ans);
763 //printf("ans = %d\n",lwork);
764 double * work = (double*)malloc(lwork*sizeof(double));
765 dgelss_ (&m,&n,&nrhs,
766 AC,&m,
767 xp,&ldb,
768 S,
769 &rcond,&rank,
770 work,&lwork,
771 &res);
772 if(res>0) {
773 return NOCONVER;
774 }
775 CHECK(res,res);
776 free(work);
777 free(S);
778 free(AC);
779 OK
780}
781
782//////////////////// least squares complex linear system using SVD ////////////
783
784// not in clapack.h
785
786int zgelss_(integer *m, integer *n, integer *nhrs,
787 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s,
788 doublereal *rcond, integer* rank,
789 doublecomplex *work, integer* lwork, doublereal* rwork,
790 integer *info);
791
792int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
793 integer m = ar;
794 integer n = ac;
795 integer nrhs = bc;
796 integer ldb = xr;
797 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE);
798 DEBUGMSG("linearSolveSVDC_l");
799 doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
800 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
801 double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));
802 memcpy(AC,ap,m*n*sizeof(doublecomplex));
803 if (m>=n) {
804 memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
805 } else {
806 int k;
807 for(k = 0; k<nrhs; k++) {
808 memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
809 }
810 }
811 integer res;
812 integer lwork = -1;
813 integer rank;
814 doublecomplex ans;
815 zgelss_ (&m,&n,&nrhs,
816 AC,&m,
817 xp,&ldb,
818 S,
819 &rcond,&rank,
820 &ans,&lwork,
821 RWORK,
822 &res);
823 lwork = ceil(ans.r);
824 //printf("ans = %d\n",lwork);
825 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
826 zgelss_ (&m,&n,&nrhs,
827 AC,&m,
828 xp,&ldb,
829 S,
830 &rcond,&rank,
831 work,&lwork,
832 RWORK,
833 &res);
834 if(res>0) {
835 return NOCONVER;
836 }
837 CHECK(res,res);
838 free(work);
839 free(RWORK);
840 free(S);
841 free(AC);
842 OK
843}
844
845//////////////////// Cholesky factorization /////////////////////////
846
847/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a,
848 integer *lda, integer *info);
849
850int chol_l_H(KCMAT(a),CMAT(l)) {
851 integer n = ar;
852 REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
853 DEBUGMSG("chol_l_H");
854 memcpy(lp,ap,n*n*sizeof(doublecomplex));
855 char uplo = 'U';
856 integer res;
857 zpotrf_ (&uplo,&n,lp,&n,&res);
858 CHECK(res>0,NODEFPOS);
859 CHECK(res,res);
860 doublecomplex zero = {0.,0.};
861 int r,c;
862 for (r=0; r<lr-1; r++) {
863 for(c=r+1; c<lc; c++) {
864 lp[r*lc+c] = zero;
865 }
866 }
867 OK
868}
869
870
871/* Subroutine */ int dpotrf_(char *uplo, integer *n, doublereal *a, integer *
872 lda, integer *info);
873
874int chol_l_S(KDMAT(a),DMAT(l)) {
875 integer n = ar;
876 REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
877 DEBUGMSG("chol_l_S");
878 memcpy(lp,ap,n*n*sizeof(double));
879 char uplo = 'U';
880 integer res;
881 dpotrf_ (&uplo,&n,lp,&n,&res);
882 CHECK(res>0,NODEFPOS);
883 CHECK(res,res);
884 int r,c;
885 for (r=0; r<lr-1; r++) {
886 for(c=r+1; c<lc; c++) {
887 lp[r*lc+c] = 0.;
888 }
889 }
890 OK
891}
892
893//////////////////// QR factorization /////////////////////////
894
895/* Subroutine */ int dgeqr2_(integer *m, integer *n, doublereal *a, integer *
896 lda, doublereal *tau, doublereal *work, integer *info);
897
898int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
899 integer m = ar;
900 integer n = ac;
901 integer mn = MIN(m,n);
902 REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
903 DEBUGMSG("qr_l_R");
904 double *WORK = (double*)malloc(n*sizeof(double));
905 CHECK(!WORK,MEM);
906 memcpy(rp,ap,m*n*sizeof(double));
907 integer res;
908 dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
909 CHECK(res,res);
910 free(WORK);
911 OK
912}
913
914/* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a,
915 integer *lda, doublecomplex *tau, doublecomplex *work, integer *info);
916
917int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
918 integer m = ar;
919 integer n = ac;
920 integer mn = MIN(m,n);
921 REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE);
922 DEBUGMSG("qr_l_C");
923 doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex));
924 CHECK(!WORK,MEM);
925 memcpy(rp,ap,m*n*sizeof(doublecomplex));
926 integer res;
927 zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res);
928 CHECK(res,res);
929 free(WORK);
930 OK
931}
932
933/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal *
934 a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
935 integer *info);
936
937int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) {
938 integer m = ar;
939 integer n = MIN(ac,ar);
940 integer k = taun;
941 DEBUGMSG("c_dorgqr");
942 integer lwork = 8*n; // FIXME
943 double *WORK = (double*)malloc(lwork*sizeof(double));
944 CHECK(!WORK,MEM);
945 memcpy(rp,ap,m*k*sizeof(double));
946 integer res;
947 dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res);
948 CHECK(res,res);
949 free(WORK);
950 OK
951}
952
953/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k,
954 doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
955 work, integer *lwork, integer *info);
956
957int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) {
958 integer m = ar;
959 integer n = MIN(ac,ar);
960 integer k = taun;
961 DEBUGMSG("z_ungqr");
962 integer lwork = 8*n; // FIXME
963 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
964 CHECK(!WORK,MEM);
965 memcpy(rp,ap,m*k*sizeof(doublecomplex));
966 integer res;
967 zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res);
968 CHECK(res,res);
969 free(WORK);
970 OK
971}
972
973
974//////////////////// Hessenberg factorization /////////////////////////
975
976/* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi,
977 doublereal *a, integer *lda, doublereal *tau, doublereal *work,
978 integer *lwork, integer *info);
979
980int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) {
981 integer m = ar;
982 integer n = ac;
983 integer mn = MIN(m,n);
984 REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE);
985 DEBUGMSG("hess_l_R");
986 integer lwork = 5*n; // fixme
987 double *WORK = (double*)malloc(lwork*sizeof(double));
988 CHECK(!WORK,MEM);
989 memcpy(rp,ap,m*n*sizeof(double));
990 integer res;
991 integer one = 1;
992 dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
993 CHECK(res,res);
994 free(WORK);
995 OK
996}
997
998
999/* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi,
1000 doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *
1001 work, integer *lwork, integer *info);
1002
1003int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) {
1004 integer m = ar;
1005 integer n = ac;
1006 integer mn = MIN(m,n);
1007 REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE);
1008 DEBUGMSG("hess_l_C");
1009 integer lwork = 5*n; // fixme
1010 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
1011 CHECK(!WORK,MEM);
1012 memcpy(rp,ap,m*n*sizeof(doublecomplex));
1013 integer res;
1014 integer one = 1;
1015 zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res);
1016 CHECK(res,res);
1017 free(WORK);
1018 OK
1019}
1020
1021//////////////////// Schur factorization /////////////////////////
1022
1023/* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n,
1024 doublereal *a, integer *lda, integer *sdim, doublereal *wr,
1025 doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work,
1026 integer *lwork, logical *bwork, integer *info);
1027
1028int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) {
1029 integer m = ar;
1030 integer n = ac;
1031 REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
1032 DEBUGMSG("schur_l_R");
1033 //int k;
1034 //printf("---------------------------\n");
1035 //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
1036 //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
1037 //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
1038 memcpy(sp,ap,n*n*sizeof(double));
1039 integer lwork = 6*n; // fixme
1040 double *WORK = (double*)malloc(lwork*sizeof(double));
1041 double *WR = (double*)malloc(n*sizeof(double));
1042 double *WI = (double*)malloc(n*sizeof(double));
1043 // WR and WI not really required in this call
1044 logical *BWORK = (logical*)malloc(n*sizeof(logical));
1045 integer res;
1046 integer sdim;
1047 dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res);
1048 //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n");
1049 //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n");
1050 //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n");
1051 if(res>0) {
1052 return NOCONVER;
1053 }
1054 CHECK(res,res);
1055 free(WR);
1056 free(WI);
1057 free(BWORK);
1058 free(WORK);
1059 OK
1060}
1061
1062
1063/* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n,
1064 doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w,
1065 doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork,
1066 doublereal *rwork, logical *bwork, integer *info);
1067
1068int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
1069 integer m = ar;
1070 integer n = ac;
1071 REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE);
1072 DEBUGMSG("schur_l_C");
1073 memcpy(sp,ap,n*n*sizeof(doublecomplex));
1074 integer lwork = 6*n; // fixme
1075 doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
1076 doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex));
1077 // W not really required in this call
1078 logical *BWORK = (logical*)malloc(n*sizeof(logical));
1079 double *RWORK = (double*)malloc(n*sizeof(double));
1080 integer res;
1081 integer sdim;
1082 zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W,
1083 up,&n,
1084 WORK,&lwork,RWORK,BWORK,&res);
1085 if(res>0) {
1086 return NOCONVER;
1087 }
1088 CHECK(res,res);
1089 free(W);
1090 free(BWORK);
1091 free(WORK);
1092 OK
1093}
1094
1095//////////////////// LU factorization /////////////////////////
1096
1097/* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer *
1098 lda, integer *ipiv, integer *info);
1099
1100int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) {
1101 integer m = ar;
1102 integer n = ac;
1103 integer mn = MIN(m,n);
1104 REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
1105 DEBUGMSG("lu_l_R");
1106 integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
1107 memcpy(rp,ap,m*n*sizeof(double));
1108 integer res;
1109 dgetrf_ (&m,&n,rp,&m,auxipiv,&res);
1110 if(res>0) {
1111 res = 0; // fixme
1112 }
1113 CHECK(res,res);
1114 int k;
1115 for (k=0; k<mn; k++) {
1116 ipivp[k] = auxipiv[k];
1117 }
1118 free(auxipiv);
1119 OK
1120}
1121
1122
1123/* Subroutine */ int zgetrf_(integer *m, integer *n, doublecomplex *a,
1124 integer *lda, integer *ipiv, integer *info);
1125
1126int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
1127 integer m = ar;
1128 integer n = ac;
1129 integer mn = MIN(m,n);
1130 REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
1131 DEBUGMSG("lu_l_C");
1132 integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
1133 memcpy(rp,ap,m*n*sizeof(doublecomplex));
1134 integer res;
1135 zgetrf_ (&m,&n,rp,&m,auxipiv,&res);
1136 if(res>0) {
1137 res = 0; // fixme
1138 }
1139 CHECK(res,res);
1140 int k;
1141 for (k=0; k<mn; k++) {
1142 ipivp[k] = auxipiv[k];
1143 }
1144 free(auxipiv);
1145 OK
1146}
1147
1148
1149//////////////////// LU substitution /////////////////////////
1150
1151/* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs,
1152 doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
1153 ldb, integer *info);
1154
1155int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
1156 integer m = ar;
1157 integer n = ac;
1158 integer mrhs = br;
1159 integer nrhs = bc;
1160
1161 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1162 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1163 int k;
1164 for (k=0; k<n; k++) {
1165 auxipiv[k] = (integer)ipivp[k];
1166 }
1167 integer res;
1168 memcpy(xp,bp,mrhs*nrhs*sizeof(double));
1169 dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,xp,&mrhs,&res);
1170 CHECK(res,res);
1171 free(auxipiv);
1172 OK
1173}
1174
1175
1176/* Subroutine */ int zgetrs_(char *trans, integer *n, integer *nrhs,
1177 doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b,
1178 integer *ldb, integer *info);
1179
1180int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x)) {
1181 integer m = ar;
1182 integer n = ac;
1183 integer mrhs = br;
1184 integer nrhs = bc;
1185
1186 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1187 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1188 int k;
1189 for (k=0; k<n; k++) {
1190 auxipiv[k] = (integer)ipivp[k];
1191 }
1192 integer res;
1193 memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex));
1194 zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res);
1195 CHECK(res,res);
1196 free(auxipiv);
1197 OK
1198}
1199
1200//////////////////// Matrix Product /////////////////////////
1201
1202void dgemm_(char *, char *, integer *, integer *, integer *,
1203 double *, const double *, integer *, const double *,
1204 integer *, double *, double *, integer *);
1205
1206int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) {
1207 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1208 DEBUGMSG("dgemm_");
1209 CHECKNANR(a,"NaN multR Input\n")
1210 CHECKNANR(b,"NaN multR Input\n")
1211 integer m = ta?ac:ar;
1212 integer n = tb?br:bc;
1213 integer k = ta?ar:ac;
1214 integer lda = ar;
1215 integer ldb = br;
1216 integer ldc = rr;
1217 double alpha = 1;
1218 double beta = 0;
1219 dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
1220 CHECKNANR(r,"NaN multR Output\n")
1221 OK
1222}
1223
1224void zgemm_(char *, char *, integer *, integer *, integer *,
1225 doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
1226 integer *, doublecomplex *, doublecomplex *, integer *);
1227
1228int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) {
1229 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1230 DEBUGMSG("zgemm_");
1231 CHECKNANC(a,"NaN multC Input\n")
1232 CHECKNANC(b,"NaN multC Input\n")
1233 integer m = ta?ac:ar;
1234 integer n = tb?br:bc;
1235 integer k = ta?ar:ac;
1236 integer lda = ar;
1237 integer ldb = br;
1238 integer ldc = rr;
1239 doublecomplex alpha = {1,0};
1240 doublecomplex beta = {0,0};
1241 zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
1242 ap,&lda,
1243 bp,&ldb,&beta,
1244 rp,&ldc);
1245 CHECKNANC(r,"NaN multC Output\n")
1246 OK
1247}
1248
1249void sgemm_(char *, char *, integer *, integer *, integer *,
1250 float *, const float *, integer *, const float *,
1251 integer *, float *, float *, integer *);
1252
1253int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)) {
1254 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1255 DEBUGMSG("sgemm_");
1256 integer m = ta?ac:ar;
1257 integer n = tb?br:bc;
1258 integer k = ta?ar:ac;
1259 integer lda = ar;
1260 integer ldb = br;
1261 integer ldc = rr;
1262 float alpha = 1;
1263 float beta = 0;
1264 sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc);
1265 OK
1266}
1267
1268void cgemm_(char *, char *, integer *, integer *, integer *,
1269 complex *, const complex *, integer *, const complex *,
1270 integer *, complex *, complex *, integer *);
1271
1272int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) {
1273 //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
1274 DEBUGMSG("cgemm_");
1275 integer m = ta?ac:ar;
1276 integer n = tb?br:bc;
1277 integer k = ta?ar:ac;
1278 integer lda = ar;
1279 integer ldb = br;
1280 integer ldc = rr;
1281 complex alpha = {1,0};
1282 complex beta = {0,0};
1283 cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,
1284 ap,&lda,
1285 bp,&ldb,&beta,
1286 rp,&ldc);
1287 OK
1288}
1289
1290//////////////////// transpose /////////////////////////
1291
1292int transF(KFMAT(x),FMAT(t)) {
1293 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1294 DEBUGMSG("transF");
1295 int i,j;
1296 for (i=0; i<tr; i++) {
1297 for (j=0; j<tc; j++) {
1298 tp[i*tc+j] = xp[j*xc+i];
1299 }
1300 }
1301 OK
1302}
1303
1304int transR(KDMAT(x),DMAT(t)) {
1305 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1306 DEBUGMSG("transR");
1307 int i,j;
1308 for (i=0; i<tr; i++) {
1309 for (j=0; j<tc; j++) {
1310 tp[i*tc+j] = xp[j*xc+i];
1311 }
1312 }
1313 OK
1314}
1315
1316int transQ(KQMAT(x),QMAT(t)) {
1317 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1318 DEBUGMSG("transQ");
1319 int i,j;
1320 for (i=0; i<tr; i++) {
1321 for (j=0; j<tc; j++) {
1322 tp[i*tc+j] = xp[j*xc+i];
1323 }
1324 }
1325 OK
1326}
1327
1328int transC(KCMAT(x),CMAT(t)) {
1329 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1330 DEBUGMSG("transC");
1331 int i,j;
1332 for (i=0; i<tr; i++) {
1333 for (j=0; j<tc; j++) {
1334 tp[i*tc+j] = xp[j*xc+i];
1335 }
1336 }
1337 OK
1338}
1339
1340int transP(KPMAT(x), PMAT(t)) {
1341 REQUIRES(xr==tc && xc==tr,BAD_SIZE);
1342 REQUIRES(xs==ts,NOCONVER);
1343 DEBUGMSG("transP");
1344 int i,j;
1345 for (i=0; i<tr; i++) {
1346 for (j=0; j<tc; j++) {
1347 memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs);
1348 }
1349 }
1350 OK
1351}
1352
1353//////////////////// constant /////////////////////////
1354
1355int constantF(float * pval, FVEC(r)) {
1356 DEBUGMSG("constantF")
1357 int k;
1358 double val = *pval;
1359 for(k=0;k<rn;k++) {
1360 rp[k]=val;
1361 }
1362 OK
1363}
1364
1365int constantR(double * pval, DVEC(r)) {
1366 DEBUGMSG("constantR")
1367 int k;
1368 double val = *pval;
1369 for(k=0;k<rn;k++) {
1370 rp[k]=val;
1371 }
1372 OK
1373}
1374
1375int constantQ(complex* pval, QVEC(r)) {
1376 DEBUGMSG("constantQ")
1377 int k;
1378 complex val = *pval;
1379 for(k=0;k<rn;k++) {
1380 rp[k]=val;
1381 }
1382 OK
1383}
1384
1385int constantC(doublecomplex* pval, CVEC(r)) {
1386 DEBUGMSG("constantC")
1387 int k;
1388 doublecomplex val = *pval;
1389 for(k=0;k<rn;k++) {
1390 rp[k]=val;
1391 }
1392 OK
1393}
1394
1395int constantP(void* pval, PVEC(r)) {
1396 DEBUGMSG("constantP")
1397 int k;
1398 for(k=0;k<rn;k++) {
1399 memcpy(rp+k*rs,pval,rs);
1400 }
1401 OK
1402}
1403
1404//////////////////// float-double conversion /////////////////////////
1405
1406int float2double(FVEC(x),DVEC(y)) {
1407 DEBUGMSG("float2double")
1408 int k;
1409 for(k=0;k<xn;k++) {
1410 yp[k]=xp[k];
1411 }
1412 OK
1413}
1414
1415int double2float(DVEC(x),FVEC(y)) {
1416 DEBUGMSG("double2float")
1417 int k;
1418 for(k=0;k<xn;k++) {
1419 yp[k]=xp[k];
1420 }
1421 OK
1422}
1423
1424//////////////////// conjugate /////////////////////////
1425
1426int conjugateQ(KQVEC(x),QVEC(t)) {
1427 REQUIRES(xn==tn,BAD_SIZE);
1428 DEBUGMSG("conjugateQ");
1429 int k;
1430 for(k=0;k<xn;k++) {
1431 tp[k].r = xp[k].r;
1432 tp[k].i = -xp[k].i;
1433 }
1434 OK
1435}
1436
1437int conjugateC(KCVEC(x),CVEC(t)) {
1438 REQUIRES(xn==tn,BAD_SIZE);
1439 DEBUGMSG("conjugateC");
1440 int k;
1441 for(k=0;k<xn;k++) {
1442 tp[k].r = xp[k].r;
1443 tp[k].i = -xp[k].i;
1444 }
1445 OK
1446}
1447
1448//////////////////// step /////////////////////////
1449
1450int stepF(FVEC(x),FVEC(y)) {
1451 DEBUGMSG("stepF")
1452 int k;
1453 for(k=0;k<xn;k++) {
1454 yp[k]=xp[k]>0;
1455 }
1456 OK
1457}
1458
1459int stepD(DVEC(x),DVEC(y)) {
1460 DEBUGMSG("stepD")
1461 int k;
1462 for(k=0;k<xn;k++) {
1463 yp[k]=xp[k]>0;
1464 }
1465 OK
1466}
1467
1468//////////////////// cond /////////////////////////
1469
1470int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {
1471 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
1472 DEBUGMSG("condF")
1473 int k;
1474 for(k=0;k<xn;k++) {
1475 rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);
1476 }
1477 OK
1478}
1479
1480int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {
1481 REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);
1482 DEBUGMSG("condD")
1483 int k;
1484 for(k=0;k<xn;k++) {
1485 rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);
1486 }
1487 OK
1488}
1489
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
new file mode 100644
index 0000000..a3f1899
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -0,0 +1,60 @@
1/*
2 * We have copied the definitions in f2c.h required
3 * to compile clapack.h, modified to support both
4 * 32 and 64 bit
5
6 http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib
7 http://www.ibm.com/developerworks/library/l-port64.html
8 */
9
10#ifdef _LP64
11typedef int integer;
12typedef unsigned int uinteger;
13typedef int logical;
14typedef long longint; /* system-dependent */
15typedef unsigned long ulongint; /* system-dependent */
16#else
17typedef long int integer;
18typedef unsigned long int uinteger;
19typedef long int logical;
20typedef long long longint; /* system-dependent */
21typedef unsigned long long ulongint; /* system-dependent */
22#endif
23
24typedef char *address;
25typedef short int shortint;
26typedef float real;
27typedef double doublereal;
28typedef struct { real r, i; } complex;
29typedef struct { doublereal r, i; } doublecomplex;
30typedef short int shortlogical;
31typedef char logical1;
32typedef char integer1;
33
34typedef logical (*L_fp)();
35typedef short ftnlen;
36
37/********************************************************/
38
39#define FVEC(A) int A##n, float*A##p
40#define DVEC(A) int A##n, double*A##p
41#define QVEC(A) int A##n, complex*A##p
42#define CVEC(A) int A##n, doublecomplex*A##p
43#define PVEC(A) int A##n, void* A##p, int A##s
44#define FMAT(A) int A##r, int A##c, float* A##p
45#define DMAT(A) int A##r, int A##c, double* A##p
46#define QMAT(A) int A##r, int A##c, complex* A##p
47#define CMAT(A) int A##r, int A##c, doublecomplex* A##p
48#define PMAT(A) int A##r, int A##c, void* A##p, int A##s
49
50#define KFVEC(A) int A##n, const float*A##p
51#define KDVEC(A) int A##n, const double*A##p
52#define KQVEC(A) int A##n, const complex*A##p
53#define KCVEC(A) int A##n, const doublecomplex*A##p
54#define KPVEC(A) int A##n, const void* A##p, int A##s
55#define KFMAT(A) int A##r, int A##c, const float* A##p
56#define KDMAT(A) int A##r, int A##c, const double* A##p
57#define KQMAT(A) int A##r, int A##c, const complex* A##p
58#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p
59#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s
60
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs
new file mode 100644
index 0000000..7d134bf
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs
@@ -0,0 +1,295 @@
1{-# LANGUAGE FlexibleContexts #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.LinearAlgebra.Util
5Copyright : (c) Alberto Ruiz 2013
6License : GPL
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10
11-}
12-----------------------------------------------------------------------------
13{-# OPTIONS_HADDOCK hide #-}
14
15module Numeric.LinearAlgebra.Util(
16
17 -- * Convenience functions
18 size, disp,
19 zeros, ones,
20 diagl,
21 row,
22 col,
23 (&), (¦), (——), (#),
24 (?), (¿),
25 rand, randn,
26 cross,
27 norm,
28 unitary,
29 mt,
30 pairwiseD2,
31 rowOuters,
32 null1,
33 null1sym,
34 -- * Convolution
35 -- ** 1D
36 corr, conv, corrMin,
37 -- ** 2D
38 corr2, conv2, separable,
39 -- * Tools for the Kronecker product
40 --
41 -- | (see A. Fusiello, A matter of notation: Several uses of the Kronecker product in
42 -- 3d computer vision, Pattern Recognition Letters 28 (15) (2007) 2127-2132)
43
44 --
45 -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@
46 vec,
47 vech,
48 dup,
49 vtrans,
50 -- * Plot
51 mplot,
52 plot, parametricPlot,
53 splot, mesh, meshdom,
54 matrixToPGM, imshow,
55 gnuplotX, gnuplotpdf, gnuplotWin
56) where
57
58import Numeric.Container
59import Numeric.LinearAlgebra.Algorithms hiding (i)
60import Numeric.Matrix()
61import Numeric.Vector()
62
63import System.Random(randomIO)
64import Numeric.LinearAlgebra.Util.Convolution
65import Graphics.Plot
66
67
68{- | print a real matrix with given number of digits after the decimal point
69
70>>> disp 5 $ ident 2 / 3
712x2
720.33333 0.00000
730.00000 0.33333
74
75-}
76disp :: Int -> Matrix Double -> IO ()
77
78disp n = putStrLn . dispf n
79
80-- | pseudorandom matrix with uniform elements between 0 and 1
81randm :: RandDist
82 -> Int -- ^ rows
83 -> Int -- ^ columns
84 -> IO (Matrix Double)
85randm d r c = do
86 seed <- randomIO
87 return (reshape c $ randomVector seed d (r*c))
88
89-- | pseudorandom matrix with uniform elements between 0 and 1
90rand :: Int -> Int -> IO (Matrix Double)
91rand = randm Uniform
92
93{- | pseudorandom matrix with normal elements
94
95>>> x <- randn 3 5
96>>> disp 3 x
973x5
980.386 -1.141 0.491 -0.510 1.512
990.069 -0.919 1.022 -0.181 0.745
1000.313 -0.670 -0.097 -1.575 -0.583
101
102-}
103randn :: Int -> Int -> IO (Matrix Double)
104randn = randm Gaussian
105
106{- | create a real diagonal matrix from a list
107
108>>> diagl [1,2,3]
109(3><3)
110 [ 1.0, 0.0, 0.0
111 , 0.0, 2.0, 0.0
112 , 0.0, 0.0, 3.0 ]
113
114-}
115diagl :: [Double] -> Matrix Double
116diagl = diag . fromList
117
118-- | a real matrix of zeros
119zeros :: Int -- ^ rows
120 -> Int -- ^ columns
121 -> Matrix Double
122zeros r c = konst 0 (r,c)
123
124-- | a real matrix of ones
125ones :: Int -- ^ rows
126 -> Int -- ^ columns
127 -> Matrix Double
128ones r c = konst 1 (r,c)
129
130-- | concatenation of real vectors
131infixl 3 &
132(&) :: Vector Double -> Vector Double -> Vector Double
133a & b = vjoin [a,b]
134
135{- | horizontal concatenation of real matrices
136
137 (unicode 0x00a6, broken bar)
138
139>>> ident 3 ¦ konst 7 (3,4)
140(3><7)
141 [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
142 , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
143 , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]
144
145-}
146infixl 3 ¦
147(¦) :: Matrix Double -> Matrix Double -> Matrix Double
148a ¦ b = fromBlocks [[a,b]]
149
150-- | vertical concatenation of real matrices
151--
152-- (unicode 0x2014, em dash)
153(——) :: Matrix Double -> Matrix Double -> Matrix Double
154infixl 2 ——
155a —— b = fromBlocks [[a],[b]]
156
157(#) :: Matrix Double -> Matrix Double -> Matrix Double
158infixl 2 #
159a # b = fromBlocks [[a],[b]]
160
161-- | create a single row real matrix from a list
162row :: [Double] -> Matrix Double
163row = asRow . fromList
164
165-- | create a single column real matrix from a list
166col :: [Double] -> Matrix Double
167col = asColumn . fromList
168
169{- | extract rows
170
171>>> (20><4) [1..] ? [2,1,1]
172(3><4)
173 [ 9.0, 10.0, 11.0, 12.0
174 , 5.0, 6.0, 7.0, 8.0
175 , 5.0, 6.0, 7.0, 8.0 ]
176
177-}
178infixl 9 ?
179(?) :: Element t => Matrix t -> [Int] -> Matrix t
180(?) = flip extractRows
181
182{- | extract columns
183
184(unicode 0x00bf, inverted question mark, Alt-Gr ?)
185
186>>> (3><4) [1..] ¿ [3,0]
187(3><2)
188 [ 4.0, 1.0
189 , 8.0, 5.0
190 , 12.0, 9.0 ]
191
192-}
193infixl 9 ¿
194(¿) :: Element t => Matrix t -> [Int] -> Matrix t
195(¿)= flip extractColumns
196
197
198cross :: Vector Double -> Vector Double -> Vector Double
199-- ^ cross product (for three-element real vectors)
200cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
201 | otherwise = error $ "cross ("++show x++") ("++show y++")"
202 where
203 [x1,x2,x3] = toList x
204 [y1,y2,y3] = toList y
205 z1 = x2*y3-x3*y2
206 z2 = x3*y1-x1*y3
207 z3 = x1*y2-x2*y1
208
209norm :: Vector Double -> Double
210-- ^ 2-norm of real vector
211norm = pnorm PNorm2
212
213
214-- | Obtains a vector in the same direction with 2-norm=1
215unitary :: Vector Double -> Vector Double
216unitary v = v / scalar (norm v)
217
218-- | ('rows' &&& 'cols')
219size :: Matrix t -> (Int, Int)
220size m = (rows m, cols m)
221
222-- | trans . inv
223mt :: Matrix Double -> Matrix Double
224mt = trans . inv
225
226----------------------------------------------------------------------
227
228-- | Matrix of pairwise squared distances of row vectors
229-- (using the matrix product trick in blog.smola.org)
230pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
231pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
232 | otherwise = error $ "pairwiseD2 with different number of columns: "
233 ++ show (size x) ++ ", " ++ show (size y)
234 where
235 ox = one (rows x)
236 oy = one (rows y)
237 oc = one (cols x)
238 one k = constant 1 k
239 x2 = x * x <> oc
240 y2 = y * y <> oc
241 ok = cols x == cols y
242
243--------------------------------------------------------------------------------
244
245-- | outer products of rows
246rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
247rowOuters a b = a' * b'
248 where
249 a' = kronecker a (ones 1 (cols b))
250 b' = kronecker (ones 1 (cols a)) b
251
252--------------------------------------------------------------------------------
253
254-- | solution of overconstrained homogeneous linear system
255null1 :: Matrix Double -> Vector Double
256null1 = last . toColumns . snd . rightSV
257
258-- | solution of overconstrained homogeneous symmetric linear system
259null1sym :: Matrix Double -> Vector Double
260null1sym = last . toColumns . snd . eigSH'
261
262--------------------------------------------------------------------------------
263
264vec :: Element t => Matrix t -> Vector t
265-- ^ stacking of columns
266vec = flatten . trans
267
268
269vech :: Element t => Matrix t -> Vector t
270-- ^ half-vectorization (of the lower triangular part)
271vech m = vjoin . zipWith f [0..] . toColumns $ m
272 where
273 f k v = subVector k (dim v - k) v
274
275
276dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t
277-- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k)
278dup k = trans $ fromRows $ map f es
279 where
280 rs = zip [0..] (toRows (ident (k^(2::Int))))
281 es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ]
282 f (i,j) | i == j = g (k*j + i)
283 | otherwise = g (k*j + i) + g (k*i + j)
284 g j = v
285 where
286 Just v = lookup j rs
287
288
289vtrans :: Element t => Int -> Matrix t -> Matrix t
290-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@
291vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m
292 | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows"
293 where
294 (q,r) = divMod (rows m) p
295
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs
new file mode 100644
index 0000000..82de476
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -0,0 +1,114 @@
1{-# LANGUAGE FlexibleContexts #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.LinearAlgebra.Util.Convolution
5Copyright : (c) Alberto Ruiz 2012
6License : GPL
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10
11-}
12-----------------------------------------------------------------------------
13
14module Numeric.LinearAlgebra.Util.Convolution(
15 corr, conv, corrMin,
16 corr2, conv2, separable
17) where
18
19import Numeric.LinearAlgebra
20
21
22vectSS :: Element t => Int -> Vector t -> Matrix t
23vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ]
24
25
26corr :: Product t => Vector t -- ^ kernel
27 -> Vector t -- ^ source
28 -> Vector t
29{- ^ correlation
30
31>>> corr (fromList[1,2,3]) (fromList [1..10])
32fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]
33
34-}
35corr ker v | dim ker <= dim v = vectSS (dim ker) v <> ker
36 | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")"
37
38
39conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t
40{- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product)
41
42>>> conv (fromList[1,1]) (fromList [-1,1])
43fromList [-1.0,0.0,1.0]
44
45-}
46conv ker v = corr ker' v'
47 where
48 ker' = (flatten.fliprl.asRow) ker
49 v' | dim ker > 1 = vjoin [z,v,z]
50 | otherwise = v
51 z = constant 0 (dim ker -1)
52
53corrMin :: (Container Vector t, RealElement t, Product t)
54 => Vector t
55 -> Vector t
56 -> Vector t
57-- ^ similar to 'corr', using 'min' instead of (*)
58corrMin ker v = minEvery ss (asRow ker) <> ones
59 where
60 minEvery a b = cond a b a a b
61 ss = vectSS (dim ker) v
62 ones = konst 1 (dim ker)
63
64
65
66matSS :: Element t => Int -> Matrix t -> [Matrix t]
67matSS dr m = map (reshape c) [ subVector (k*c) n v | k <- [0 .. r - dr] ]
68 where
69 v = flatten m
70 c = cols m
71 r = rows m
72 n = dr*c
73
74
75corr2 :: Product a => Matrix a -> Matrix a -> Matrix a
76-- ^ 2D correlation
77corr2 ker mat = dims
78 . concatMap (map (udot ker' . flatten) . matSS c . trans)
79 . matSS r $ mat
80 where
81 r = rows ker
82 c = cols ker
83 ker' = flatten (trans ker)
84 rr = rows mat - r + 1
85 rc = cols mat - c + 1
86 dims | rr > 0 && rc > 0 = (rr >< rc)
87 | otherwise = error $ "corr2: dim kernel ("++sz ker++") > dim matrix ("++sz mat++")"
88 sz m = show (rows m)++"x"++show (cols m)
89
90conv2 :: (Num a, Product a, Container Vector a) => Matrix a -> Matrix a -> Matrix a
91-- ^ 2D convolution
92conv2 k m = corr2 (fliprl . flipud $ k) pm
93 where
94 pm | r == 0 && c == 0 = m
95 | r == 0 = fromBlocks [[z3,m,z3]]
96 | c == 0 = fromBlocks [[z2],[m],[z2]]
97 | otherwise = fromBlocks [[z1,z2,z1]
98 ,[z3, m,z3]
99 ,[z1,z2,z1]]
100 r = rows k - 1
101 c = cols k - 1
102 h = rows m
103 w = cols m
104 z1 = konst 0 (r,c)
105 z2 = konst 0 (r,w)
106 z3 = konst 0 (h,c)
107
108-- TODO: could be simplified using future empty arrays
109
110
111separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t
112-- ^ matrix computation implemented as separated vector operations by rows and columns.
113separable f = fromColumns . map f . toColumns . fromRows . map f . toRows
114
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-20 11:57:51 +0000