Still problem with ZGETRF and ZGETRI

[Dan4]

Dear all,

I had posted another topic regarding my problem with ?GETRF and ?GETRI and ?GETRS. The problems seemed to be solve, but I still get very strange results, sometimes. The results are totally different from results I get from another library just in order to verify that I am using MKL correctly. I am using Linux, Intel C++ Compiler 11.1 and MKL. The code I used is below:



dense_matrix RL2(nCells, nCells);

double *MKL_RL2;
MKL_RL2=new double[nCells*nCells];

double x=1.0;

for (int i=0; i for (int j=0; j {
RL2(i, j)=x;
MKL_RL2[j*nCells+i]=x;

x+=1.0;
}

try
{
gmm::lu_inverse(RL2);

// MKL solution to perform "RL" matrix inversion
int info;
int ipiv[nCells];
int lwork=-1;
double optwork;

dgetrf(&nCells, &nCells, MKL_RL2, &nCells, ipiv, &info);
dgetri(&nCells, MKL_RL2, &nCells, ipiv, &optwork, &lwork, &info);

double *work;
lwork=(int)optwork;
work=new double[lwork];

dgetri(&nCells, MKL_RL2, &nCells, ipiv, work, &lwork, &info);

delete[] work;

}
catch (dal::failure_error& e)
{
DAL_THROW(failure_error, "Failed to invert RL matrix\\n" << e.what());
}

delete[] MKL_RL2;



As you see, I am comparing "RL2" from GMM++ library with "MKL_RL2" with MKL and as I stated, the results are 100% different. The results from GMM++ are correct (I verified with MATLAB) but MKL, by some reason doesn't work correctly. I have same problem with I use MKL function with complex numbers.

Is there anything wrong with the code that I call MKL functions? First I factorize the matrix, then query to get optimal value for the work space, and the call getri to invert the matrix. This is exactly as MKL example in FORTRAN does (except querying and they used 64*N which I tried also without any success).

I would appreciate your helps.

Regards,

Dan

[mecej4]

The matrix that you use is singular. For example, with nCells=3, it is

[1 4 7]

[2 5 8]

[3 6 9]

Please give a specific example, preferably involving a small matrix such as this one, but one that has an inverse, where MKL gives you an incorrect result for the inverse.

[Dan4]

Yes you are right. The previous code can generate singular matrices in some cases. So I changed it that the matrices are filled randomly but there is a huge difference between results:

dense_matrix RL2(nCells, nCells);

double *MKL_RL2;
MKL_RL2=new double[nCells*nCells];

srand((unsigned)time(0));

for (int i=0; i for (int j=0; j {
RL2(i, j)=(double)(rand()%100)+1;
MKL_RL2[j*nCells+i]=(double)(rand()%100)+1;
}


try
{
gmm::lu_inverse(RL2);

// MKL solution to perform "RL" matrix inversion
int info;
int ipiv[nCells];
int lwork=-1;
double optwork;

dgetrf(&nCells, &nCells, MKL_RL2, &nCells, ipiv, &info);
dgetri(&nCells, MKL_RL2, &nCells, ipiv, &optwork, &lwork, &info);

double *work;
lwork=(int)optwork;
work=new double[lwork];

dgetri(&nCells, MKL_RL2, &nCells, ipiv, work, &lwork, &info);

delete[] work;

}
catch (dal::failure_error& e)
{
DAL_THROW(failure_error, "Failed to invert RL matrix\n" << e.what());
}

delete[] MKL_RL2;


Do I call MKL functions incorrectly? Is there something else that I need to take care but I haven't?

Thanks,

Dan

[mecej4]

Using a random matrix without reporting the values of the matrix elements practically ensures that nobody can reproduce your results.

One major problem in your code: you do not seem to realize that random numbers are, if may I say it, random!

Since you make different calls to rand() for generating elements of matrices RL2 and MKL_RL2, you are calling GMM and MKL with different matrices. In this situation, if you obtained identical inverse matrices from MKL and some other package, that would be a miracle, or rand() would have to possess the highly undesirable property of returning numbers in identical pairs!

Use, instead, a well-behaved and well-defined matrix such as (in Matlab notation) A =[2 -1 0; -1 2 -1; 0 -1 2], with inverse [3/4 1/2 1/4; 1/2 1 1/2; 1/4 1/2 3/4]. Report what you get from MKL for the inverse of this matrix with your code.

With such a simple matrix, one can deal with MKL and how to call it correctly, without bringing in yet another package (such as GMM, which I have never used) into the picture.

[Todd_R_Intel]

Dan,

Are you still seeing problems? Have the comments and suggestions made by mecej4 been helpful? Let us know if you still have problems with a test case that is reproducible.

Todd