tcov and ldl python

SPEX/Config/Tcov_Makefile.in

@@ -199,6 +199,7 @@ OBJ_Cholesky = \
     SPEX_ldl_analyze.o \
     SPEX_ldl_factorize.o \
     SPEX_ldl_solve.o \
+    SPEX_ldl_backslash.o \
 
 $(OBJ_Cholesky): ../Include/SPEX.h ../SPEX_Cholesky/Source/spex_cholesky_internal.h
 

SPEX/Python/SPEXpy/Source/spex_python_connect.c

@@ -121,6 +121,9 @@ SPEX_info spex_python
         case 3:
             SPEX_CHECK( SPEX_cholesky_backslash(&x, SPEX_MPQ, A, b, option));
             break;
+        case 4:
+            SPEX_CHECK( SPEX_ldl_backslash(&x, SPEX_MPQ, A, b, option));
+            break;
         default:
             return SPEX_INCORRECT_INPUT;
     }

SPEX/Python/SPEXpy/ldl_backslash.py

@@ -0,0 +1,49 @@
+#-------------------------------------------------------------------------------
+# SPEX/Python/SPEXpy/cholesky_backslash.py: solve Ax=b using Cholesky factorization
+#-------------------------------------------------------------------------------
+
+# SPEX: (c) 2022-2024, Christopher Lourenco, Jinhao Chen,
+# Lorena Mejia Domenzain, Erick Moreno-Centeno, and Timothy A. Davis.
+# All Rights Reserved.
+# SPDX-License-Identifier: GPL-2.0-or-later or LGPL-3.0-or-later
+
+#------------------------------------------------------------------------------
+
+from .Options import Options
+from .SPEX_error import *
+from .spex_connect import spex_connect
+
+import scipy
+from scipy.sparse import isspmatrix, isspmatrix_csc, linalg
+
+def cholesky_backslash( A, b, options=Options('double', 'amd')):
+    ## A is a scipy.sparse(data must be float64) #technically it only needs to be numerical
+    ## b is a numpy.array (data must be float64)
+    ## options is a dictionary that specifies what tipe the solution should be, this by default is double
+
+    ##--------------------------------------------------------------------------
+    ## Verify inputs
+    ##--------------------------------------------------------------------------
+    if not isspmatrix(A):
+        raise SPEXerror(determine_error(3))
+    ## If the sparse input matrix is not in csc form, convert it into csc form
+    if not isspmatrix_csc(A):
+        A.tocsc()
+    ## Check symmetry
+    tol=1e-8
+    if linalg.norm(A-A.T, scipy.Inf) > tol:
+        raise SPEX_error(determine_error(-4))
+    # Check input shape
+    if A.shape[1]!=b.shape[0]:
+        raise SPEX_error(determine_error(-3))
+
+    if options.ordering==None:
+        options.default_chol()
+
+    ##--------------------------------------------------------------------------
+    ## Call SPEX
+    ##--------------------------------------------------------------------------
+    x=spex_connect(A,b,options.order(),options.charOut(),4) #4 calls the LDL factorization
+
+    return x
+

SPEX/SPEX_Cholesky/Source/SPEX_cholesky_solve.c

