Seeking help for an abnormal result.

[D-teg]

Hello, Professor Tim Davis. I'm trying to solve a liner system for Ax = b, where A is a positively definite sparse matrix. I write the following code,and the result x is abnormal [1,1,1,1]. When AAx = b, the result holds. I should have an error in the call of the function and I would like to seek help with this. Looking forward for your message.

  cholmod_sparse *sparse = nullptr;
  cholmod_dense *x = nullptr;
  cholmod_dense *b = nullptr;
  cholmod_factor *factor = nullptr;
  cholmod_common common;
  cholmod_start(&common);
  unsigned n = 4;
  size_t nonZeros = 6;

  double A[n*n] = {
                    3.0, 0.0, 0.0, 1.0, 
                    0.0, 1.0, 0.0, 0.0,
                    0.0, 0.0, 2.0, 0.0,
                    1.0, 0.0, 0.0, 4.0,
                  };
  double B[n] = {17.0, 1.0, 4.0, 24.0};

  // Dense vector b
  b = cholmod_allocate_dense(n, 1, n, CHOLMOD_REAL, &common);
  memcpy(static_cast<double*> (b->x), B, sizeof(double) * n);

  // Sparse matrix A
  size_t index = 0;
  int Tp[n+1];
  int Ti[nonZeros];
  double Tx[nonZeros];
  for(size_t i  = 0; i < n;++i){
    Tp[i] = index;
    for(size_t j = 0; j< n; ++j){
      size_t id = i*n+j;
      if(A[id] >1e-20){
        Ti[index] = j;
        Tx[index] = A[id];
        index++;
      }
    }
  }
  Tp[n] = index;

  sparse = cholmod_allocate_sparse(n, n, nonZeros, true, true, 0, CHOLMOD_REAL, &common);
  if(sparse){
    memcpy(static_cast<int*> (sparse->p), Tp, sizeof(int) * (n+1));
    memcpy(static_cast<int*> (sparse->i), Ti, sizeof(int) * nonZeros);
    memcpy(static_cast<double*> (sparse->x), Tx, sizeof(double) * nonZeros);
  }
  cholmod_print_sparse(sparse, "A", &common);
  factor = cholmod_analyze(sparse, &common);
  cholmod_factorize(sparse, factor, &common);
  x = cholmod_solve(CHOLMOD_A, factor, b, &common);

  int size = n; 
  std::vector<double> result(size);
  std::copy_n(static_cast<double*>(x->x),size,result.data());
  for(int j = 0;j < size;++j)
    std::cout<<result[j]<<std::endl;
  if(factor)  cholmod_free_factor( &factor, &common);
  if(sparse)  cholmod_free_sparse( &sparse, &common);
  cholmod_free_dense( &x, &common);
  cholmod_free_dense( &b, &common);
  cholmod_finish( &common);

[D-teg]

my version is suitesparse-4.0.2-10.el7.x86_64.

[DrTimothyAldenDavis]

When you created the matrix, you gave it an "stype" of zero, which means "unsymmetric, with both lower and triangular parts stored". That is a parameter in cholmod_allocate_sparse.

If you ask for a Cholesky factorization of such a matrix, I factorize A*A'. I use that feature for handling linear programming problems.

if you want to factorize and solve Ax=b and if A is symmetric, you have to use a different A->stype.

[D-teg]

Thank you for your suggestion, I misunderstood the meaning of this parameter.Thanks again. : )