lapacke_dpptrf could contain a bug

lapacke_dpptrf could contain a bug

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

5 帖子 / 0 全新
最新文章
如需更全面地了解编译器优化,请参阅优化注意事项

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.

 

@ mecej4

Yes, in lapack packed arrays there is no place for the ordering. But LAPACKE ( and the final E is not an error) is a C interface to LAPACK supported by intel and the first argument of LAPACKE_dpptrf is  the ordering  (RowMajor, ColMajor) of the Matrix to be factorized. If you look at the source code in lapacke_dpptrf_work.c  you will find how RowMajor matrices are treated differently from the ColMajor, the latter direcly passed to the lapack dpptrf, the former first transposed and than transposed again after dpptrf return. 

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

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.

 

登陆并发表评论。