Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
FPGA community forums and blogs have moved to the Altera Community. Existing Intel Community members can sign in with their current credentials.

MKL function dgetrf wrong results

dimdol10
New Contributor I
1,838 Views

Hello!

I am using Intel oneAPI MKL version 2022.2, and visual studio 2022.

I am currently trying a LU decomposition of a matrix using the function DGETRF provided by MKL. 

 

As I found a wrong result in a simple example, I would like ask you if something is incorrect. The example code is attached in the below. As far as I understood, in this example, the array ipiv should be 

ipiv(1) = 2,

ipiv(2) = 1,

such that 

A_PLU = P*L*U. 

However, the function gives

ipiv(1) = 2,

ipiv(2) = 2.

 

---------------------------------------------------------

PROGRAM TEST_MAT_LU_FAC_INV

DOUBLE PRECISION A(2,2)
DOUBLE PRECISION A_L(2,2), A_U(2,2)
DOUBLE PRECISION A_LU(2,2)
DOUBLE PRECISION A_PLU(2,2), A_PLU_TEST(2,2)
integer info
integer ipiv(2)
DOUBLE PRECISION P(2,2) ! permutation matrix

A = 0.D0
A(1,1) = 1.D0
A(1,2) = 2.D0
A(2,1) = 3.D0
A(2,2) = 4.D0

CALL DGETRF(2,2,A,2,ipiv,info)
A_L = 0.d0
A_L(1,1) = 1.D0 !A(1,1)
A_L(2,1) = A(2,1)
A_L(2,2) = 1.D0 !A(2,2)

A_U = 0.D0
A_U(1,1) = A(1,1)
A_U(1,2) = A(1,2)
A_U(2,2) = A(2,2)

A_LU = 0.D0
DO I = 1,2
DO J = 1,2
DO K = 1,2
A_LU(I,J) = A_LU(I,J) + A_L(I,K)*A_U(K,J)
END DO
END DO
END DO

P = 0.D0
DO I = 1,2
P(I,IPIV(I)) = 1.D0
END DO

A_PLU_TEST = 0.D0
DO I = 1,2
DO J = 1,2
DO K = 1,2
A_PLU_TEST(I,J) = A_PLU_TEST(I,J) + A_L(I,K)*A_U(K,J)
END DO
END DO
END DO

A_PLU = 0.D0
DO I = 1,2
DO J = 1,2
DO K = 1,2
A_PLU(I,J) = A_PLU(I,J) + P(I,K)*A_PLU_TEST(K,J)
END DO
END DO
END DO


END

0 Kudos
5 Replies
ShanmukhS_Intel
Moderator
1,805 Views

Hi Myung,


Thanks for posting on Intel Communities.


Could you please cross-verify the shared source code? We tried compiling it and we couldn't see any output.


Best Regards,

Shanmukh.SS


0 Kudos
dimdol10
New Contributor I
1,796 Views

Dear Shanmukh,

 

Thanks for the reply. 

Yes, the code print out nothing. But, what I meant is that the array A_PLU is supposed to be the same with the array A. However, if you run the code, you can find that they are different. 

 

A_PLU(1,1) = 1.D0

A_PLU(2,1) = 1.D0

A_PLU(1,2) = 2.D0

A_PLU(2,2) = 2.D0

 

I suspect the output array ipiv from DGETRF is wrong. Could you please take a look at this point?

 

Thank you very much.

 

0 Kudos
ShanmukhS_Intel
Moderator
1,759 Views

Hi Myung,


Please find the below link for more details and usage regarding the DGETRF functionality.


https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/lapack-routines/lapack-linear-equation-routines/lapack-linear-equation-computational-routines/matrix-factorization-lapack-computational-routines/getrf.html


The function describes the input parameters to be used and the proper usage of the routine.


Best Regards,

Shanmukh.SS


0 Kudos
ShanmukhS_Intel
Moderator
1,718 Views

Hi,


A gentle reminder:

We have not heard back from you for a while. Could you please let us know if any updates regarding your issue?


Best Regards,

Shanmukh.SS


0 Kudos
ShanmukhS_Intel
Moderator
1,678 Views

Hi Myung,


We assume that your issue is resolved. If you need any other information, please post a new question as this thread will no longer be monitored by Intel.


Best Regards,

Shanmukh.SS


0 Kudos
Reply