problem with the C interface to Lapack LU factorization

[Fabio_G_]

Hi,

I've written the following C code that should perform the LU factorization

of a general n x n matrix

#include <mkl.h>

#include "mkl_lapack.h"

#define N 3

double* A;
lapack_int n=N;
lapack_int* ipiv;

int main(){
    int i=0,j=0,ierr=0;
    A=malloc(n*n*sizeof(double));
    ipiv=malloc(n*sizeof(lapack_int));
    A[0]=1e0;A[1]=0e0;A[2]=0e0;
    A[3]=2e0;A[4]=3e0;A[5]=0e0;
    A[6]=4e0;A[7]=5e0;A[8]=6e0;
    ierr=LAPACKE_dgetrf(LAPACK_ROW_MAJOR,n,n,A,n,ipiv);
    for(i=0;i<n;i++){
        for(j=0;j<n;j++)
          printf("%f ",A[i*n+j]);
        printf("\n");
    }
    printf("ierr=%d\n",ierr);
    free(A);
    free(ipiv);
    return 0;
}

the result (computed with mathematica) should give for U the 3x3 diagonal matrix {1,3,6}

and for L={{1., 0., 0.}, {2., 1., 0.}, {4., 1.66667, 1.}}. However, the code seems to fail.

Am I doing some mistake?

Thank you!

F

[Fabio_G_]

I guess I am misunderstanding the storage of the result. I have tried

to solve a linear system with b={1,1,1} and it seems that it is working.

How to interpret the result from LU stored in A?

Thanks

[mecej4]

I think that the issue is how to interpret the results, since your program works correctly. Its output agrees with that Matlab prints when given the command lu(A).

Different packages  (e.g., MKL, Matlab, Mathematica) may use different conventions since, as with finding the factors of composite integers, the P,L,U factors of a matrix A, A = P.L.U,  are not unique. The convention used in MKL is that L has a unit diagonal.

Your program prints the factors of a permutation of A, not those of A itself, and you need to take this into account, as well, when comparing results.

Commonly, finding the L and U factors is the first of a two step solution of the problem of solving linear equations, and the second step, in which you call DGETRS, handles the permutation vector properly but transparently. Furthermore, the solution is independent of the convention used for the L-U factors.

[Fabio_G_]

I agree, indeed by using  dgetres gave me the right solution for a system. However, there is no reference in the MKL guide how to interpret this result and I think should be.

Thank you,

F

[mecej4]

Please see the MKL reference manual, https://software.intel.com/en-us/node/520877, where it says about DGETRF:

The routine computes the LU factorization of a general m-by-n matrix A as A = P*L*U, where P is a permutation matrix, L is lower triangular with unit diagonal elements (lower trapezoidal if m > n) and U is upper triangular (upper trapezoidal if m < n). The routine uses partial pivoting, with row interchanges.