Why the function “LAPACKE_spbtrf” of Intel MKL library keeps returning 1?

[charlesqhappy]

guys, I'm using the function "LAPACKE_spbtrs" of MKL library to solve symmetric-positive-definite-band matrix equations. The 1st step is factorizing the matrix by "LAPACKE_spbtrf" according to the office website: https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-c/top/lapack-routines/lapack-linear-equation-routines/lapack-linear-equation-computational-routines/solving-systems-of-linear-equations-lapack-computational-routines/pbtrs.html#pbtrs.

But when I run the "LAPACKE_spbtrf", it always return 1, which indicates that the input matrix is not positive-definite. But I'm pretty sure it is.

Here is the full code:

/* Parameters */

#define N 5

#define NRHS 1

#define LDA N

#define LDB NRHS

/* Main program */ int main() {

//convert a symmetric band matrix A to lower band storage mode LA
/*A =  1 1
       1 3
*/

float ab[4] = { 0,1,
                1,3 };
const int n = 2;
int kd = 1;
int lda = kd + 1;
int info = LAPACKE_spbtrf(LAPACK_ROW_MAJOR, 'L', n, kd, ab, lda);
std::cout << "info:" << info << std::endl;

}

What am I missing?

[charlesqhappy]

Intersting, when I set the lower mode  "L" to upper mode "U", the function will return 0, which indicate sucess. But I suppose I do store the lower part of the original matrix, right?

[mecej4 · Accepted solution]

If you wanted to provide the lower triangle of the symmetric matrix, following the conventions at Netlib , you should set the values into the matrix using

float ab[4] = { 1,3,
                1,0 };

That is, the main diagonal should be in the first row, followed by the sub-diagonal in the second row. The zero is just a placeholder, since the subdiagonal has only one element, and that zero corresponds to the '*' in the diagram at the link that I listed.

[charlesqhappy]

Hey Mecej, thanks very very much, it works! This problem has annoyed me 3 days. But I still have a confusion. It is that why I should use the band storage mode of Netlib instead of the MKL's?  After all, I'm calling the function of MKL. And the MKl website gives a different storage mode of the band matrix, of which the diagonal elements are in a column :(I use row major layout as you can see in my funcion "LAPACKE_spbtrf")

But in Netlib, they are in a row:

[mecej4]

Let me answer in a broad way since, if I try to analyze what went wrong with your original attempt, it is easy to get confused by C versus Fortran, row versus column versus band storage, and the version of MKL that you use versus the version of the MKL documentation that I refer to, etc.

Lapack is open source, and is maintained and contributed to by a large number of people from around the world. That fact can tend to increase inconsistencies between the software and the documentation, but there is a counter-influence: people worldwide use Lapack, and discrepancies get reported and fixed quickly.

MKL is different. It is closed source, and belongs to Intel. Lapack is a small subset of what MKL covers, and the MKL documentation has to cover much more than Lapack. Secondly, C users will probably use the LapackE interface to Lapack, and LapackE is somewhat of a step-child.

In the specific case of SPBTRF, the MKL  documentation and the Netlib documentation differ regarding the argument AB. In such cases, my inclination is to follow the Netlib documentation and see if that works. You also made LDAB = KD+1, which makes the 1-D array that you passed for AB completely equivalent to the packed 2-D array that should have been used.

During the last few years, it did not help that at times the MKL documentation for C and Fortran usage was merged into a single document. The result was, at least for the BLAS and Lapack parts, that usually half the monolingual (i.e., C/Fortran) readers were left confused. This problem persists in some places in the current MKL documentation, even though the C and Fortran versions have been "separated". One could say, "they have taken the C out of the Fortran, but they failed to take the Fortran out of the C". There are parts of MKL where the opposite is true.

[charlesqhappy]

Hi , Mecej, thanks! But I got the same problem with another matrix using the Netlib band matrix storage mode. Please help me with it. The link is :

https://community.intel.com/t5/Intel-oneAPI-Math-Kernel-Library/Why-the-quot-LAPACKE-spbtrf-quot-always-return-quot-6-quot/m-p/1212720#M30115