Data error in triangular factorization

[jrk1982]

Seems like LU factorization works only for full rank (row or column rank) matrices...

Is this a limitation of the beta, or is this an unimplemented feature?
Below is the example code that I use to run... 32Bit FC10, GTX295 Cuda 2.3

=============================================================

// This example is meant to show a typical program flow when using CULA.
int main(int argc, char** argv)
{
#ifdef NDEBUG
int M = 3;
#else
int M = 3;
#endif
int N = M;
int LWORK = 0;
culaStatus status;

float* A = NULL;
int* TAU = NULL;

A = (float*)malloc(M*N*sizeof(float));
TAU = (int*)malloc((M<N?M:N)*sizeof(int));
uint i,j;
for(i=0;i<M;i++)
for(j=0;j<N;j++)
{
*(A+j*M+i) = (float)(i+j+((int)rand())%10);
printf("%f\t",*(A+j*M+i));
}

if(!A || !TAU)
{
free(A);
free(TAU);
return EXIT_FAILURE;
}

printf("Initializing CULA\n");
status = culaInitialize();
checkStatus(status);

memset(A, 1, M*N*sizeof(float));

printf("Calling SGETRF\n");
status = culaSgetrf(M, N, A, M, TAU);
checkStatus(status);
culaShutdown();

free(A);
free(TAU);

return EXIT_SUCCESS;
}

[john]

The LAPACK-style LU factorization doesn't work on rank deficient matrices; it considers them singular and rather than diving by zero returns a code to indicate singularity. Receiving a culaDataError is equivalent to a "info" code > 0 in other implementations. In this respect, we should be providing identical functionality.

For more information, please see: http://www.intel.com/software/products/ ... getrf.html

It's worth noting that getrf() is valid on non-square matrices (both over- and under-determined) as long as they're still full rank.

I hope this answers your question. If I have misunderstood, can you please reply back with more information.

Thanks,
John

[jrk1982]

Sure. Thanks for the clarification.

But I still get Data Error with the code snipped I had attached. The matrix input is full rank I believe...

Also, MATLAB behaves differently(well) for singular matrix inputs?!

Thanks for looking into this matter!

[john]

I see a quick error in your code:
memset(A, 1, M*N*sizeof(float));

This line sets your entire matrix to 2.36942e-38 for every value, which is rank deficient.

[jrk1982]

Oops. Good catch. I ought to have to have removed that :) Thanks.

[jrk1982]

Thanks. That worked for the LU factors. However, 1 last issue which I think is genuine.

The permutation vector returned does not seem right... I get [3 3 3] for input, [3 7 9; 6 5 8; 8 5 13]. Matlab returns [3 1 2]...

[john]

That's just a notational difference. Matlab's vector is meant to be evaluated in parallel (ie A(p,: )) and LAPACK's is procedural (apply this swap, THEN do this swap, THEN do this one, etc.) In this case they're equivalent.

[jrk1982]

Cool. Thanks!

[kolonel]

Hi,
I run this code (I removed the memset), but I cannot get the permutation vector as [3 3 3]. TAU vector here returns weird digits as 4.204e-045#DEN, any idea?!
Thanks a lot