title: CHOLMOD supernodal factorization fails when sparse matrix becomes too dense
body: <p>I am using CHOLMOD in SuiteSparse to factor an <code>N</code> by <code>N</code> large band-diagonal matrix that is relatively sparse, ie, it contains only a few diagonals that are non-zero. The sparsity of the matrix is set by a covariance length parameter <code>l</code>. The larger <code>l</code> the greater number of off-diagonal elements that are non-zero. </p>
  
  <p>When <code>l</code> becomes large and many elements are non-zero, the supernodal CHOLMOD factorization suddenly begins to fail with the error message that "CHOLMOD warning: matrix not positive definite." From checking the math with a separate Python implementation, I know that the matrix should be positive definite. Moreover, when I change from  </p>
  
  <pre><code>Common-&gt;.supernodal = CHOLMOD_SUPERNODAL; 
  </code></pre>
  
  <p>to</p>
  
  <pre><code>Common-&gt;supernodal = CHOLMOD_SIMPLICIAL;
  </code></pre>
  
  <p>the factorization will then succeed. For the following code example below, the supernodal factorization will succeed as long as <code>l &lt; 2.5</code>. If I increase <code>l &gt;= 3.0</code>, I get the error that the matrix is not positive definite. Yet, if I then choose CHOMOD_SIMPLICIAL, the factorization will succeed. Can anyone help me identify why CHOLMOD_SUPERNODAL suddenly fails above a certain sparsity/density? Thank you!</p>
  
  <pre><code>//file is cov.c
  //Compile with $ gcc -Wall -o cov -Icholmod cov.c -lm -lcholmod -lamd -lcolamd -lblas -llapack -lsuitesparseconfig
  #include &lt;stdio.h&gt;
  #include &lt;math.h&gt;
  #include "cholmod.h"
  
  #define PI 3.14159265
  
  float k_3_2 (float r, float a, float l, float r0)
  {
      return (0.5 + 0.5 * cos(PI * r/r0)) * pow(a, 2.) * (1 + sqrt(3) * r/l) * exp(-sqrt(3) * r/l) ;
  }
  
  // function to initialize an array of increasing wavelengths
  void linspace (double *wl, int N, double start, double end)
  {
      double j; //double index
      double Ndist = (double) N;
      double increment = (end - start)/(Ndist -1.);
  
      int i;
      for (i = 0; i &lt; N; i++) {
          j = (double) i;
      wl[i] = start + j * increment;
      }
  }
  
  // create and return a sparse matrix using a wavelength array and parameters
  // for a covariance kernel.
  cholmod_sparse *create_sparse(double *wl, int N, double a, double l, cholmod_common *c)
  {
      double r0 = 6.0 * l; //Beyond r0, all entries will be 0
      //Pairwise calculate all of the r distances
      int i = 0, j = 0;
      double r;
  
      //First loop to determine the number non-zero elements
      int M = 0; //number of non-zero elements
      for (i = 0; i &lt; N; i++)
      {
          for (j = 0; j &lt; N; j++)
          {
          r = fabs(wl[i] - wl[j]);
          if (r &lt; r0) //Is the separation below our cutoff?
              M++;
          }
      }
  
      /* Initialize a cholmod_triplet matrix, which we will subsequently fill with
       * values. This matrix is NxN sparse with M total non-zero elements. 1 means we
       * want a square and symmetric matrix. */
      cholmod_triplet *T = cholmod_allocate_triplet(N, N, M, 1, CHOLMOD_REAL, c);
  
      if (T == NULL || T-&gt;stype == 0)         /* T must be symmetric */
      {
        cholmod_free_triplet (&amp;T, c) ;
        cholmod_finish (c) ;
        return (0) ;
      }
  
      //Do the loop again, this time to fill in the matrix
      int * Ti = T-&gt;i;
      int * Tj = T-&gt;j;
      double * Tx = T-&gt;x;
      int k = 0;
  
      //This time, only fill in the lower entries (and diagonal). 
      for (i = 0; i &lt; N; i++)
      {
          for (j = 0; j &lt;= i; j++)
          {
          r = fabs(wl[i] - wl[j]);
  
              if (r &lt; r0) //If the distance is below our cutoff, initialize
              {
                  Ti[k] = i;
                  Tj[k] = j;
                  Tx[k] = k_3_2(r, a, l, r0);
                  k++;
              }
          }
      }
      T-&gt;nnz = k;
  
      //The conversion will transpose the entries and add to the upper half.
      cholmod_sparse *A = cholmod_triplet_to_sparse(T, k, c);
      cholmod_free_triplet(&amp;T, c);
      return A;
  
  }
  
  int main(void)
  {
      //Create a sample array of wavelengths for testing purposes. 
      int N = 3000; 
      double wl[N];
      linspace(wl, N, 5100., 5200.); //initialize with wavelength values
  
      cholmod_common c ; //c is actually a struct, not a pointer to it.
      cholmod_start (&amp;c) ; // start CHOLMOD
      c.print = 5;
  
      //c.supernodal = CHOLMOD_SIMPLICIAL;
      c.supernodal = CHOLMOD_SUPERNODAL;
  
      float l = 2.5;
      cholmod_sparse *A = create_sparse(wl, N, 1.0, l, &amp;c);
  
      cholmod_factor *L ;
      L = cholmod_analyze (A, &amp;c) ;           
  
      cholmod_factorize (A, L, &amp;c) ;          
      printf("L-&gt;minor = %d\n", (int) L-&gt;minor);
  
      cholmod_dense *b, *x, *r;
  
      //Create a vector with the same number of rows as A
      b = cholmod_ones (A-&gt;nrow, 1, A-&gt;xtype, &amp;c) ;   // b = ones(n,1) 
      x = cholmod_solve (CHOLMOD_A, L, b, &amp;c) ;       // solve Ax=b    
      //the reason these are length two is because they can be complex
      double alpha [2] = {1,0}, beta [2] = {0,0} ;        // basic scalars 
  
      r = cholmod_copy_dense (b, &amp;c) ;            // r = b 
      cholmod_sdmult (A, 0, alpha, beta, x, r, &amp;c) ;      // r = Ax 
  
      cholmod_print_dense(r, "r", &amp;c); //This should be equal to b
  
      cholmod_free_sparse(&amp;A, &amp;c); // free all of the variables
      cholmod_free_factor(&amp;L, &amp;c);
      cholmod_free_dense(&amp;b, &amp;c);
      cholmod_free_dense(&amp;x, &amp;c);
      cholmod_free_dense(&amp;r, &amp;c);
      cholmod_finish (&amp;c) ;  // finish CHOLMOD 
  
      return (0) ;
  }
  </code></pre>