2 files changed, 49 insertions, 10 deletions

diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 8392feb..310f6ee 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -768,3 +768,49 @@ int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) {
     OK
     #endif
 }
+
+//////////////////// LU factorization /////////////////////////
+
+int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
+    DEBUGMSG("lu_l_R");
+    integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
+    memcpy(rp,ap,m*n*sizeof(double));
+    integer res;
+    dgetrf_ (&m,&n,rp,&m,auxipiv,&res);
+    if(res>0) {
+        res = 0; // fixme
+    }
+    CHECK(res,res);
+    int k;
+    for (k=0; k<mn; k++) {
+        ipivp[k] = auxipiv[k];
+    }
+    free(auxipiv);
+    OK
+}
+
+int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
+    integer m = ar;
+    integer n = ac;
+    integer mn = MIN(m,n);
+    REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE);
+    DEBUGMSG("lu_l_C");
+    integer* auxipiv = (integer*)malloc(mn*sizeof(integer));
+    memcpy(rp,ap,m*n*sizeof(doublecomplex));
+    integer res;
+    zgetrf_ (&m,&n,(doublecomplex*)rp,&m,auxipiv,&res);
+    if(res>0) {
+        res = 0; // fixme
+    }
+    CHECK(res,res);
+    int k;
+    for (k=0; k<mn; k++) {
+        ipivp[k] = auxipiv[k];
+    }
+    free(auxipiv);
+    OK
+}
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index e5d74d7..bccd4b8 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -14,41 +14,34 @@
 
 int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
 int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v));
-
 int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));
 
 int eig_l_C(KCMAT(a),CMAT(u),CVEC(s),CMAT(v));
-
 int eig_l_R(KDMAT(a),DMAT(u),CVEC(s),DMAT(v));
 
 int eig_l_S(KDMAT(a),DVEC(s),DMAT(v));
-
 int eig_l_H(KCMAT(a),DVEC(s),CMAT(v));
 
 int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x));
-
 int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x));
 
 int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x));
-
 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));
 
 int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r));
-
 int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r));
 
 int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r));
-
 int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r));
 
 int schur_l_R(KDMAT(a), DMAT(u), DMAT(s));
-
 int schur_l_C(KCMAT(a), CMAT(u), CMAT(s));
+
+int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r));
+int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r));