CHOLMOD sparse dense multiplication issue

242 Views Asked by At

I am testing out cholmod to solve a sparse matrix system. As b of the Ax=b equation I want to use something that will lead to a known result, so I choose to declare as xe a vector of ones and b=A*xe.

I tried some different approaches to do the sparse-matrix dense-vector multiplication without any success. As first test I am using the ex15 matrix from SuitSparse matrix collection

#include "cholmod.h"
#include <iostream>
#pragma warning(disable:4996)
int main()
{
    const char* matFile = "C:/www/CholeskyMatGPU/ex15.mtx";
    const char* bfile = "C:/www/CholeskyMatGPU/b.mtx";
    FILE* bp = fopen(bfile, "w");
    FILE* fp = fopen(matFile, "r");
    cholmod_sparse *A;
    cholmod_dense *x, *xe, *b;
    cholmod_factor *L;

    cholmod_common* c = (cholmod_common*)malloc(sizeof(cholmod_common));
    cholmod_l_start(c); /* start CHOLMOD */
    c->useGPU = 0;
    c->supernodal = CHOLMOD_SUPERNODAL;

    A = cholmod_l_read_sparse(fp, c); /* read in a matrix */
    cholmod_l_print_sparse(A, "A", c); /* print the matrix */
    fclose(fp);

    double alpha[2] = { 1., 0. };
    double beta[2] = { 0., 0. };
    xe = cholmod_l_ones(A->nrow, 1, A->xtype, c); 
    b = cholmod_l_allocate_dense(A->nrow, 1, A->nrow, CHOLMOD_LONG, c); //attempt1
    b = cholmod_copy_dense(xe, c); //attempt2
    b = cholmod_l_zeros(A->nrow, 1, A->xtype, c); //attempt3
    cholmod_sdmult(A, 0, alpha, beta, xe, b, c); //b = A*xe
    cholmod_l_write_dense(bp, b, "", c);
    L = cholmod_l_analyze(A, c); /* analyze */
    cholmod_l_factorize(A, L, c); /* factorize */
    x = cholmod_l_solve(CHOLMOD_A, L, b, c); /* solve Ax=b */
}

I expect b to be the result of the A*xe operation, so something like

-4487487.55116591
228083519.121975
0.0561846979753665
18213301.5299736
...
30383411.1062435
8073.91311889887
0.00602254606925790
-7101038.89577146
378419792.843878
0.0931721845637848

(checked from matlab) but the shown code will print a vector of zeros, there is something wrong with what I am doing with sdmult.

A = Problem.A;
n = length(A);
xe(1:n,1) = 1;
b=A*xe;
R=chol(A);
tmp=(R'\b);
x=R\tmp;
1

There are 1 best solutions below

0
David On BEST ANSWER

The issue was an increadibly stupid mistake, there isn't the long prefix in the method name! This won't throw any error but it won't execute the multiplication. With cholmod_l_sdmult(A, 0, alpha, beta, xe, b, c); instead of cholmod_sdmult(A, 0, alpha, beta, xe, b, c); everything works fine..