diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-06 13:10:08 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-06 13:10:08 +0200 |
commit | c04b342324001dc74baaa5e74264e61a76937f88 (patch) | |
tree | 1155bdcce361167dc62f76b537fb0f34b8e6aab5 /lib | |
parent | 8ecec15cee88462e46ec9d0ce361224c0fcdba31 (diff) |
added qrgr
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Numeric/HMatrix.hs | 2 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 21 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK.hs | 26 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 41 |
4 files changed, 80 insertions, 10 deletions
diff --git a/lib/Numeric/HMatrix.hs b/lib/Numeric/HMatrix.hs index f49ea53..09ad518 100644 --- a/lib/Numeric/HMatrix.hs +++ b/lib/Numeric/HMatrix.hs | |||
@@ -86,7 +86,7 @@ module Numeric.HMatrix ( | |||
86 | geigSH', | 86 | geigSH', |
87 | 87 | ||
88 | -- * QR | 88 | -- * QR |
89 | qr, rq, | 89 | qr, rq, qrRaw, qrgr, |
90 | 90 | ||
91 | -- * Cholesky | 91 | -- * Cholesky |
92 | chol, cholSH, mbCholSH, | 92 | chol, cholSH, mbCholSH, |
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index de9cb00..8c4b610 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -49,7 +49,7 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
49 | eigenvalues, eigenvaluesSH, eigenvaluesSH', | 49 | eigenvalues, eigenvaluesSH, eigenvaluesSH', |
50 | geigSH', | 50 | geigSH', |
51 | -- ** QR | 51 | -- ** QR |
52 | qr, rq, | 52 | qr, rq, qrRaw, qrgr, |
53 | -- ** Cholesky | 53 | -- ** Cholesky |
54 | chol, cholSH, mbCholSH, | 54 | chol, cholSH, mbCholSH, |
55 | -- ** Hessenberg | 55 | -- ** Hessenberg |
@@ -115,7 +115,8 @@ class (Product t, | |||
115 | eigOnlySH :: Matrix t -> Vector Double | 115 | eigOnlySH :: Matrix t -> Vector Double |
116 | cholSH' :: Matrix t -> Matrix t | 116 | cholSH' :: Matrix t -> Matrix t |
117 | mbCholSH' :: Matrix t -> Maybe (Matrix t) | 117 | mbCholSH' :: Matrix t -> Maybe (Matrix t) |
118 | qr' :: Matrix t -> (Matrix t, Matrix t) | 118 | qr' :: Matrix t -> (Matrix t, Vector t) |
119 | qrgr' :: Int -> (Matrix t, Vector t) -> Matrix t | ||
119 | hess' :: Matrix t -> (Matrix t, Matrix t) | 120 | hess' :: Matrix t -> (Matrix t, Matrix t) |
120 | schur' :: Matrix t -> (Matrix t, Matrix t) | 121 | schur' :: Matrix t -> (Matrix t, Matrix t) |
121 | 122 | ||
@@ -136,7 +137,8 @@ instance Field Double where | |||
136 | eigOnlySH = eigOnlyS | 137 | eigOnlySH = eigOnlyS |
137 | cholSH' = cholS | 138 | cholSH' = cholS |
138 | mbCholSH' = mbCholS | 139 | mbCholSH' = mbCholS |
139 | qr' = unpackQR . qrR | 140 | qr' = qrR |
141 | qrgr' = qrgrR | ||
140 | hess' = unpackHess hessR | 142 | hess' = unpackHess hessR |
141 | schur' = schurR | 143 | schur' = schurR |
142 | 144 | ||
@@ -161,7 +163,8 @@ instance Field (Complex Double) where | |||
161 | eigOnlySH = eigOnlyH | 163 | eigOnlySH = eigOnlyH |
162 | cholSH' = cholH | 164 | cholSH' = cholH |
163 | mbCholSH' = mbCholH | 165 | mbCholSH' = mbCholH |
164 | qr' = unpackQR . qrC | 166 | qr' = qrC |
167 | qrgr' = qrgrC | ||
165 | hess' = unpackHess hessC | 168 | hess' = unpackHess hessC |
166 | schur' = schurC | 169 | schur' = schurC |
167 | 170 | ||
@@ -285,7 +288,15 @@ eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m | |||
285 | -- | 288 | -- |
286 | -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. | 289 | -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular. |
287 | qr :: Field t => Matrix t -> (Matrix t, Matrix t) | 290 | qr :: Field t => Matrix t -> (Matrix t, Matrix t) |
288 | qr = {-# SCC "qr" #-} qr' | 291 | qr = {-# SCC "qr" #-} unpackQR . qr' |
292 | |||
293 | qrRaw m = qr' m | ||
294 | |||
295 | {- | generate a matrix with k orthogonal columns from the output of qrRaw | ||
296 | -} | ||
297 | qrgr 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) | ||
289 | 300 | ||
290 | -- | RQ factorization. | 301 | -- | RQ factorization. |
291 | -- | 302 | -- |
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs index a12f6d8..11394a6 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs | |||
@@ -35,7 +35,7 @@ module Numeric.LinearAlgebra.LAPACK ( | |||
35 | -- * Cholesky | 35 | -- * Cholesky |
36 | cholS, cholH, mbCholS, mbCholH, | 36 | cholS, cholH, mbCholS, mbCholH, |
37 | -- * QR | 37 | -- * QR |
38 | qrR, qrC, | 38 | qrR, qrC, qrgrR, qrgrC, |
39 | -- * Hessenberg | 39 | -- * Hessenberg |
40 | hessR, hessC, | 40 | hessR, hessC, |
41 | -- * Schur | 41 | -- * Schur |
@@ -444,9 +444,27 @@ qrAux f st a = unsafePerformIO $ do | |||
444 | tau <- createVector mn | 444 | tau <- createVector mn |
445 | app3 f mat a vec tau mat r st | 445 | app3 f mat a vec tau mat r st |
446 | return (r,tau) | 446 | return (r,tau) |
447 | where m = rows a | 447 | where |
448 | n = cols a | 448 | m = rows a |
449 | mn = min m n | 449 | n = cols a |
450 | mn = min m n | ||
451 | |||
452 | foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM | ||
453 | foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM | ||
454 | |||
455 | -- | build rotation from reflectors | ||
456 | qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double | ||
457 | qrgrR = qrgrAux dorgqr "qrgrR" | ||
458 | -- | build rotation from reflectors | ||
459 | qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double) | ||
460 | qrgrC = qrgrAux zungqr "qrgrC" | ||
461 | |||
462 | qrgrAux 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] | ||
450 | 468 | ||
451 | ----------------------------------------------------------------------------------- | 469 | ----------------------------------------------------------------------------------- |
452 | foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM | 470 | foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM |
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c index 51972eb..e5e45ef 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | |||
@@ -930,6 +930,47 @@ int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { | |||
930 | OK | 930 | OK |
931 | } | 931 | } |
932 | 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 | |||
937 | int 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 | |||
957 | int 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 | |||
933 | //////////////////// Hessenberg factorization ///////////////////////// | 974 | //////////////////// Hessenberg factorization ///////////////////////// |
934 | 975 | ||
935 | /* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi, | 976 | /* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi, |