LAPACK cholesky

2 files changed, 46 insertions, 0 deletions

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));