4 files changed, 49 insertions, 1 deletions

diff --git a/README b/README
index efb32da..c034e2c 100644
--- a/README
+++ b/README
@@ -217,3 +217,5 @@ in the Haskell mailing lists for their help.
 
 - Pedro E. López de Teruel discovered the need of asm("finit") to
   avoid the wrong NaNs produced by foreign functions.
+
+- Reiner Pope added support for luSolve, based on (d|z)getrs.
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 83db901..d78b506 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -19,7 +19,7 @@ module Numeric.LinearAlgebra.LAPACK (
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
-    luR, luC,
+    luR, luC, lusR,
     cholS, cholH,
     qrR, qrC,
     hessR, hessC,
@@ -337,3 +337,18 @@ luAux f st a = unsafePerformIO $ do
     return (lu, map (pred.round) (toList piv))
   where n = rows a
         m = cols a
+
+
+-----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h luS_l_R" dgetrs :: TMVMM
+
+lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
+lusR a piv b = lusR' (fmat a) piv (fmat b)
+
+lusR' a piv b = unsafePerformIO $ do
+         x <- createMatrix ColumnMajor n m
+         app4 dgetrs mat a vec piv' mat b mat x "lusR"
+         return x
+    where n = rows b
+          m = cols b
+          piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 310f6ee..de3cc98 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -20,6 +20,9 @@
 #define OK return 0;
 #endif
 
+#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \
+                     for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");}
+
 #define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
 
 #define BAD_SIZE 2000
@@ -814,3 +817,29 @@ int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
     free(auxipiv);
     OK
 }
+
+
+//////////////////// LU substitution /////////////////////////
+char charN = 'N';
+
+int luS_l_R(DMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) {
+  integer m = ar;
+  integer n = ac;
+  integer mrhs = br;
+  integer nrhs = bc;
+
+  REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
+  integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+  int k;
+  for (k=0; k<n; k++) {
+    auxipiv[k] = (integer)ipivp[k];
+  }
+  integer res;
+  memcpy(xp,bp,mrhs*nrhs*sizeof(double));
+  integer ldb = mrhs;
+  integer lda = n;
+  dgetrs_ (&charN,&n,&nrhs,ap,&lda,auxipiv,xp,&ldb,&res);
+  CHECK(res,res);
+  free(auxipiv);
+  OK
+}
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index 79e52be..e98e9dc 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -84,3 +84,5 @@ 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));
+
+int luS_l_R(DMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x));