lapacke_dpptrf could contain a bug

[igino_p_]

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

[mecej4]

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.

[Alexander_K_Intel3]

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

[Gennady_F_Intel]

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.