Strange partial pivoting of LAPACKE_dgetrf

[yang__xiaolin]

Hello,

I'm using LAPACKE_dgetrf to compute the LU factorization of square matrices in double precision. The matrix is in column major. Here is what I am doing. The environment is MKL 2018 Update 3 for Windows + Visual studio 2017.

for(...)

{

lapack_int m = dim;    //dim is around 40 to 80

double * A = (double*)mkl_malloc(m * m * sizeof(double), 64);

memcpy(A, source, m * m * sizeof(double));

lapack_int * ipiv = (lapack_int*)mkl_malloc(m * sizeof(lapack_int), 64);

lapack_int info = 0;

int mat_layout = LAPACK_COL_MAJOR;

for(int k = 0; k < m; k++)    //flush ipiv before calling dgetrf

    ipiv = -1;

info = LAPACKE_dgetrf(mat_layout, m, m, A, m, ipiv);

if(info != 0)

{clean up memory;

    break;}

}

I checked the info that was returned by LAPACKE_dgetrf and it was always 0. However, I found out I got duplicate items in ipiv array after each dgetrf call. 

I mean I got ipiv == ipiv  (0 < i , j < m), which doesn't make sense for partial pivoting. I also flushed the ipiv array before calling dgetrf.

What is the possible reason for getting duplicate values in partial pivoting array?

Thank you.

[Ying_H_Intel]

Hi Xiaolin,

could you please attach one test case to show the problem?

​for example, below

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <mkl.h>

int main(void)
{

 int len = 198;
 char buf[198];
 mkl_get_version_string(buf, len);
 printf("%s\n", buf);
 printf("\n");

 lapack_int m = 3; //dim is around 40 to 80
 double  A[9] = { 0, 5, 5, 2, 9, 0,6, 8, 8 };
 //double * A = (double*)mkl_malloc(m * m * sizeof(double), 64);
 //memcpy(A, source, m * m * sizeof(double));
 lapack_int * ipiv = (lapack_int*)mkl_malloc(m * sizeof(lapack_int), 64);
 lapack_int info = 0;
 int mat_layout = LAPACK_COL_MAJOR;
 for (int k = 0; k < m; k++)  //flush ipiv before calling dgetrf
  ipiv = -1;
 info = LAPACKE_dgetrf(mat_layout, m, m, A, m, ipiv);

 if (info != 0)
 {
  printf("info error %d", info);
  //break;
 }

printf("\nLU_A\n");
for (int i = 0; i <m; i++) {
 for (int j = 0; j < m; j++) {
  printf("%f ", A[i + j * 3]);
 }
 printf("\n");
}

printf("ipiv:\n");
for (int i = 0; i < m; i++) {
 printf("%d ", ipiv);
}


return 0;

}

i did a quick test with MKL 2018 u3.  the result looks fine with small size input 

Intel(R) Math Kernel Library Version 2018.0.3 Product Build 20180406 for Intel(R) 64 architecture applications

LU_A
5.000000 9.000000 8.000000
1.000000 -9.000000 0.000000
0.000000 -0.222222 6.000000
ipiv:
2 3 3 Press any key to continue . . .

Best Regards,
Ying
https://software.intel.com/en-us/articles/intel-math-kernel-library-intel-mkl-2018-bug-fixes-list

[Kirana_B_Intel]

The ipiv array is an array that lists row interchanges, it is not a true pivot vector. So if you see an ipiv array for a 3x3 case that looks like [2,2,2], that means that row 1 was interchanged with row 2, then row 2 (which is the original row 1) stayed the same, then row 3 was interchanged with row 2 (which was the original row 1). 

You get find more information about this in the description for ipiv in the getrf documentation https://software.intel.com/en-us/mkl-developer-reference-c-getrf. The routines getrs and getri are provided to work with pivot arrays of this form.

[yang__xiaolin]

Thank you.