topic Re:incorrect pivot index of cholesky decomposition in Intel® oneAPI Math Kernel Library

incorrect pivot index of cholesky decomposition

weishenme — Wed, 09 Sep 2020 12:26:26 GMT

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);
}

Re:incorrect pivot index of cholesky decomposition

Gennady_F_Intel — Thu, 10 Sep 2020 10:42:40 GMT

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,

Re:incorrect pivot index of cholesky decomposition

Gennady_F_Intel — Thu, 01 Oct 2020 03:20:03 GMT

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.