jluttine / cholmod-extra

Extra functions for SuiteSparse CHOLMOD module.
GNU General Public License v2.0
0 stars 2 forks source link

Problem with getting the inverse #1

Closed LindleyLentati closed 11 years ago

LindleyLentati commented 11 years ago

Hi there,

I'm trying to use your code to calculate the inverse of a matrix in combination with Cholmod, i have the following:

int main(void) { int N = 10 ; int i, j, n ; //Int _I, J ; //double X ; int nz = 0; double Ax ; double x, error ; cholmod_dense A, invK, spinvK, I ; cholmod_sparse K, V ; cholmod_factor L ; cholmod_common Common,_cm ; clock_t start, end; double cpu_time_used; cm=&Common; // Start using CHOLMOD cholmod_start(cm) ; cm->print=5; /* SPARSE COVARIANCE MATRIX CONSTRUCTION */

// Generate random symmetric positive (semi)definite matrix
A = cholmod_zeros(N, N, CHOLMOD_REAL, &Common) ;
Ax =(double*) A->x ;
nz = N ;
// Make positive-definite by adding something positive to the
// diagonal
for (n = 0; n < N; n++)
{
    Ax[n+n*N] += 5;
}
cholmod_print_dense(A,"A",&Common);
// Make the matrix sparse
K = cholmod_dense_to_sparse(A, TRUE, &Common) ;

// Identity matrix
I = cholmod_eye(N,N,CHOLMOD_REAL,&Common) ;

/* SIMPLICIAL */

// Factorize
Common.supernodal = CHOLMOD_SIMPLICIAL ;
L = cholmod_analyze(K, &Common) ;
cholmod_factorize(K, L, &Common) ;

// Compute the sparse inverse and the full inverse
start = clock();
V = cholmod_spinv(L, &Common) ;

end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
invK = cholmod_solve(CHOLMOD_A, L, I, &Common) ;
spinvK = cholmod_sparse_to_dense(V, &Common) ;
cholmod_print_dense(spinvK,"V",&Common);

}

I.e. just a diagonal matrix, with 5 along the diagonal, which should give me 1/5, but it gives me 1/5^2 and I have no idea why! Have you come across this kind of issue before? Am i just doing something stupid.

Cheers Lindley

jluttine commented 11 years ago

Thanks for reporting this issue! If the sparse matrix K has stype==0, then it is considered as unsymmetric matrix. From CHOLMOD manual: "In the symmetric case, A or A(p,p) is analyzed, where p is the fill-reducing ordering. In the unsymmetric case, A_A’ or A(p,:)_A(p,:)’ is analyzed." Thus, your code computes the inverse of K*K'. To fix your code, it should sufficient to have a line K->stype = 1; // Use upper triangular part at some point between cholmod_dense_to_sparse and cholmod_analyze. I'll fix the similar code in the test function to set the stype to a positive (non-zero) value as well. Thanks!

jluttine commented 11 years ago

There's also issue1.c in Source folder (and compiled by Makefile to Build folder), which should behave correctly.