7 files changed, 46 insertions, 15 deletions

diff --git a/CHANGES b/CHANGES
index f4d634e..69902ef 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,7 +1,7 @@
 0.9.2.0
 =======
 
-- cholSolve
+- cholSolve, mbCholSH
 
 - GSL Nonlinear Least-Squares fitting using Levenberg-Marquardt.
 
diff --git a/examples/error.hs b/examples/error.hs
index f13a48e..5efae7c 100644
--- a/examples/error.hs
+++ b/examples/error.hs
@@ -1,4 +1,5 @@
 import Numeric.GSL
+import Numeric.GSL.Special
 import Numeric.LinearAlgebra
 import Prelude hiding (catch)
 import Control.Exception
diff --git a/examples/plot.hs b/examples/plot.hs
index 1583dfe..f950aa5 100644
--- a/examples/plot.hs
+++ b/examples/plot.hs
@@ -1,6 +1,6 @@
 import Numeric.LinearAlgebra
 import Graphics.Plot
-import Numeric.GSL(erf_Z, erf)
+import Numeric.GSL.Special(erf_Z, erf)
 
 sombrero n = f x y where 
     (x,y) = meshdom range range
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index bfa63f1..455b176 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -17,7 +17,7 @@
 module Data.Packed.Internal.Common(
   Adapt,
   app1, app2, app3, app4,
-  (//), check,
+  (//), check, mbCatch,
   splitEvery, common, compatdim,
   fi,
   table
@@ -29,6 +29,7 @@ import Foreign.C.String(peekCString)
 import Foreign.C.Types
 import Foreign.Storable.Complex()
 import Data.List(transpose,intersperse)
+import Control.Exception as E
 
 -- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
 splitEvery :: Int -> [a] -> [[a]]
@@ -151,3 +152,9 @@ check msg f = do
 
 -- | description of GSL error codes
 foreign import ccall "auxi.h gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
+
+-- | Error capture and conversion to Maybe
+mbCatch :: IO x -> IO (Maybe x)
+mbCatch act = E.catch (Just `fmap` act) f
+    where f :: SomeException -> IO (Maybe x)
+          f _ = return Nothing
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 ()