diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/LAPACK.hs | 31 | ||||
-rw-r--r-- | lib/LAPACK/lapack-aux.c | 35 | ||||
-rw-r--r-- | lib/LAPACK/lapack-aux.h | 4 |
3 files changed, 69 insertions, 1 deletions
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs index e84647b..43a313b 100644 --- a/lib/LAPACK.hs +++ b/lib/LAPACK.hs | |||
@@ -19,7 +19,8 @@ module LAPACK ( | |||
19 | linearSolveR, linearSolveC, | 19 | linearSolveR, linearSolveC, |
20 | linearSolveLSR, linearSolveLSC, | 20 | linearSolveLSR, linearSolveLSC, |
21 | linearSolveSVDR, linearSolveSVDC, | 21 | linearSolveSVDR, linearSolveSVDC, |
22 | cholS, cholH | 22 | cholS, cholH, |
23 | qrR, qrC | ||
23 | ) where | 24 | ) where |
24 | 25 | ||
25 | import Data.Packed.Internal | 26 | import Data.Packed.Internal |
@@ -304,3 +305,31 @@ cholS a = unsafePerformIO $ do | |||
304 | return r | 305 | return r |
305 | where n = rows a | 306 | where n = rows a |
306 | 307 | ||
308 | ----------------------------------------------------------------------------------- | ||
309 | foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM | ||
310 | |||
311 | -- | Wrapper for LAPACK's /dgeqr2/,which computes a QR factorization of a real matrix. | ||
312 | qrR :: Matrix Double -> (Matrix Double, Vector Double) | ||
313 | qrR a = unsafePerformIO $ do | ||
314 | r <- createMatrix ColumnMajor m n | ||
315 | tau <- createVector mn | ||
316 | dgeqr2 // mat fdat a // vec tau // mat dat r // check "qrR" [fdat a] | ||
317 | return (r,tau) | ||
318 | where m = rows a | ||
319 | n = cols a | ||
320 | mn = min m n | ||
321 | |||
322 | ----------------------------------------------------------------------------------- | ||
323 | foreign import ccall "LAPACK/lapack-aux.h qr_l_C" zgeqr2 :: TCMCVCM | ||
324 | |||
325 | -- | Wrapper for LAPACK's /zgeqr2/,which computes a QR factorization of a complex matrix. | ||
326 | qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) | ||
327 | qrC a = unsafePerformIO $ do | ||
328 | r <- createMatrix ColumnMajor m n | ||
329 | tau <- createVector mn | ||
330 | zgeqr2 // mat fdat a // vec tau // mat dat r // check "qrC" [fdat a] | ||
331 | return (r,tau) | ||
332 | where m = rows a | ||
333 | n = cols a | ||
334 | mn = min m n | ||
335 | |||
diff --git a/lib/LAPACK/lapack-aux.c b/lib/LAPACK/lapack-aux.c index 5781cd1..4ef9a6e 100644 --- a/lib/LAPACK/lapack-aux.c +++ b/lib/LAPACK/lapack-aux.c | |||
@@ -619,3 +619,38 @@ int chol_l_S(KDMAT(a),DMAT(l)) { | |||
619 | } | 619 | } |
620 | OK | 620 | OK |
621 | } | 621 | } |
622 | |||
623 | //////////////////// QR factorization ///////////////////////// | ||
624 | // TO DO: unpack | ||
625 | |||
626 | int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) { | ||
627 | integer m = ar; | ||
628 | integer n = ac; | ||
629 | integer mn = MIN(m,n); | ||
630 | REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); | ||
631 | DEBUGMSG("qr_l_R"); | ||
632 | double *WORK = (double*)malloc(m*sizeof(double)); | ||
633 | CHECK(!WORK,MEM); | ||
634 | memcpy(rp,ap,m*n*sizeof(double)); | ||
635 | integer res; | ||
636 | dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); | ||
637 | CHECK(res,res); | ||
638 | free(WORK); | ||
639 | OK | ||
640 | } | ||
641 | |||
642 | int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { | ||
643 | integer m = ar; | ||
644 | integer n = ac; | ||
645 | integer mn = MIN(m,n); | ||
646 | REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); | ||
647 | DEBUGMSG("qr_l_C"); | ||
648 | doublecomplex *WORK = (doublecomplex*)malloc(m*sizeof(doublecomplex)); | ||
649 | CHECK(!WORK,MEM); | ||
650 | memcpy(rp,ap,m*n*sizeof(doublecomplex)); | ||
651 | integer res; | ||
652 | zgeqr2_ (&m,&n,(doublecomplex*)rp,&m,(doublecomplex*)taup,WORK,&res); | ||
653 | CHECK(res,res); | ||
654 | free(WORK); | ||
655 | OK | ||
656 | } | ||
diff --git a/lib/LAPACK/lapack-aux.h b/lib/LAPACK/lapack-aux.h index c5bf3cc..ea71847 100644 --- a/lib/LAPACK/lapack-aux.h +++ b/lib/LAPACK/lapack-aux.h | |||
@@ -40,3 +40,7 @@ int linearSolveSVDC_l(double,KCMAT(a),KCMAT(b),CMAT(x)); | |||
40 | int chol_l_H(KCMAT(a),CMAT(r)); | 40 | int chol_l_H(KCMAT(a),CMAT(r)); |
41 | 41 | ||
42 | int chol_l_S(KDMAT(a),DMAT(r)); | 42 | int chol_l_S(KDMAT(a),DMAT(r)); |
43 | |||
44 | int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)); | ||
45 | |||
46 | int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)); | ||