LAPACK cholesky

5 files changed, 75 insertions, 3 deletions

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"
 errorCode 2003 = "bad file"
 errorCode 2004 = "singular"
 errorCode 2005 = "didn't converge"
+errorCode 2006 = "the input matrix is not positive definite"
 errorCode n    = "code "++show n
 
 {- | 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 (
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
+    cholS, cholH
 ) where
 
 import Data.Packed.Internal
@@ -279,3 +280,27 @@ linearSolveSVDC_l rcond  a b = unsafePerformIO $ do
         n = cols a
         nrhs = cols b
 
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h chol_l_H" zpotrf :: TCMCM
+
+-- | Wrapper for LAPACK's /zpotrf/,which computes the Cholesky factorization of a
+-- complex Hermitian positive definite matrix.
+cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
+cholH a = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor n n
+    zpotrf // mat fdat a // mat dat r // check "cholH" [fdat a]
+    return r
+  where n = rows a
+
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h chol_l_S" dpotrf :: TMM
+
+-- | Wrapper for LAPACK's /dpotrf/,which computes the Cholesky factorization of a
+-- real symmetric positive definite matrix.
+cholS :: Matrix Double -> Matrix Double
+cholS a = unsafePerformIO $ do
+    r <- createMatrix ColumnMajor n n
+    dpotrf // mat fdat a // mat dat r // check "cholS" [fdat a]
+    return r
+  where n = rows a
+
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 @@
 #define BAD_FILE 2003
 #define SINGULAR 2004
 #define NOCONVER 2005
+#define NODEFPOS 2006
 
 //////////////////// real svd ////////////////////////////////////
 
@@ -577,3 +578,44 @@ int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) {
     free(AC);
     OK
 }
+
+//////////////////// Cholesky factorization /////////////////////////
+
+int chol_l_H(KCMAT(a),CMAT(l)) {
+    integer n = ar;
+    REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
+    DEBUGMSG("chol_l_H");
+    memcpy(lp,ap,n*n*sizeof(doublecomplex));
+    char uplo = 'U';
+    integer res;
+    zpotrf_ (&uplo,&n,(doublecomplex*)lp,&n,&res);
+    CHECK(res>0,NODEFPOS);
+    CHECK(res,res);
+    doublecomplex zero = {0.,0.};
+    int r,c;
+    for (r=0; r<lr-1; r++) {
+        for(c=r+1; c<lc; c++) {
+            ((doublecomplex*)lp)[r*lc+c] = zero;
+        }
+    }
+    OK
+}
+
+int chol_l_S(KDMAT(a),DMAT(l)) {
+    integer n = ar;
+    REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE);
+    DEBUGMSG("chol_l_S");
+    memcpy(lp,ap,n*n*sizeof(double));
+    char uplo = 'U';
+    integer res;
+    dpotrf_ (&uplo,&n,lp,&n,&res);
+    CHECK(res>0,NODEFPOS);
+    CHECK(res,res);
+    int r,c;
+    for (r=0; r<lr-1; r++) {
+        for(c=r+1; c<lc; c++) {
+            lp[r*lc+c] = 0.;
+        }
+    }
+    OK
+}
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));
 int linearSolveSVDR_l(double,KDMAT(a),KDMAT(b),DMAT(x));
 
 int linearSolveSVDC_l(double,KCMAT(a),KCMAT(b),CMAT(x));
+
+int chol_l_H(KCMAT(a),CMAT(r));
+
+int 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 (
 
 import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
 import Data.Packed
-import GSL.Matrix
+import GSL.Matrix(luR,luC,qr)
 import GSL.Vector
 import LAPACK
 import Complex
@@ -70,7 +70,7 @@ instance GenMat Double where
     ctrans = trans
     eig = eigR
     eigSH = LAPACK.eigS
-    chol = cholR
+    chol = cholS
 
 instance GenMat (Complex Double) where
     svd = svdC
@@ -80,7 +80,7 @@ instance GenMat (Complex Double) where
     ctrans = conjTrans
     eig = eigC
     eigSH =  LAPACK.eigH
-    chol = error "cholC not yet implemented" -- waiting for GSL-1.10
+    chol = cholH
 
 -- | eigensystem of a symmetric matrix
 eigS :: Matrix Double -> (Vector Double, Matrix Double)