summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Numeric/HMatrix.hs2
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs21
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK.hs26
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c41
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.
287qr :: Field t => Matrix t -> (Matrix t, Matrix t) 290qr :: Field t => Matrix t -> (Matrix t, Matrix t)
288qr = {-# SCC "qr" #-} qr' 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)
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
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]
450 468
451----------------------------------------------------------------------------------- 469-----------------------------------------------------------------------------------
452foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM 470foreign 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
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
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,