Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
The Intel sign-in experience has changed to support enhanced security controls. If you sign in, click here for more information.
6743 Discussions

incorrect pivot index of cholesky decomposition

weishenme
Beginner
367 Views

Can anyone confirm the output of LAPACKE_dpstrf on matrix

a = [ 30, 11, 13,
        11, 10, 13,
       13, 13, 17]?

I got U = [5.47723,2.37346,2.00832

                0,3.37145,2.44208

                 0,0,0.054153]

 

and U'*U =[30,13,11

                   13,17,13

                    11,13,10].

 

Therefore, the permutation matrix must be equal to (switching row/col 2 and 3)

P = [ 1 0 0
0 0 1
0 1 0]

However, the piv I got is [1,0,3], which makes no sense. So either I got the wrong piv or I made wrong interpretation of [1,0,3].

 

Here is the code

#include <stdlib.h>
#include <stdio.h>
#include "mkl_lapacke.h"
#include "mkl.h"

 

int main() {

double a[3*3] = {
30, 11, 13,
11, 10, 13,
13, 13, 17
};

int m=3;

lapack_int* piv = (lapack_int*)mkl_malloc(m*sizeof(lapack_int), 64);
long tmp = (long)piv;

int rank;

int info = LAPACKE_dpstrf(LAPACK_COL_MAJOR, 'U', m, a, m, piv, &rank, -1);

for(int i=1;i<=m;i++){
for(int j=1;j<=m;j++)
std::cout << a[(i-1)+(j-1)*m] << ", ";
std::cout << "\n";
}

for(int i=1;i<=m;i++)
std::cout << piv[i-1] << "\n";

mkl_free(piv);
}

 

 

0 Kudos
2 Replies
Gennady_F_Intel
Moderator
357 Views

Checking with the MKL 2020 u2, i see the ipiv == 1,3,2


MKL_VERBOSE Intel(R) MKL 2020.0 Update 2 Product build 20200624 for Intel(R) 64 architecture Intel(R) Advanced Vector Extensions 2 (Intel(R) AVX2) enabled processors, Win 2.60GHz cdecl intel_thread

MKL_VERBOSE DPSTRF(U,3,0000002BD3B4FA90,3,000001DC6F54E640,3,0000002BD3B4F9E8,000001DC6F56F180,0) 100.79us CNR:OFF Dyn:1 FastMM:1 TID:0 NThr:2

...LAPACKE_dpstrf results == 0

5.47723, 2.37346, 2.00832,

11, 3.37145, 2.44208,

13, 13, 0.054153,

1

3

2




Gennady_F_Intel
Moderator
323 Views

The issue is closing and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.


Reply