Add QDLDL_symbolic and QDLDL_numeric for separate factorization phases

PTNobel · claude · PTNobel · commit 8c47410613b6 · 2026-03-31T23:28:00.000-07:00
Split QDLDL_factor into QDLDL_symbolic (computes Lp/Li sparsity pattern
without numeric values) and QDLDL_numeric (computes Lx/D/Dinv given
pre-computed pattern). QDLDL_factor now calls both internally, preserving
backward compatibility. This enables symbolic-only use and efficient
refactorization when the sparsity pattern is unchanged.

Co-Authored-By: Claude Opus 4.6 (1M context) &lt;noreply@anthropic.com&gt;
diff --git a/include/qdldl.h b/include/qdldl.h
@@ -94,6 +94,79 @@ QDLDL_API QDLDL_int QDLDL_etree(const QDLDL_int n, const QDLDL_int* Ap, const QD
                                 QDLDL_int* work, QDLDL_int* Lnz, QDLDL_int* etree);
 
 
+/**
+ * Compute the symbolic factorization (sparsity pattern) of the L factor
+ * for a quasidefinite matrix in compressed sparse column form, where
+ * the input matrix is assumed to contain data for the upper triangular
+ * part of A only, and there are no duplicate indices.
+ *
+ * Returns the column pointers Lp and row indices Li of L, which together
+ * define the sparsity pattern of L in CSC format.  This function uses
+ * only the sparsity structure of A (not the numeric values), so the
+ * result can be reused across multiple numeric factorizations with
+ * QDLDL_numeric when the sparsity pattern of A remains the same.
+ *
+ * Does not use MALLOC.  It is assumed that Lp is allocated with n+1
+ * elements, Li is allocated with sumLnz elements (as returned by
+ * QDLDL_etree), and iwork is allocated with 3*n elements.
+ *
+ * @param  n      number of columns in L and A (both square)
+ * @param  Ap     column pointers (size n+1) for columns of A (not modified)
+ * @param  Ai     row indices of A.  Has Ap[n] elements (not modified)
+ * @param  Lp     column pointers (size n+1) for columns of L
+ * @param  Li     row indices of L.  Has sumLnz elements
+ * @param  Lnz    count of nonzeros in each column of L below diagonal,
+ *                as given by QDLDL_etree (not modified)
+ * @param  etree  elimination tree as given by QDLDL_etree (not modified)
+ * @param  bwork  working array of bools. Length is n
+ * @param  iwork  working array of integers. Length is 3*n
+ *
+ */
+QDLDL_API void QDLDL_symbolic(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* Ai,
+                              QDLDL_int* Lp, QDLDL_int* Li, const QDLDL_int* Lnz,
+                              const QDLDL_int* etree, QDLDL_bool* bwork, QDLDL_int* iwork);
+
+
+/**
+ * Compute the numeric LDL decomposition for a quasidefinite matrix,
+ * given a pre-computed symbolic factorization (sparsity pattern of L).
+ *
+ * This function computes the numeric values Lx, D, and Dinv using
+ * the sparsity pattern (Lp, Li) previously computed by QDLDL_symbolic.
+ * It can be called repeatedly with different Ax values (but the same
+ * sparsity pattern) to efficiently refactorize.
+ *
+ * Does not use MALLOC.
+ *
+ * @param  n      number of columns in L and A (both square)
+ * @param  Ap     column pointers (size n+1) for columns of A (not modified)
+ * @param  Ai     row indices of A.  Has Ap[n] elements (not modified)
+ * @param  Ax     data of A.  Has Ap[n] elements (not modified)
+ * @param  Lp     column pointers (size n+1) for columns of L (not modified,
+ *                as computed by QDLDL_symbolic)
+ * @param  Li     row indices of L (not modified, as computed by QDLDL_symbolic)
+ * @param  Lx     data of L.  Has Lp[n] elements
+ * @param  D      vectorized factor D.  Length is n
+ * @param  Dinv   reciprocal of D.  Length is n
+ * @param  Lnz    count of nonzeros in each column of L below diagonal,
+ *                as given by QDLDL_etree (not modified)
+ * @param  etree  elimination tree as given by QDLDL_etree (not modified)
+ * @param  bwork  working array of bools. Length is n
+ * @param  iwork  working array of integers. Length is 3*n
+ * @param  fwork  working array of floats. Length is n
+ * @return        Returns a count of the number of positive elements
+ *                in D.  Returns -1 and exits immediately if any element
+ *                of D evaluates exactly to zero (matrix is not quasidefinite
+ *                or otherwise LDL factorisable)
+ *
+ */
+QDLDL_API QDLDL_int QDLDL_numeric(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* Ai,
+                                  const QDLDL_float* Ax, const QDLDL_int* Lp, const QDLDL_int* Li,
+                                  QDLDL_float* Lx, QDLDL_float* D, QDLDL_float* Dinv,
+                                  const QDLDL_int* Lnz, const QDLDL_int* etree, QDLDL_bool* bwork,
+                                  QDLDL_int* iwork, QDLDL_float* fwork);
+
+
 /**
  * Compute an LDL decomposition for a quasidefinite matrix
  * in compressed sparse column form, where the input matrix is
@@ -102,6 +175,10 @@ QDLDL_API QDLDL_int QDLDL_etree(const QDLDL_int n, const QDLDL_int* Ap, const QD
  *
  * Returns factors L, D and Dinv = 1./D.
  *
+ * This is equivalent to calling QDLDL_symbolic followed by QDLDL_numeric.
+ * If you need to refactorize with different numeric values but the same
+ * sparsity pattern, use QDLDL_symbolic once and QDLDL_numeric repeatedly.
+ *
  * Does not use MALLOC.  It is assumed that L will be a compressed
  * sparse column matrix with data (n,Lp,Li,Lx)  with sufficient space
  * allocated, with a number of nonzeros equal to the count given
diff --git a/src/qdldl.c b/src/qdldl.c
@@ -91,11 +91,99 @@ QDLDL_int QDLDL_etree(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* A
 }
 
 
-QDLDL_int QDLDL_factor(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* Ai,
-                       const QDLDL_float* Ax, QDLDL_int* Lp, QDLDL_int* Li, QDLDL_float* Lx,
-                       QDLDL_float* D, QDLDL_float* Dinv, const QDLDL_int* Lnz,
-                       const QDLDL_int* etree, QDLDL_bool* bwork, QDLDL_int* iwork,
-                       QDLDL_float* fwork) {
+void QDLDL_symbolic(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* Ai, QDLDL_int* Lp,
+                    QDLDL_int* Li, const QDLDL_int* Lnz, const QDLDL_int* etree,
+                    QDLDL_bool* bwork, QDLDL_int* iwork) {
+    QDLDL_int   i = 0;
+    QDLDL_int   k = 0;
+    QDLDL_int   bidx = 0;
+    QDLDL_int   nextIdx = 0;
+    QDLDL_int   nnzE = 0;
+    QDLDL_int   nnzY = 0;
+    QDLDL_int   tmpIdx = 0;
+    QDLDL_int*  yIdx;
+    QDLDL_int*  elimBuffer;
+    QDLDL_int*  LNextSpaceInCol;
+    QDLDL_bool* yMarkers;
+
+    // Partition working memory into pieces
+    yMarkers = bwork;
+    yIdx = iwork;
+    elimBuffer = iwork + n;
+    LNextSpaceInCol = iwork + n * 2;
+
+    Lp[0] = 0; // First column starts at index zero
+
+    for(i = 0; i < n; i++) {
+        // Compute L column pointers (cumulative sum of Lnz)
+        Lp[i + 1] = Lp[i] + Lnz[i];
+
+        // Initialize markers and next-space trackers
+        yMarkers[i] = QDLDL_UNUSED;
+        LNextSpaceInCol[i] = Lp[i];
+    }
+
+    // Start from 1: column 0 of L has no subdiagonal entries
+    for(k = 1; k < n; k++) {
+        nnzY = 0;
+
+        tmpIdx = Ap[k + 1];
+
+        for(i = Ap[k]; i < tmpIdx; i++) {
+            bidx = Ai[i];
+
+            // Skip the diagonal entry
+            if(bidx == k) {
+                continue;
+            }
+
+            // Walk the elimination tree to find the nonzero pattern
+            nextIdx = bidx;
+
+            if(yMarkers[nextIdx] == QDLDL_UNUSED) {
+
+                yMarkers[nextIdx] = QDLDL_USED;
+                elimBuffer[0] = nextIdx;
+                nnzE = 1;
+
+                nextIdx = etree[bidx];
+
+                while(nextIdx != QDLDL_UNKNOWN && nextIdx < k) {
+                    if(yMarkers[nextIdx] == QDLDL_USED)
+                        break;
+
+                    yMarkers[nextIdx] = QDLDL_USED;
+                    elimBuffer[nnzE] = nextIdx;
+                    nnzE++;
+                    nextIdx = etree[nextIdx];
+                }
+
+                while(nnzE) {
+                    yIdx[nnzY++] = elimBuffer[--nnzE];
+                }
+            }
+        }
+
+        // Place the row indices into Li
+        for(i = (nnzY - 1); i >= 0; i--) {
+            QDLDL_int cidx = yIdx[i];
+            tmpIdx = LNextSpaceInCol[cidx];
+
+            Li[tmpIdx] = k;
+            LNextSpaceInCol[cidx]++;
+
+            // Reset marker
+            yMarkers[cidx] = QDLDL_UNUSED;
+        }
+    }
+}
+
+
+QDLDL_int QDLDL_numeric(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* Ai,
+                        const QDLDL_float* Ax, const QDLDL_int* Lp, const QDLDL_int* Li,
+                        QDLDL_float* Lx, QDLDL_float* D, QDLDL_float* Dinv, const QDLDL_int* Lnz,
+                        const QDLDL_int* etree, QDLDL_bool* bwork, QDLDL_int* iwork,
+                        QDLDL_float* fwork) {
     QDLDL_int    i = 0;
     QDLDL_int    j = 0;
     QDLDL_int    k = 0;
@@ -120,16 +208,7 @@ QDLDL_int QDLDL_factor(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int*
     LNextSpaceInCol = iwork + n * 2;
     yVals = fwork;
 
-
-    Lp[0] = 0; // First column starts at index zero
-
     for(i = 0; i < n; i++) {
-        // Compute L column indices
-        Lp[i + 1] = Lp[i] + Lnz[i]; // cumsum, total at the end
-
-        // Set all Yidx to be 'unused' initially
-        // in each column of L, the next available space
-        // to start is just the first space in the column
         yMarkers[i] = QDLDL_UNUSED;
         yVals[i] = 0.0;
         D[i] = 0.0;
@@ -228,7 +307,6 @@ QDLDL_int QDLDL_factor(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int*
             // Now I have the cidx^th element of y = L\b.
             // so compute the corresponding element of
             // this row of L and put it into the right place
-            Li[tmpIdx] = k;
             Lx[tmpIdx] = yVals_cidx * Dinv[cidx];
 
             // D[k] -= yVals[cidx]*yVals[cidx]*Dinv[cidx];
@@ -261,6 +339,18 @@ QDLDL_int QDLDL_factor(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int*
     return positiveValuesInD;
 }
 
+
+QDLDL_int QDLDL_factor(const QDLDL_int n, const QDLDL_int* Ap, const QDLDL_int* Ai,
+                       const QDLDL_float* Ax, QDLDL_int* Lp, QDLDL_int* Li, QDLDL_float* Lx,
+                       QDLDL_float* D, QDLDL_float* Dinv, const QDLDL_int* Lnz,
+                       const QDLDL_int* etree, QDLDL_bool* bwork, QDLDL_int* iwork,
+                       QDLDL_float* fwork) {
+
+    QDLDL_symbolic(n, Ap, Ai, Lp, Li, Lnz, etree, bwork, iwork);
+
+    return QDLDL_numeric(n, Ap, Ai, Ax, Lp, Li, Lx, D, Dinv, Lnz, etree, bwork, iwork, fwork);
+}
+
 // Solves (L+I)x = b
 void QDLDL_Lsolve(const QDLDL_int n, const QDLDL_int* Lp, const QDLDL_int* Li,
                   const QDLDL_float* Lx, QDLDL_float* x) {
diff --git a/tests/qdldl_tester.c b/tests/qdldl_tester.c
@@ -45,6 +45,7 @@ int ldl_factor_solve(QDLDL_int An, QDLDL_int* Ap, QDLDL_int* Ai, QDLDL_float* Ax
 #include "test_two_by_two.h"
 #include "test_zero_on_diag.h"
 #include "test_osqp_kkt.h"
+#include "test_symbolic.h"
 
 
 int tests_run = 0;
@@ -60,6 +61,8 @@ static char* all_tests() {
     mu_run_test(test_two_by_two);
     mu_run_test(test_zero_on_diag);
     mu_run_test(test_osqp_kkt);
+    mu_run_test(test_symbolic);
+    mu_run_test(test_symbolic_refactor);
 
     return 0;
 }
diff --git a/tests/test_symbolic.h b/tests/test_symbolic.h
