mbCholSH

3 files changed, 35 insertions, 12 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index aa18d37..60f5971 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -46,7 +46,7 @@ module Numeric.LinearAlgebra.Algorithms (
 -- ** QR
     qr, rq,
 -- ** Cholesky
-    chol, cholSH,
+    chol, cholSH, mbCholSH,
 -- ** Hessenberg
     hess,
 -- ** Schur
@@ -100,6 +100,7 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     eigOnly      :: Matrix t -> Vector (Complex Double)
     eigOnlySH    :: Matrix t -> Vector Double
     cholSH'      :: Matrix t -> Matrix t
+    mbCholSH'    :: Matrix t -> Maybe (Matrix t)
     qr'          :: Matrix t -> (Matrix t, Matrix t)
     hess'        :: Matrix t -> (Matrix t, Matrix t)
     schur'       :: Matrix t -> (Matrix t, Matrix t)
@@ -123,6 +124,7 @@ instance Field Double where
     eigOnly = eigOnlyR
     eigOnlySH = eigOnlyS
     cholSH' = cholS
+    mbCholSH' = mbCholS
     qr' = unpackQR . qrR
     hess' = unpackHess hessR
     schur' = schurR
@@ -149,6 +151,7 @@ instance Field (Complex Double) where
     eigSH'' = eigH
     eigOnlySH = eigOnlyH
     cholSH' = cholH
+    mbCholSH' = mbCholH
     qr' = unpackQR . qrC
     hess' = unpackHess hessC
     schur' = schurC
@@ -263,11 +266,11 @@ eig = eig'
 eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
 eigenvalues = eigOnly
 
--- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric.
+-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
 eigSH'      :: Field t => Matrix t -> (Vector Double, Matrix t)
 eigSH' = eigSH''
 
--- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric.
+-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
 eigenvaluesSH' :: Field t => Matrix t -> Vector Double
 eigenvaluesSH' = eigOnlySH
 
@@ -328,8 +331,11 @@ ctrans = ctrans'
 multiply :: Field t => Matrix t -> Matrix t -> Matrix t
 multiply = multiply'
 
+-- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'.
+mbCholSH :: Field t => Matrix t -> Maybe (Matrix t)
+mbCholSH = mbCholSH'
 
--- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
+-- | Similar to 'chol', without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
 cholSH      :: Field t => Matrix t -> Matrix t
 cholSH = cholSH'
 
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 5d4eb0d..539ffb9 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -32,7 +32,7 @@ module Numeric.LinearAlgebra.LAPACK (
     -- * LU
     luR, luC,
     -- * Cholesky
-    cholS, cholH,
+    cholS, cholH, mbCholS, mbCholH,
     -- * QR
     qrR, qrC,
     -- * Hessenberg
@@ -392,19 +392,27 @@ linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b)
 foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM
 foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
 
+cholAux f st a = do
+    r <- createMatrix ColumnMajor n n
+    app2 f mat a mat r st
+    return r
+  where n = rows a
+
 -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
 cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
-cholH = cholAux zpotrf "cholH" . fmat
+cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat
 
 -- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
 cholS :: Matrix Double -> Matrix Double
-cholS = cholAux dpotrf "cholS" . fmat
+cholS =  unsafePerformIO . cholAux dpotrf "cholS" . fmat
 
-cholAux f st a = unsafePerformIO $ do
-    r <- createMatrix ColumnMajor n n
-    app2 f mat a mat r st
-    return r
-  where n = rows a
+-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version).
+mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
+mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat
+
+-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/  ('Maybe' version).
+mbCholS :: Matrix Double -> Maybe (Matrix Double)
+mbCholS =  unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat
 
 -----------------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h qr_l_R" dgeqr2 :: TMVM
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 60c60e6..ac26466 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -161,6 +161,14 @@ fittingTest = utest "levmar" (ok1 && ok2)
     ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d
     ok2 = pnorm PNorm2 (fromList (map fst sols) - fromList sol) < 1E-5
 
+-----------------------------------------------------
+
+mbCholTest = utest "mbCholTest" (ok1 && ok2) where
+    m1 = (2><2) [2,5,5,8 :: Double]
+    m2 = (2><2) [3,5,5,9 :: Complex Double]
+    ok1 = mbCholSH m1 == Nothing
+    ok2 = mbCholSH m2 == Just (chol m2)
+
 ---------------------------------------------------------------------
 
 randomTestGaussian = c :~1~: snd (meanCov dat) where
@@ -325,6 +333,7 @@ runTests n = do
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
         , odeTest
         , fittingTest
+        , mbCholTest
         ]
     return ()