@@ -47,7 +47,23 @@ SPEX_info SPEX_cholesky_solve
 {
     // Just need to call the symmetric solve with chol = true
     SPEX_info info;
+        // Ensure SPEX is initialized
+    if (!spex_initialized())
+    {
+        return SPEX_PANIC;
+    }
+
+    // Check the inputs
+    if (!x_handle || b->type != SPEX_MPZ || b->kind != SPEX_DENSE)
+    {
+        return SPEX_INCORRECT_INPUT;
+    }
 
+    if (F->kind != SPEX_CHOLESKY_FACTORIZATION)
+    {
+        return SPEX_INCORRECT_INPUT;
+    }
+
     info = spex_symmetric_solve(x_handle, F, b, true, option);
 
     return info;

SPEX/SPEX_Cholesky/Source/SPEX_ldl_solve.c

@@ -48,6 +48,23 @@ SPEX_info SPEX_ldl_solve
     // Just need to call the symmetric solve with chol = false
     SPEX_info info;
 
+        // Ensure SPEX is initialized
+    if (!spex_initialized())
+    {
+        return SPEX_PANIC;
+    }
+
+    // Check the inputs
+    if (!x_handle || b->type != SPEX_MPZ || b->kind != SPEX_DENSE)
+    {
+        return SPEX_INCORRECT_INPUT;
+    }
+
+    if (F->kind != SPEX_LDL_FACTORIZATION)
+    {
+        return SPEX_INCORRECT_INPUT;
+    }
+
     info = spex_symmetric_solve(x_handle, F, b, false, option);
 
     return info;

SPEX/SPEX_Cholesky/Source/spex_symmetric_solve.c

@@ -62,27 +62,7 @@ SPEX_info spex_symmetric_solve
 
     SPEX_info info;
 
-    // Ensure SPEX is initialized
-    if (!spex_initialized())
-    {
-        return SPEX_PANIC;
-    }
-
-    // Check the inputs
-    if (!x_handle || b->type != SPEX_MPZ || b->kind != SPEX_DENSE)
-    {
-        return SPEX_INCORRECT_INPUT;
-    }
-
-    if (chol && F->kind != SPEX_CHOLESKY_FACTORIZATION)
-    {
-        return SPEX_INCORRECT_INPUT;
-    }
-
-    if (!chol && F->kind != SPEX_LDL_FACTORIZATION)
-    {
-        return SPEX_INCORRECT_INPUT;
-    }
+    SPEX_REQUIRE(b, SPEX_DENSE, SPEX_MPZ);
 
     // det is the determinant of the PAP matrix. It is obtained for free
     // from the SPEX Cholesky factorization det = rhos[n-1] = L[n,n]

SPEX/Tcov/tcov_for_cholesky.c

@@ -113,6 +113,24 @@ SPEX_info spex_test_chol_backslash (SPEX_matrix A, SPEX_matrix b,
     return (SPEX_OK) ;
 }
 
+SPEX_info spex_test_ldl_backslash (SPEX_matrix A, SPEX_matrix b,
+    SPEX_options option);
+
+SPEX_info spex_test_ldl_backslash (SPEX_matrix A, SPEX_matrix b,
+    SPEX_options option)
+{
+    SPEX_matrix x = NULL ;
+    // solve Ax=b
+    OK2 (SPEX_ldl_backslash (&x, SPEX_MPQ, A, b, option));
+    // disable memory testing when checking the solution
+    int64_t save = malloc_count ; malloc_count = INT64_MAX ;
+    OK (spex_demo_check_solution (A, x, b, option));
+    // re-enable memory testing
+    malloc_count = save ;
+    SPEX_FREE_ALL;
+    return (SPEX_OK) ;
+}
+
 //------------------------------------------------------------------------------
 // spex_test_cdiv_qr: test SPEX_cdiv_qr
 //------------------------------------------------------------------------------
@@ -181,6 +199,31 @@ SPEX_info spex_test_chol_afs (SPEX_matrix A, SPEX_matrix b, SPEX_options option)
     return (SPEX_OK);
 }
 
+SPEX_info spex_test_ldl_afs
+(
+    SPEX_matrix A,
+    SPEX_matrix b,
+    SPEX_options option
+) ;
+
+SPEX_info spex_test_ldl_afs (SPEX_matrix A, SPEX_matrix b, SPEX_options option)
+{
+    SPEX_symbolic_analysis S = NULL ;
+    SPEX_factorization F = NULL ;
+    SPEX_matrix x = NULL ;
+    // solve Ax=b
+    OK2 (SPEX_ldl_analyze (&S, A, option));
+    OK2 (SPEX_ldl_factorize (&F, A, S, option));
+    OK2 (SPEX_ldl_solve (&x, F, b, option));
+    // disable memory testing when checking the solution
+    int64_t save = malloc_count ; malloc_count = INT64_MAX ;
+    OK (spex_demo_check_solution (A, x, b, option));
+    // re-enable memory testing
+    malloc_count = save ;
+    SPEX_FREE_ALL;
+    return (SPEX_OK);
+}
+
 //------------------------------------------------------------------------------
 // tcov_for_cholesky: main program
 //------------------------------------------------------------------------------
