title: CHOLMOD sparse matrix cholesky decomposition: incorrect factor?
body: <p>I have been using CHOLMOD to factorise the matrix A and solve the system Ax = b, for A being the Hessian matrix (printed below) and b = [1, 1, 1] created by the <code>cholmod_ones</code> function.</p>
  
  <p>Unfortunately, the solution for x is incorrect (should be [1.5, 2.0, 1.5]) and to confirm I then multiplied A and x back together and don't get [1, 1, 1]. I don't quite understand what I am doing wrong.</p>
  
  <p>Additionally, I've looked at the factor and the values of the matrix elements don't make sense either.</p>
  
  <p><strong>Output</strong></p>
  
  <pre><code>Hessian:
     2.000   -1.000    0.000 
    -1.000    2.000   -1.000 
     0.000   -1.000    2.000 
  Solution:
     2.500    0.000    0.000 
     3.500    0.000    0.000 
     2.500    0.000    0.000 
  B vector:
     1.500    0.000    0.000 
     2.000    0.000    0.000 
     1.500    0.000    0.000 
  </code></pre>
  
  <p><strong>Code</strong></p>
  
  <p><code>iterate_hessian()</code> is an external function that returns doubles which are read into the CHOLMOD hessian matrix.</p>
  
  <p>The entry point for the code is <code>cholesky_determinant</code> which is called with an argument which gives the dimension of the (square) matrix.</p>
  
  <pre class="lang-c prettyprint-override"><code>#include &lt;cholmod.h&gt;
  #include &lt;string.h&gt;
  
  // Function prototype that gives the next value of the Hessian
  double iterate_hessian();
  
  cholmod_sparse *cholmod_hessian(double *hessian, size_t dimension, cholmod_common *common) {
    // This function assigns the Hessian matrix from OPTIM to a dense matrix for CHOLMOD to use.
    // Allocate a dense cholmod matrix of appropriate size
    cholmod_triplet *triplet_hessian;
    triplet_hessian = cholmod_allocate_triplet(dimension, dimension, dimension*dimension, 0, CHOLMOD_REAL, common);
    // Loop through values of hessian and assign their row/column index and values to triplet_hessian.
    size_t loop;
    for (loop = 0; loop &lt; (dimension * dimension); loop++) {
      if (hessian[loop] == 0) {
        continue;
      }
      ((int*)triplet_hessian-&gt;i)[triplet_hessian-&gt;nnz] = loop / dimension;
      ((int*)triplet_hessian-&gt;j)[triplet_hessian-&gt;nnz] = loop % dimension;
      ((double*)triplet_hessian-&gt;x)[triplet_hessian-&gt;nnz] = hessian[loop];
      triplet_hessian-&gt;nnz++;
    }
    // Convert the triplet to a sparse matrix and return.
    cholmod_sparse *sparse_hessian;
    sparse_hessian = cholmod_triplet_to_sparse(triplet_hessian, (dimension * dimension), common);
    return sparse_hessian;
  }
  
  void print_matrix(cholmod_dense *matrix, size_t dimension) {
    // matrix-&gt;x is a void pointer, so first copy it to a double pointer
    // of an appropriate size
    double *y = malloc(sizeof(matrix-&gt;x));
    y = matrix-&gt;x;
    // Loop variables
    size_t i, j;
    // Row
    for(i = 0; i &lt; dimension; i++) {
    // Column
      for(j = 0; j &lt; dimension; j++) {
        printf("% 8.3f ", y[i + j * dimension]);
      }
      printf("\n");
    }
  }
  
  cholmod_dense *factorized(cholmod_sparse *sparse_hessian, cholmod_common *common) {
    cholmod_factor *factor;
    factor = cholmod_analyze(sparse_hessian, common);
    cholmod_factorize(sparse_hessian, factor, common);
    cholmod_dense *b, *x;
    b = cholmod_ones(sparse_hessian-&gt;nrow, 1, sparse_hessian-&gt;xtype, common);
    x = cholmod_solve(CHOLMOD_LDLt, factor, b, common);
    cholmod_free_factor(&amp;factor, common);
    // Return the solution, x
    return x;
  }   
  
  double cholesky_determinant(int *dimension) {
    // Declare variables
    double determinant;
    cholmod_sparse *A;
    cholmod_dense *B, *Y;
    cholmod_common common;
    // Start CHOLMOD
    cholmod_start (&amp;common);
    // Allocate storage for the hessian (we want to copy it)
    double *hessian = malloc(*dimension * *dimension * sizeof(hessian)); 
    // Get the hessian from OPTIM
    int i = 0;
    for (i = 0; i &lt; (*dimension * *dimension); i++) {
      hessian[i] = iterate_hessian();
    }
    A = cholmod_hessian(hessian, *dimension, &amp;common);
    printf("Hessian:\n");
    print_matrix(cholmod_sparse_to_dense(A, &amp;common), *dimension);
    B = factorized(A, &amp;common);
    printf("Solution:\n");
    print_matrix(B, *dimension);
    double alpha[] = {1, 0};
    double beta[] = {0, 0};
    Y = cholmod_allocate_dense(*dimension, 1, *dimension, CHOLMOD_REAL, &amp;common);
    cholmod_sdmult(A, 0, alpha, beta, B, Y, &amp;common);
    printf("B vector:\n");
    print_matrix(Y, *dimension);
    determinant = 0.0;
    // Free up memory and finish CHOLMOD
    cholmod_free_sparse (&amp;A, &amp;common);
    cholmod_free_dense (&amp;B, &amp;common);
    cholmod_finish (&amp;common);           
    return determinant;
  }
  </code></pre>