- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page