@@ -232,23 +275,23 @@ int main (int argc, char *argv [])
     printf ("Cholesky: error handling for unsymmetric matrix (1)\n");
     read_test_matrix (&A, "../ExampleMats/test1.mat.txt");
     create_test_rhs (&b, A->n);
-    ERR (SPEX_cholesky_backslash (&x, SPEX_MPQ, A, b, option), SPEX_NOTSPD);
+    ERR (SPEX_cholesky_backslash (&x, SPEX_MPQ, A, b, option), SPEX_UNSYMMETRIC);
     OK (SPEX_matrix_free (&A, option));
     OK (SPEX_matrix_free (&b, option));
 
     // unsymmetric matrix (unsymmetric pattern)
     printf ("Cholesky: error handling for unsymmetric matrix (2)\n");
     read_test_matrix (&A, "../ExampleMats/test2.mat.txt");
     create_test_rhs (&b, A->n);
-    ERR (SPEX_cholesky_backslash (&x, SPEX_MPQ, A, b, option), SPEX_NOTSPD);
+    ERR (SPEX_cholesky_backslash (&x, SPEX_MPQ, A, b, option), SPEX_UNSYMMETRIC);
     OK (SPEX_matrix_free (&A, option));
     OK (SPEX_matrix_free (&b, option));
 
     // unsymmetric matrix (unsymmetric values)
     printf ("Cholesky: error handling for unsymmetric matrix (3)\n");
     read_test_matrix (&A, "../ExampleMats/test3.mat.txt");
     create_test_rhs (&b, A->n);
-    ERR (SPEX_cholesky_backslash (&x, SPEX_MPQ, A, b, option), SPEX_NOTSPD);
+    ERR (SPEX_cholesky_backslash (&x, SPEX_MPQ, A, b, option), SPEX_UNSYMMETRIC);
     OK (SPEX_matrix_free (&A, option));
     OK (SPEX_matrix_free (&b, option));
 
@@ -264,6 +307,17 @@ int main (int argc, char *argv [])
     OK (SPEX_matrix_free (&A, option));
     OK (SPEX_matrix_free (&b, option));
 
+    // symmetric indefinite matrix
+    printf ("Cholesky: error handling for symmetric matrix with zero in diagonal\n");
+    read_test_matrix (&A, "../ExampleMats/test4.mat.txt");
+    create_test_rhs (&b, A->n);
+    option->algo = SPEX_CHOL_UP ;
+    ERR (SPEX_ldl_backslash (&x, SPEX_MPQ, A, b, option), SPEX_ZERODIAG);
+    option->algo = SPEX_CHOL_LEFT ;
+    ERR (SPEX_ldl_backslash (&x, SPEX_MPQ, A, b, option), SPEX_ZERODIAG);
+    OK (SPEX_matrix_free (&A, option));
+    OK (SPEX_matrix_free (&b, option));
+
     //--------------------------------------------------------------------------
     // symetry check
     //--------------------------------------------------------------------------
@@ -479,6 +533,10 @@ int main (int argc, char *argv [])
     ERR (SPEX_cholesky_solve (&x, F, b, option),
         SPEX_INCORRECT_INPUT);
     b->type = SPEX_MPZ ;
+    F->kind = SPEX_LDL_FACTORIZATION;
+    ERR (SPEX_cholesky_solve (&x, F, b, option),
+        SPEX_INCORRECT_INPUT);
+    F->kind = SPEX_CHOLESKY_FACTORIZATION;
     OK (SPEX_symbolic_analysis_free (&S, option));
     OK (SPEX_factorization_free (&F, option));
 
