matrixfunctions / GraphMatFun.jl

Computation graphs for matrix functions
MIT License
13 stars 0 forks source link

Code gen LangC identity matrix lincomb #17

Closed jarlebring closed 3 years ago

jarlebring commented 3 years ago

When using lincomb with an identity matrix it seems LangC causes problems.

Considering that ONE is declared as

    const double ZERO = 0.0;

I believe, this use of ONE should be something different:

https://github.com/jarlebring/GraphMatFun.jl/blob/2558b55495b7bafeaa0b6082105546407b83352f/src/gen_code.jl#L584-L589

which generates code like

    memcpy(memslots+(1-1)*n*n, A, n*n*sizeof(*memslots));
    cblas_dscal(n*n, coeff2, memslots+(1-1)*n*n, 1);
    cblas_daxpby(n*n, coeff1, &ONE, 0,
             ONE, memslots+(1-1)*n*n, n+1);

It segfaults.

The third parameter in axbpy should be a matrix/vector.

I wouldn't mind looking it myself with a hint of what was the plan

jarlebring commented 3 years ago

Maybe &ONE should be a &onevec which is preallocated as a vector of ones?

jarlebring commented 3 years ago

If so, probably the incx (fourth argument) should be 1 and the size (first argument) should be n.

jarlebring commented 3 years ago

Wait! It was much more clever than I thought. incx=0 is very clever. The bug is probably just that n*n should be n (and keeping everything as before).

    cblas_daxpby(n*n, coeff1, &ONE, 0,
             ONE, memslots+(1-1)*n*n, n+1);

should be

    cblas_daxpby(n, coeff1, &ONE, 0,
             ONE, memslots+(1-1)*n*n, n+1);

Do you agree @mfasi ?

mfasi commented 3 years ago

Do you agree @mfasi ?

Yes, there was definitely a problem there. I changed n*n to n in a bunch of places, hopefully that solved the segmentation fault. The matrices I was using in my tests were probably too small to hit that issue.

jarlebring commented 3 years ago

Yes. Problem solved. :+1: