4 files changed, 41 insertions, 23 deletions

diff --git a/lib/GSL/Matrix.hs b/lib/GSL/Matrix.hs
index ec8ceea..c1bce37 100644
--- a/lib/GSL/Matrix.hs
+++ b/lib/GSL/Matrix.hs
@@ -16,7 +16,7 @@ module GSL.Matrix(
     eigSg, eigHg,
     svdg,
     qr,
-    chol,
+    cholR, -- cholC,
     luSolveR, luSolveC,
     luR, luC,
     fromFile, extractRows
@@ -153,24 +153,29 @@ foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
 
 {- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
 
-@\> let c = chol $ 'fromLists' [[5,4],[4,5]]
-\
-\> c
-2.236    0.
-1.789 1.342
-\
-\> c \<\> 'trans' c
-5.000 4.000
-4.000 5.000@
+@\> chol $ (2><2) [1,2,
+                   2,9::Double]
+(2><2)
+ [ 1.0,              0.0
+ , 2.0, 2.23606797749979 ]@
 
 -}
-chol :: Matrix Double -> Matrix Double
-chol x = unsafePerformIO $ do
+cholR :: Matrix Double -> Matrix Double
+cholR x = unsafePerformIO $ do
     res <- createMatrix RowMajor r r
-    c_chol // mat cdat x // mat dat res // check "chol" [cdat x]
+    c_cholR // mat cdat x // mat dat res // check "cholR" [cdat x]
     return res
   where r = rows x
-foreign import ccall "gsl-aux.h chol" c_chol :: TMM
+foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
+
+cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
+cholC x = unsafePerformIO $ do
+    res <- createMatrix RowMajor r r
+    c_cholC // mat cdat x // mat dat res // check "cholC" [cdat x]
+    return res
+  where r = rows x
+foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
+
 
 --------------------------------------------------------
 
diff --git a/lib/GSL/Special.hs b/lib/GSL/Special.hs
index bdee1b2..fa002b9 100644
--- a/lib/GSL/Special.hs
+++ b/lib/GSL/Special.hs
@@ -41,7 +41,6 @@ module GSL.Special (
 , module GSL.Special.Synchrotron
 , module GSL.Special.Trig
 , module GSL.Special.Zeta
-, setErrorHandlerOff
 )
 where
 
@@ -73,8 +72,3 @@ import GSL.Special.Psi
 import GSL.Special.Synchrotron
 import GSL.Special.Trig
 import GSL.Special.Zeta
-
-
--- | This action removes the GSL default error handler which aborts the program, so that
--- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
-foreign import ccall "gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
diff --git a/lib/GSL/gsl-aux.c b/lib/GSL/gsl-aux.c
index 0e8406c..c602d5e 100644
--- a/lib/GSL/gsl-aux.c
+++ b/lib/GSL/gsl-aux.c
@@ -471,9 +471,9 @@ int QR(KRMAT(x),RMAT(q),RMAT(r)) {
     OK
 }
 
-int chol(KRMAT(x),RMAT(l)) {
+int cholR(KRMAT(x),RMAT(l)) {
     REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
-    DEBUGMSG("chol");
+    DEBUGMSG("cholR");
     KDMVIEW(x);
     DMVIEW(l);
     gsl_matrix_memcpy(M(l),M(x));
@@ -488,6 +488,23 @@ int chol(KRMAT(x),RMAT(l)) {
     OK
 }
 
+int cholC(KCMAT(x),CMAT(l)) {
+    REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
+    DEBUGMSG("cholC");
+    KCMVIEW(x);
+    CMVIEW(l);
+    gsl_matrix_complex_memcpy(M(l),M(x));
+    int res = 0; // gsl_linalg_complex_cholesky_decomp(M(l));
+    CHECK(res,res);
+    gsl_complex zero = {0.,0.};
+    int r,c;
+    for (r=0; r<xr-1; r++) {
+        for(c=r+1; c<xc; c++) {
+            lp[r*lc+c] = zero;
+        }
+    }
+    OK
+}
 
 int fft(int code, KCVEC(X), CVEC(R)) {
     REQUIRES(Xn == Rn,BAD_SIZE);
diff --git a/lib/GSL/gsl-aux.h b/lib/GSL/gsl-aux.h
index 8f8618a..dd1a247 100644
--- a/lib/GSL/gsl-aux.h
+++ b/lib/GSL/gsl-aux.h
@@ -38,7 +38,9 @@ int eigensystemC(KCMAT(x),RVEC(l),CMAT(v));
 
 int QR(KRMAT(x),RMAT(q),RMAT(r));
 
-int chol(KRMAT(x),RMAT(l));
+int cholR(KRMAT(x),RMAT(l));
+
+int cholC(KCMAT(x),CMAT(l));
 
 int fft(int code, KCVEC(a), CVEC(b));