Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- problem with the C interface to Lapack LU factorization

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

Fabio_G_

Beginner

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

05-06-2015
11:34 AM

58 Views

problem with the C interface to Lapack LU factorization

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

4 Replies

Fabio_G_

Beginner

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

05-06-2015
11:53 AM

58 Views

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

Black Belt

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

05-06-2015
02:33 PM

58 Views

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_

Beginner

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

05-06-2015
10:38 PM

58 Views

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

Black Belt

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

05-07-2015
05:03 AM

58 Views

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

For more complete information about compiler optimizations, see our Optimization Notice.