Showing results for

- Intel Community
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library & Intel® Math Kernel Library
- lapacke_dpptrf could contain a bug

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

Highlighted
##

Hello,
using the intel MKL 11.1 for windows, i discovered what seems a bug in LAPACKE C interface to Lapack: calling LAPACKE_dpptrf (Cholesky decomposition) with a packed LT row major matrix does not produce the expected result. Results are correct only passing a column major matrix and thus the LAPACK dpptrf is right and i suspect the LAPACKE interface is not. Searching on the web i found the source file lapacke_dpptrf.c on the netlib site (with intel copyright) and the related lapacke_dpptrf_work.c which, i think, is the source of the bug. (see http://www.netlib.org/lapack/explore-html/de/ddd/lapacke_8h.html)
In the lapacke_dpptrf_work function a copy ap_t of the working matrix ap is transposed before and after the call of the factorization routine. The double transposition should not really necessary since a packed row major LT matrix has the same linear array representation of a packed column major UT matrix (and viceversa) and thus it would suffice to call LAPACK_dpptrf(&uplo,6n, &ap,&info) with uplo='U' if ap is a row major LT matrix and with uplo='L' if ap is a row major UT matrix.
The correct code is very simple and does not require temporaries (ap_t) nor transpositions:
char tmp_uplo=uplo;
if(matrix_order==RowMajor)
tmp_uplo=(uplo=='L')?'U':'L';
LAPACK_dpptrf( &tmp_uplo, &n, &ap_t, &info );
Note that after the factorization it isnt necessary to transpose the factor since it is already correct.
I tested the above code and the factor is the correct one.
Clearly also the companion solver routine LAPACKE_dpptrs_work should be affected by the same problem and there you make 3 transpositions! one for the working matrix and two for the rhs, while only a pre and a post rhs transposition would suffice (may be a better solution exists).
best regards

igino_p_

Beginner

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

01-18-2014
07:54 AM

15 Views

lapacke_dpptrf could contain a bug

3 Replies

Highlighted
##

mecej4

Black Belt

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

01-18-2014
10:53 AM

15 Views

A careful reading of the MKL documentation on ?potrf and packed storage shows that there is no place for row-major or column-major packed storage issues. The concepts of row- or column-major apply only to arrays of rank greater than 1.

There is only one convention for packed storage, and that is **by columns**. If one is storing a lower triangular matrix, the packed 1-D array contains n elements from col-1, then n-1 elements from col-2,..., 1 element from col-n. If storing an upper triangular matrix, the packed 1-D array contains 1 element from col-1, then 2 from col-2,..., n elements from col-n.

Highlighted
##

Alexander_K_Intel3

Employee

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

05-22-2014
11:10 PM

15 Views

Hi,

It seems there is indeed a bug in LAPACKE transposition for symmetric packed format if it's stored initially in Row major.

W.B.R., Alexander

Highlighted
##

igino, please check the latest version 11.2 update 1 where the fix of that problem has been released and let us know the results.

Gennady_F_Intel

Moderator

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

11-30-2014
08:59 PM

15 Views

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