@@ -488,6 +546,33 @@ int main (int argc, char *argv [])
         SPEX_INCORRECT_ALGORITHM);
     option->algo = SPEX_CHOL_UP ;
 
+
+    ERR (SPEX_ldl_factorize (NULL, NULL, NULL, NULL),
+        SPEX_INCORRECT_INPUT);
+
+    // valid analysis, but break the factorization
+    OK (SPEX_ldl_analyze (&S, A, option));
+    A->type = SPEX_INT64 ;
+    ERR (SPEX_ldl_factorize (&F, A, S, option),
+        SPEX_INCORRECT_INPUT);
+    A->type = SPEX_MPZ ;
+    OK (SPEX_symbolic_analysis_free (&S, option));
+
+    // valid analysis and factorization, but break the solve
+    OK (SPEX_ldl_analyze (&S, A, option));
+    OK (SPEX_ldl_factorize (&F, A, S, option));
+    b->type = SPEX_INT64 ;
+    ERR (SPEX_ldl_solve (&x, F, b, option),
+        SPEX_INCORRECT_INPUT);
+    b->type = SPEX_MPZ ;
+    F->kind = SPEX_CHOLESKY_FACTORIZATION;
+    ERR (SPEX_ldl_solve (&x, F, b, option),
+        SPEX_INCORRECT_INPUT);
+    F->kind = SPEX_LDL_FACTORIZATION;
+    OK (SPEX_symbolic_analysis_free (&S, option));
+    OK (SPEX_factorization_free (&F, option));
+
+
     //--------------------------------------------------------------------------
     // solve Ax=b with SPEX_cholesky_backslash and check the solution
     //--------------------------------------------------------------------------
@@ -549,6 +634,33 @@ int main (int argc, char *argv [])
     OK (SPEX_mpz_set_ui (b->x.mpz [0], 0));
     BRUTAL (spex_test_chol_afs (A, b, option));
 
+
+    //--------------------------------------------------------------------------
+    // solve Ax=b with SPEX_ldl_backslash and check the solution
+    //--------------------------------------------------------------------------
+
+    option->order = SPEX_AMD ;
+    option->algo = SPEX_CHOL_UP ;
+    option->print_level = 3 ;
+    printf ("LDL backslash, up-looking, no malloc testing:\n");
+    OK (spex_test_ldl_backslash (A, b, option));
+    option->print_level = 0 ;
+
+    printf ("LDL backslash, left-looking with malloc testing:\n");
+    option->algo = SPEX_CHOL_LEFT ;
+    BRUTAL (spex_test_ldl_backslash (A, b, option));
+
+    //--------------------------------------------------------------------------
+    // solve Ax=b with SPEX_cholesky_[analyze,factorize,solve]; check solution
+    //--------------------------------------------------------------------------
+
+    option->algo = SPEX_CHOL_UP ;
+
+    printf ("LDL analyze/factorize/solve, no malloc testing:\n");
+    spex_set_gmp_ntrials (INT64_MAX) ;
+    malloc_count = INT64_MAX ;
+    OK (spex_test_ldl_afs (A, b, option));
+
     //--------------------------------------------------------------------------
     // error handling
     //--------------------------------------------------------------------------
@@ -560,6 +672,10 @@ int main (int argc, char *argv [])
     ERR (SPEX_cholesky_solve (NULL, NULL, NULL, NULL), SPEX_PANIC);
     ERR (SPEX_cholesky_backslash (NULL, SPEX_MPQ, NULL, NULL, NULL),
         SPEX_PANIC);
+
+    ERR (SPEX_ldl_factorize (&F2, A, S, option), SPEX_PANIC);
+    ERR (SPEX_ldl_solve (NULL, NULL, NULL, NULL), SPEX_PANIC);
+
     spex_set_initialized (true);
     SPEX_FREE_ALL;
     SPEX_finalize ( ) ;