summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--examples/tests.hs9
-rw-r--r--lib/Data/Packed/Internal/Common.hs1
-rw-r--r--lib/LAPACK.hs25
-rw-r--r--lib/LAPACK/lapack-aux.c42
-rw-r--r--lib/LAPACK/lapack-aux.h4
-rw-r--r--lib/LinearAlgebra/Algorithms.hs6
6 files changed, 82 insertions, 5 deletions
diff --git a/examples/tests.hs b/examples/tests.hs
index 5c14f96..80070fd 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -266,7 +266,8 @@ gammaTest = test "gamma" (gamma 5 == 24.0)
266 266
267--------------------------------------------------------------------- 267---------------------------------------------------------------------
268 268
269cholRTest = chol ((2><2) [1,2,2,9::Double]) == (2><2) [1,0,2,2.23606797749979] 269cholRTest = chol ((2><2) [1,2,2,9::Double]) == (2><2) [1,2,0,2.23606797749979]
270cholCTest = chol ((2><2) [1,2,2,9::Complex Double]) == (2><2) [1,2,0,2.23606797749979]
270 271
271--------------------------------------------------------------------- 272---------------------------------------------------------------------
272 273
@@ -316,6 +317,11 @@ tests = do
316 putStrLn "--------- pinv ------" 317 putStrLn "--------- pinv ------"
317 quickCheck (pinvTest . sqm :: SqM Double -> Bool) 318 quickCheck (pinvTest . sqm :: SqM Double -> Bool)
318 quickCheck (pinvTest . sqm :: SqM (Complex Double) -> Bool) 319 quickCheck (pinvTest . sqm :: SqM (Complex Double) -> Bool)
320 putStrLn "--------- chol ------"
321 runTestTT $ TestList
322 [ test "cholR" cholRTest
323 , test "cholC" cholRTest
324 ]
319 putStrLn "--------- nullspace ------" 325 putStrLn "--------- nullspace ------"
320 quickCheck (nullspaceTest :: RM -> Bool) 326 quickCheck (nullspaceTest :: RM -> Bool)
321 quickCheck (nullspaceTest :: CM -> Bool) 327 quickCheck (nullspaceTest :: CM -> Bool)
@@ -339,7 +345,6 @@ tests = do
339 , integrateTest 345 , integrateTest
340 , polySolveTest 346 , polySolveTest
341 , test "det" detTest 347 , test "det" detTest
342 , test "cholR" cholRTest
343 ] 348 ]
344 349
345bigtests = do 350bigtests = do
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index 218fb6b..75a5b1e 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -69,6 +69,7 @@ errorCode 2002 = "memory problem"
69errorCode 2003 = "bad file" 69errorCode 2003 = "bad file"
70errorCode 2004 = "singular" 70errorCode 2004 = "singular"
71errorCode 2005 = "didn't converge" 71errorCode 2005 = "didn't converge"
72errorCode 2006 = "the input matrix is not positive definite"
72errorCode n = "code "++show n 73errorCode n = "code "++show n
73 74
74{- | conversion of Haskell functions into function pointers that can be used in the C side 75{- | conversion of Haskell functions into function pointers that can be used in the C side
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index 54eea8a..e84647b 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -19,6 +19,7 @@ module LAPACK (
19 linearSolveR, linearSolveC, 19 linearSolveR, linearSolveC,
20 linearSolveLSR, linearSolveLSC, 20 linearSolveLSR, linearSolveLSC,
21 linearSolveSVDR, linearSolveSVDC, 21 linearSolveSVDR, linearSolveSVDC,
22 cholS, cholH
22) where 23) where
23 24
24import Data.Packed.Internal 25import Data.Packed.Internal
@@ -279,3 +280,27 @@ linearSolveSVDC_l rcond a b = unsafePerformIO $ do
279 n = cols a 280 n = cols a
280 nrhs = cols b 281 nrhs = cols b
281 282
283-----------------------------------------------------------------------------------
284foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM
285
286-- | Wrapper for LAPACK's /zpotrf/,which computes the Cholesky factorization of a
287-- complex Hermitian positive definite matrix.
288cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
289cholH a = unsafePerformIO $ do
290 r <- createMatrix ColumnMajor n n
291 zpotrf // mat fdat a // mat dat r // check "cholH" [fdat a]
292 return r
293 where n = rows a
294
295-----------------------------------------------------------------------------------
296foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
297
298-- | Wrapper for LAPACK's /dpotrf/,which computes the Cholesky factorization of a
299-- real symmetric positive definite matrix.
300cholS :: Matrix Double -> Matrix Double
301cholS a = unsafePerformIO $ do
302 r <- createMatrix ColumnMajor n n
303 dpotrf // mat fdat a // mat dat r // check "cholS" [fdat a]
304 return r
305 where n = rows a
306
diff --git a/lib/LAPACK/lapack-aux.c b/lib/LAPACK/lapack-aux.c
index c7d6d86..5781cd1 100644
--- a/lib/LAPACK/lapack-aux.c
+++ b/lib/LAPACK/lapack-aux.c
@@ -28,6 +28,7 @@
28#define BAD_FILE 2003 28#define BAD_FILE 2003
29#define SINGULAR 2004 29#define SINGULAR 2004
30#define NOCONVER 2005 30#define NOCONVER 2005
31#define NODEFPOS 2006
31 32
32//////////////////// real svd //////////////////////////////////// 33//////////////////// real svd ////////////////////////////////////
33 34
@@ -577,3 +578,44 @@ int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
577 free(AC); 578 free(AC);
578 OK 579 OK
579} 580}
581
582//////////////////// Cholesky factorization /////////////////////////
583
584int chol_l_H(KCMAT(a),CMAT(l)) {
585 integer n = ar;
586 REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
587 DEBUGMSG("chol_l_H");
588 memcpy(lp,ap,n*n*sizeof(doublecomplex));
589 char uplo = 'U';
590 integer res;
591 zpotrf_ (&uplo,&n,(doublecomplex*)lp,&n,&res);
592 CHECK(res>0,NODEFPOS);
593 CHECK(res,res);
594 doublecomplex zero = {0.,0.};
595 int r,c;
596 for (r=0; r<lr-1; r++) {
597 for(c=r+1; c<lc; c++) {
598 ((doublecomplex*)lp)[r*lc+c] = zero;
599 }
600 }
601 OK
602}
603
604int chol_l_S(KDMAT(a),DMAT(l)) {
605 integer n = ar;
606 REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
607 DEBUGMSG("chol_l_S");
608 memcpy(lp,ap,n*n*sizeof(double));
609 char uplo = 'U';
610 integer res;
611 dpotrf_ (&uplo,&n,lp,&n,&res);
612 CHECK(res>0,NODEFPOS);
613 CHECK(res,res);
614 int r,c;
615 for (r=0; r<lr-1; r++) {
616 for(c=r+1; c<lc; c++) {
617 lp[r*lc+c] = 0.;
618 }
619 }
620 OK
621}
diff --git a/lib/LAPACK/lapack-aux.h b/lib/LAPACK/lapack-aux.h
index 869eabb..c5bf3cc 100644
--- a/lib/LAPACK/lapack-aux.h
+++ b/lib/LAPACK/lapack-aux.h
@@ -36,3 +36,7 @@ int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x));
36int linearSolveSVDR_l(double,KDMAT(a),KDMAT(b),DMAT(x)); 36int linearSolveSVDR_l(double,KDMAT(a),KDMAT(b),DMAT(x));
37 37
38int linearSolveSVDC_l(double,KCMAT(a),KCMAT(b),CMAT(x)); 38int linearSolveSVDC_l(double,KCMAT(a),KCMAT(b),CMAT(x));
39
40int chol_l_H(KCMAT(a),CMAT(r));
41
42int chol_l_S(KDMAT(a),DMAT(r));
diff --git a/lib/LinearAlgebra/Algorithms.hs b/lib/LinearAlgebra/Algorithms.hs
index a67f822..007067f 100644
--- a/lib/LinearAlgebra/Algorithms.hs
+++ b/lib/LinearAlgebra/Algorithms.hs
@@ -44,7 +44,7 @@ module LinearAlgebra.Algorithms (
44 44
45import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj) 45import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
46import Data.Packed 46import Data.Packed
47import GSL.Matrix 47import GSL.Matrix(luR,luC,qr)
48import GSL.Vector 48import GSL.Vector
49import LAPACK 49import LAPACK
50import Complex 50import Complex
@@ -70,7 +70,7 @@ instance GenMat Double where
70 ctrans = trans 70 ctrans = trans
71 eig = eigR 71 eig = eigR
72 eigSH = LAPACK.eigS 72 eigSH = LAPACK.eigS
73 chol = cholR 73 chol = cholS
74 74
75instance GenMat (Complex Double) where 75instance GenMat (Complex Double) where
76 svd = svdC 76 svd = svdC
@@ -80,7 +80,7 @@ instance GenMat (Complex Double) where
80 ctrans = conjTrans 80 ctrans = conjTrans
81 eig = eigC 81 eig = eigC
82 eigSH = LAPACK.eigH 82 eigSH = LAPACK.eigH
83 chol = error "cholC not yet implemented" -- waiting for GSL-1.10 83 chol = cholH
84 84
85-- | eigensystem of a symmetric matrix 85-- | eigensystem of a symmetric matrix
86eigS :: Matrix Double -> (Vector Double, Matrix Double) 86eigS :: Matrix Double -> (Vector Double, Matrix Double)