ldl factorization

6 files changed, 181 insertions, 8 deletions

diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index 6ce1830..c8b2d3e 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -55,6 +55,8 @@ class (Product t,
     mbLinearSolve' :: Matrix t -> Matrix t -> Maybe (Matrix t)
     linearSolve' :: Matrix t -> Matrix t -> Matrix t
     cholSolve'   :: Matrix t -> Matrix t -> Matrix t
+    ldlPacked'   :: Matrix t -> (Matrix t, [Int])
+    ldlSolve'    :: (Matrix t, [Int]) -> Matrix t -> Matrix t
     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
     linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t
     eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
@@ -90,6 +92,8 @@ instance Field Double where
     qrgr' = qrgrR
     hess' = unpackHess hessR
     schur' = schurR
+    ldlPacked' = ldlR
+    ldlSolve'= uncurry ldlsR
 
 instance Field (Complex Double) where
 #ifdef NOZGESDD
@@ -117,6 +121,8 @@ instance Field (Complex Double) where
     qrgr' = qrgrC
     hess' = unpackHess hessC
     schur' = schurC
+    ldlPacked' = ldlC
+    ldlSolve' = uncurry ldlsC
 
 --------------------------------------------------------------
 
@@ -333,6 +339,23 @@ linearSolveSVD = {-# SCC "linearSolveSVD" #-} linearSolveSVD'
 linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolveLS = {-# SCC "linearSolveLS" #-} linearSolveLS'
 
+--------------------------------------------------------------------------------
+
+-- | Similar to 'ldlPacked', without checking that the input matrix is hermitian or symmetric. It works with the lower triangular part.
+ldlPackedSH :: Field t => Matrix t -> (Matrix t, [Int])
+ldlPackedSH = {-# SCC "ldlPacked" #-} ldlPacked'
+
+-- | Obtains the LDL decomposition of a matrix in a compact data structure suitable for 'ldlSolve'.
+ldlPacked :: Field t => Matrix t -> (Matrix t, [Int])
+ldlPacked m
+    | exactHermitian m = {-# SCC "ldlPacked" #-} ldlPackedSH m
+    | otherwise = error "ldlPacked requires complex Hermitian or real symmetrix matrix"
+
+
+-- | Solution of a linear system (for several right hand sides) from the precomputed LDL factorization obtained by 'ldlPacked'.
+ldlSolve :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
+ldlSolve = {-# SCC "ldlSolve" #-} ldlSolve'
+
 --------------------------------------------------------------
 
 {- | Eigenvalues (not ordered) and eigenvectors (as columns) of a general square matrix.
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 30689bf..177d373 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1086,6 +1086,113 @@ int luS_l_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {
     OK
 }
 
+
+//////////////////// LDL factorization /////////////////////////
+
+int dsytrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *ipiv,
+            doublereal *work, integer *lwork, integer *info);
+
+int ldl_R(DVEC(ipiv), ODMAT(r)) {
+    integer n = rr;
+    REQUIRES(n>=1 && rc==n && ipivn == n, BAD_SIZE);
+    DEBUGMSG("ldl_R");
+    integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+    integer res;
+    integer lda = rXc;
+    integer lwork = -1;
+    doublereal ans;
+    dsytrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);
+    lwork = ceil(ans);
+    doublereal* work = (doublereal*)malloc(lwork*sizeof(doublereal));
+    dsytrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);
+    CHECK(res,res);
+    int k;
+    for (k=0; k<n; k++) {
+        ipivp[k] = auxipiv[k];
+    }
+    free(auxipiv);
+    free(work);
+    OK
+}
+
+
+int zhetrf_(char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv,
+            doublecomplex *work, integer *lwork, integer *info);
+
+int ldl_C(DVEC(ipiv), OCMAT(r)) {
+    integer n = rr;
+    REQUIRES(n>=1 && rc==n && ipivn == n, BAD_SIZE);
+    DEBUGMSG("ldl_R");
+    integer* auxipiv = (integer*)malloc(n*sizeof(integer));
+    integer res;
+    integer lda = rXc;
+    integer lwork = -1;
+    doublecomplex ans;
+    zhetrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);
+    lwork = ceil(ans.r);
+    doublecomplex* work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    zhetrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);
+    CHECK(res,res);
+    int k;
+    for (k=0; k<n; k++) {
+        ipivp[k] = auxipiv[k];
+    }
+    free(auxipiv);
+    free(work);
+    OK
+
+}
+
+//////////////////// LDL solve /////////////////////////
+
+int dsytrs_(char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda,
+            integer *ipiv, doublereal *b, integer *ldb, integer *info);
+
+int ldl_S_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) {
+  integer m = ar;
+  integer n = ac;
+  integer lda = aXc;
+  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;
+  dsytrs_ ("L",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);
+  CHECK(res,res);
+  free(auxipiv);
+  OK
+}
+
+
+int zhetrs_(char *uplo, integer *n, integer *nrhs, doublecomplex *a, integer *lda,
+            integer *ipiv, doublecomplex *b, integer *ldb, integer *info);
+
+int ldl_S_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {
+    integer m = ar;
+    integer n = ac;
+    integer lda = aXc;
+    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;
+    zhetrs_ ("L",&n,&nrhs,(doublecomplex*)ap,&lda,auxipiv,bp,&mrhs,&res);
+    CHECK(res,res);
+    free(auxipiv);
+    OK
+}
+
+
 //////////////////// Matrix Product /////////////////////////
 
 void dgemm_(char *, char *, integer *, integer *, integer *,
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index f2fc68d..c2c140b 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -591,7 +591,7 @@ foreign import ccall unsafe "luS_l_C" zgetrs :: C ::> R :> C ::> Ok
 lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
 lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv b
 
--- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.
+-- | Solve a complex linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.
 lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
 lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv b
 
@@ -600,10 +600,41 @@ lusAux f st a piv b
          x <- copy ColumnMajor b
          f # a # piv' # x #| st
          return x
-    | otherwise = error $ st ++ " on LU factorization of nonsquare matrix"
+    | otherwise = error st
   where
     n1 = rows a
     n2 = cols a
     n = rows b
     piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double
 
+-----------------------------------------------------------------------------------
+foreign import ccall unsafe "ldl_R" dsytrf :: R :> R ::> Ok
+foreign import ccall unsafe "ldl_C" zhetrf :: R :> C ::> Ok
+
+-- | LDL factorization of a symmetric real matrix, using LAPACK's /dsytrf/.
+ldlR :: Matrix Double -> (Matrix Double, [Int])
+ldlR = ldlAux dsytrf "ldlR"
+
+-- | LDL factorization of a hermitian complex matrix, using LAPACK's /zhetrf/.
+ldlC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
+ldlC = ldlAux zhetrf "ldlC"
+
+ldlAux f st a = unsafePerformIO $ do
+    ldl <- copy ColumnMajor a
+    piv <- createVector (rows a)
+    f # piv # ldl #| st
+    return (ldl, map (pred.round) (toList piv))
+
+-----------------------------------------------------------------------------------
+
+foreign import ccall unsafe "ldl_S_R" dsytrs :: R ::> R :> R ::> Ok
+foreign import ccall unsafe "ldl_S_C" zsytrs :: C ::> R :> C ::> Ok
+
+-- | Solve a real linear system from a precomputed LDL decomposition ('ldlR'), using LAPACK's /dsytrs/.
+ldlsR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
+ldlsR a piv b = lusAux dsytrs "ldlsR" (fmat a) piv b
+
+-- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/.
+ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
+ldlsC a piv b = lusAux zsytrs "ldlsC" (fmat a) piv b
+
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 0b8abbb..dd4cc67 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -15,7 +15,7 @@ module Numeric.LinearAlgebra (
     -- * Basic types and data processing
     module Numeric.LinearAlgebra.Data,
 
-    -- * Arithmetic and numeric classes
+    -- * Numeric classes
     -- |
     -- The standard numeric classes are defined elementwise:
     --
@@ -27,7 +27,9 @@ module Numeric.LinearAlgebra (
     --  [ 1.0, 0.0, 0.0
     --  , 0.0, 5.0, 0.0
     --  , 0.0, 0.0, 9.0 ]
-    --
+
+    -- * Autoconformable dimensions
+    -- |
     -- In arithmetic operations single-element vectors and matrices
     -- (created from numeric literals or using 'scalar') automatically
     -- expand to match the dimensions of the other operand:
@@ -79,6 +81,7 @@ module Numeric.LinearAlgebra (
     luSolve,
     luSolve',
     cholSolve,
+    ldlSolve,
     cgSolve,
     cgSolve',
 
@@ -115,15 +118,18 @@ module Numeric.LinearAlgebra (
     -- * Cholesky
     chol, cholSH, mbCholSH,
 
+    -- * LU
+    lu, luPacked, luPacked', luFact,
+
+    -- * LDL
+    ldlPacked, ldlPackedSH,
+
     -- * Hessenberg
     hess,
 
     -- * Schur
     schur,
 
-    -- * LU
-    lu, luPacked, luPacked', luFact,
-
     -- * Matrix functions
     expm,
     sqrtm,
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index d9bc9a0..2ff1580 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -587,6 +587,9 @@ runTests n = do
     putStrLn "------ luSolve"
     test (linearSolveProp (luSolve.luPacked) . rSqWC)
     test (linearSolveProp (luSolve.luPacked) . cSqWC)
+    putStrLn "------ ldlSolve"
+    test (linearSolveProp (ldlSolve.ldlPacked) . rSymWC)
+    test (linearSolveProp (ldlSolve.ldlPacked) . cSymWC)
     putStrLn "------ cholSolve"
     test (linearSolveProp (cholSolve.chol) . rPosDef)
     test (linearSolveProp (cholSolve.chol) . cPosDef)
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 904ae05..7c54535 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -16,7 +16,7 @@ module Numeric.LinearAlgebra.Tests.Instances(
     Rot(..),    rRot,cRot,
     Her(..),    rHer,cHer,
     WC(..),     rWC,cWC,
-    SqWC(..),   rSqWC, cSqWC,
+    SqWC(..),   rSqWC, cSqWC, rSymWC, cSymWC,
     PosDef(..), rPosDef, cPosDef,
     Consistent(..), rConsist, cConsist,
     RM,CM, rM,cM,
@@ -176,6 +176,9 @@ cWC (WC m) = m :: CM
 rSqWC (SqWC m) = m :: RM
 cSqWC (SqWC m) = m :: CM
 
+rSymWC (SqWC m) = m + tr m :: RM
+cSymWC (SqWC m) = m + tr m :: CM
+
 rPosDef (PosDef m) = m :: RM
 cPosDef (PosDef m) = m :: CM