4 files changed, 120 insertions, 18 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 069d9a3..b19c0ec 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -37,6 +37,8 @@ module Numeric.LinearAlgebra.Algorithms (
     hess,
 -- ** Schur
     schur,
+-- ** LU
+    lu,
 -- * Matrix functions
     expm,
     sqrtm,
@@ -52,11 +54,11 @@ module Numeric.LinearAlgebra.Algorithms (
 -- * Util
     haussholder,
     unpackQR, unpackHess,
-    Field(linearSolveSVD,lu,eigSH',cholSH)
+    Field(linearSolveSVD,eigSH',cholSH)
 ) where
 
 
-import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj)
+import Data.Packed.Internal hiding (fromComplex, toComplex, comp, conj, (//))
 import Data.Packed
 import qualified Numeric.GSL.Matrix as GSL
 import Numeric.GSL.Vector
@@ -64,12 +66,13 @@ import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Complex
 import Numeric.LinearAlgebra.Linear
 import Data.List(foldl1')
+import Data.Array
 
 -- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
 class (Normed (Matrix t), Linear Matrix t) => Field t where
     -- | Singular value decomposition using lapack's dgesvd or zgesvd.
     svd         :: Matrix t -> (Matrix t, Vector Double, Matrix t)
-    lu          :: Matrix t -> (Matrix t, Matrix t, [Int], t)
+    luPacked    :: Matrix t -> (Matrix t, [Int])
     -- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv and zgesv.
     --  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
     linearSolve :: Matrix t -> Matrix t -> Matrix t
@@ -106,7 +109,7 @@ class (Normed (Matrix t), Linear Matrix t) => Field t where
 
 instance Field Double where
     svd = svdR
-    lu  = GSL.luR
+    luPacked = luR
     linearSolve = linearSolveR
     linearSolveSVD = linearSolveSVDR Nothing
     ctrans = trans
@@ -119,7 +122,7 @@ instance Field Double where
 
 instance Field (Complex Double) where
     svd = svdC
-    lu  = GSL.luC
+    luPacked = luC
     linearSolve = linearSolveC
     linearSolveSVD = linearSolveSVDC Nothing
     ctrans = conj . trans
@@ -146,10 +149,19 @@ chol m | m `equal` ctrans m = cholSH m
 
 square m = rows m == cols m
 
+-- | determinant of a square matrix, computed from the LU decomposition.
 det :: Field t => Matrix t -> t
-det m | square m = s * (product $ toList $ takeDiag $ u)
+det m | square m = s * (product $ toList $ takeDiag $ lu)
       | otherwise = error "det of nonsquare matrix"
-    where (_,u,_,s) = lu m
+    where (lu,perm) = luPacked m
+          s = signlp (rows m) perm
+
+-- | LU factorization of a general matrix using lapack's dgetrf or zgetrf.
+--
+-- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
+-- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
+lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
+lu = luFact . luPacked
 
 -- | Inverse of a square matrix using lapacks' dgesv and zgesv.
 inv :: Field t => Matrix t -> Matrix t
@@ -457,3 +469,35 @@ sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
           (.*) = scale
           (|+|) = add
           (|-|) = sub
+
+------------------------------------------------------------------
+
+signlp r vals = foldl f 1 (zip [0..r-1] vals)
+    where f s (a,b) | a /= b    = -s
+                    | otherwise =  s
+
+swap (arr,s) (a,b) | a /= b    = (arr // [(a, arr!b),(b,arr!a)],-s)
+                   | otherwise = (arr,s)
+
+fixPerm r vals = (fromColumns $ elems res, sign)
+    where v = [0..r-1]
+          s = toColumns (ident r)
+          (res,sign) = foldl swap (listArray (0,r-1) s, 1) (zip v vals)
+
+triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]
+    where el i j = if j-i>=h then v else 1 - v
+
+luFact (lu,perm) | r <= c    = (l ,u ,p, s)
+                 | otherwise = (l',u',p, s)
+  where
+    r = rows lu
+    c = cols lu
+    tu = triang r c 0 1
+    tl = triang r c 0 0
+    l = takeColumns r (lu |*| tl) |+| diagRect (constant 1 r) r r
+    u = lu |*| tu
+    (p,s) = fixPerm r perm
+    l' = (lu |*| tl) |+| diagRect (constant 1 c) r c
+    u' = takeRows c (lu |*| tu)
+    (|+|) = add
+    (|*|) = mul
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index cacad87..83db901 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -19,6 +19,7 @@ module Numeric.LinearAlgebra.LAPACK (
     linearSolveR, linearSolveC,
     linearSolveLSR, linearSolveLSC,
     linearSolveSVDR, linearSolveSVDC,
+    luR, luC,
     cholS, cholH,
     qrR, qrC,
     hessR, hessC,
@@ -299,7 +300,7 @@ hessAux f st a = unsafePerformIO $ do
         mn = min m n
 
 -----------------------------------------------------------------------------------
-foreign import ccall safe "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
+foreign import ccall "LAPACK/lapack-aux.h schur_l_R" dgees :: TMMM
 foreign import ccall "LAPACK/lapack-aux.h schur_l_C" zgees :: TCMCMCM
 
 -- | Wrapper for LAPACK's /dgees/, which computes a Schur factorization of a square real matrix.
@@ -318,3 +319,21 @@ schurAux f st a = unsafePerformIO $ do
   where n = rows a
 
 -----------------------------------------------------------------------------------
+foreign import ccall "LAPACK/lapack-aux.h lu_l_R" dgetrf :: TMVM
+foreign import ccall "LAPACK/lapack-aux.h lu_l_C" zgetrf :: TCMVCM
+
+-- | Wrapper for LAPACK's /dgetrf/, which computes a LU factorization of a general real matrix.
+luR :: Matrix Double -> (Matrix Double, [Int])
+luR = luAux dgetrf "luR" . fmat
+
+-- | Wrapper for LAPACK's /zgees/, which computes a Schur factorization of a square complex matrix.
+luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
+luC = luAux zgetrf "luC" . fmat
+
+luAux f st a = unsafePerformIO $ do
+    lu <- createMatrix ColumnMajor n m
+    piv <- createVector (min n m)
+    app3 f mat a vec piv mat lu st
+    return (lu, map (pred.round) (toList piv))
+  where n = rows a
+        m = cols a
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));