Fix for perm col case

include/qdldl.h

@@ -10,7 +10,7 @@
 #ifndef QDLDL_TYPES_DEFINED
   #include <stdbool.h>
   #define QDLDL_TYPES_DEFINED
-  typedef long long      QDLDL_bool; //DEBUG
+  typedef bool      QDLDL_bool;
   typedef long long QDLDL_int;
   typedef double    QDLDL_float;
 # endif
@@ -36,15 +36,23 @@ extern "C" {
   * The data in (n,Ap,Ai) are from a square matrix A in CSC format, and
   * should include the upper triangular part of A only.
   *
+  * This function is only intended for factorisation of QD matrices specified
+  * by their upper triangular part.  An error is returned if any column has
+  * data below the diagonal or s completely empty.
+  *
+  * For matrices with a non-empty column but a zero on the corresponding diagonal,
+  * this function will *not* return an error, as it may still be possible to factor
+  * such a matrix in LDL form.   No promises are made in this case though...
+  *
   * @param   n     number of columns in CSC matrix A (assumed square)
   * @param  Ap     column pointers (size n+1) for columns of A
   * @param  Ai     row indices of A.  Has Ap[n] elements
   * @param  work   work vector (size n) (no meaning on return)
   * @param  Lnz    count of nonzeros in each column of L (size n) below diagonal
   * @param  etree  elimination tree (size n)
   * @return total  sum of Lnz (i.e. total nonzeros in L below diagonal). Returns
-  *                -1 if the input does not have triu structure or is missing entries
-  *                on the diagonal
+  *                -1 if the input does not have triu structure or has an empty
+  *                column.
   *
   *
 */
@@ -85,7 +93,8 @@ extern "C" {
   * @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)
+  *                of D evaluates exactly to zero (matrix is not quasidefinite
+  *                or otherwise LDL factorisable)
   *
 */
 

src/qdldl.c

@@ -5,32 +5,6 @@
 #define UNUSED 0
 
 
-//DEBUG
-#include <stdio.h>
-void qdprint_arrayi(const QDLDL_int* data, QDLDL_int n,char* varName){
-
-  QDLDL_int i;
-  printf("%s = [",varName);
-  for(i=0; i< n; i++){
-    printf("%lli,",data[i]);
-  }
-  printf("]\n");
-
-}
-
-void qdprint_arrayf(const QDLDL_float* data, QDLDL_int n, char* varName){
-
-  QDLDL_int i;
-  printf("%s = [",varName);
-  for(i=0; i< n; i++){
-    printf("%.3g,",data[i]);
-  }
-  printf("]\n");
-
-}
-// END DEBUG
-
-
 /* Compute the elimination tree for a quasidefinite matrix
    in compressed sparse column form.
 */
@@ -112,11 +86,6 @@ QDLDL_int QDLDL_factor(const QDLDL_int    n,
   yVals           = fwork;
 
 
-  qdprint_arrayi(Ap,n,"Ap");
-  qdprint_arrayi(Ai,Ap[n],"Ai");
-  qdprint_arrayf(Ax,Ap[n],"Ax");
-
-
   Lp[0] = 0; //first column starts at index zero
 
   for(i = 0; i < n; i++){
@@ -149,28 +118,22 @@ QDLDL_int QDLDL_factor(const QDLDL_int    n,
     //This loop determines where nonzeros
     //will go in the kth row of L, but doesn't
     //compute the actual values
-    printf("\n\nIn loop : k = %lli\n",k);
     tmpIdx = Ap[k+1];
 
-    printf("loop range : [%lli,%lli)\n",Ap[k],tmpIdx);
-
     for(i = Ap[k]; i < tmpIdx; i++){
 
-      printf("Inner loop : i = %lli\n",i);
+      bidx        = Ai[i];   // we are working on this element of b
 
       //Initialize D[k] as the element of this column
-      //corresponding to the diagonal place
-      if(Ai[i] == k){
+      //corresponding to the diagonal place.  Don't use
+      //this element as part of the elimination step
+      //that computes the k^th row of L
+      if(bidx == k){
         D[k] = Ax[i];
-        printf("Continuing loop : D[k] = %f\n",Ax[i]);
-        if(D[k] == 0.0){return -1;}
         continue;
       }
 
-      bidx        = Ai[i];   // we are working on this element of b
       yVals[bidx] = Ax[i];   // initialise y(bidx) = b(bidx)
-      printf("i = %lli : bidx = %lli : ",i,bidx);
-      qdprint_arrayf(yVals,n,"yVals");
 
       // use the forward elimination tree to figure
       // out which elements must be eliminated after
@@ -204,13 +167,6 @@ QDLDL_int QDLDL_factor(const QDLDL_int    n,
 
     } //end for i
 
-    //DEBUG
-    printf("k = %lli : ",k);
-    qdprint_arrayf(D,n,"D");
-    qdprint_arrayf(yVals,n,"yVals");
-    qdprint_arrayi(yMarkers,n,"yMarkers");
-    //END DEBUG
-
     //This for loop places nonzeros values in the k^th row
     for(i = (nnzY-1); i >=0; i--){
 
@@ -245,11 +201,7 @@ QDLDL_int QDLDL_factor(const QDLDL_int    n,
     //Maintain a count of the positive entries
     //in D.  If we hit a zero, we can't factor
     //this matrix, so abort
-    if(D[k] == 0.0){
-      printf("Abort case: k = %lli : ",k);
-      qdprint_arrayf(D,n,"D");
-      return -1;
-    }
+    if(D[k] == 0.0){return -1;}
     if(D[k]  > 0.0){positiveValuesInD++;}
 
     //compute the inverse of the diagonal

tests/julia/test.jl

@@ -124,10 +124,6 @@ A = sparse([0.0 5.0; 5.0 1.0])
 b = randn(c_float,A.n);
 println("Zero at A[0,0]         : ",  LDLsolve(triu(A),b))
 
-A = sparse([1.0 5.0; 5.0 0.0])
-b = randn(c_float,A.n);
-println("Zero at A[end,end]     : ",  LDLsolve(triu(A),b))
-
 A = sparse([1.0 1.0; 1.0 1.0])
 b = randn(c_float,A.n);
 println("Rank deficient matrix  : ",  LDLsolve(triu(A),b))
@@ -147,6 +143,14 @@ A = qdMatExample();
 b = randn(c_float,A.n);
 println("Basic example                : ",  norm(A\b-LDLsolve(triu(A),b)))
 
+A = sparse([4.0 1.0 2.0; 1.0 0.0 1.0;2.0 1.0 -3.0])
+b = [6.0,9.0,12.0];
+println("Zero mid matrix              : ",  norm(A\b-LDLsolve(triu(A),b)))
+
+A = sparse([1.0 5.0; 5.0 0.0])
+b = randn(c_float,A.n);
+println("Zero at A[end,end]           : ",  norm(A\b-LDLsolve(triu(A),b)))
+
 A = sparse([1.0 1.0;1.0 -1.0]);
 b = cumsum(ones(c_float(A.n)));
 println("2x2                          : ", norm(A\b- LDLsolve(triu(A),b)))
@@ -197,6 +201,3 @@ Ax = [-0.25; -0.25; 1; 0.513578; 0.529142; -0.25; -0.25; 1.10274; 0.15538; 1.258
 A = sparse(Ai, Aj, Ax)
 b = randn(c_float,A.n);
 println("tricky permuted osqp matrix  : ", norm((A + A' - diagm(diag(A)))\b- LDLsolve(A, b)))
-
-
-

tests/qdldl_tester.c

@@ -36,7 +36,7 @@ static char* all_tests() {
   mu_run_test(test_sym_structure);
   mu_run_test(test_tril_structure);
   mu_run_test(test_two_by_two);
-  //mu_run_test(test_zero_on_diag);
+  mu_run_test(test_zero_on_diag);
   mu_run_test(test_osqp_kkt);
 
   return 0;

tests/test_zero_on_diag.h

@@ -7,14 +7,15 @@ static char* test_zero_on_diag()
   QDLDL_int An = 3;
 
   // RHS and solution to Ax = b
-  QDLDL_float b[]    = {1,2,3};
-  QDLDL_float xsol[] = {4,-11,-2};
+  QDLDL_float b[]    = {6,9,12};
+  QDLDL_float xsol[] = {17,-46,-8};
 
   //x replaces b during solve (should fill due to zero in middle)
   //NB : this system is solvable, but not by LDL
   int status = ldl_factor_solve(An,Ap,Ai,Ax,b);
 
-  mu_assert("Zero diag Factorisation failed", status < 0);
+  mu_assert("Factorisation failed", status >= 0);
+  mu_assert("Solve accuracy failed", vec_diff_norm(b,xsol,An) < QDLDL_TESTS_TOL);
 
   return 0;
 }