Cholesky Decomposition for Positive Semidefinite Matrices (LAPACKE_dpstrf)

Paul_Buettiker — Tue, 13 Feb 2024 23:11:39 GMT

I would like to revisit an older question (see incorrect pivot index of cholesky decomposition) which, from my perspective, has not been conclusively answered or for which I may be interpreting the previous answer incorrectly:

For real, symmetric, and positive semidefinite matrices, a Cholesky decomposition can be performed using the LAPACKE_dpstrf(...) function, for example, for the matrix

m = {{30, 11, 13},
{11, 10, 13},
{13, 13, 17}}

(see code at the end of this message and incorrect pivot index of cholesky decomposition).

Even with the latest Intel compilers and the corresponding MKL (2024.0) on Ubuntu 20.04, the LAPACKE_dpstrf(...) function appears to produce peculiar values for the permutation matrix, namely piv = (1, 0, 3) in my view. In particular, piv = 0 seems strange to me.

The results computed by LAPACKE_dpstrf(...) would fit with piv = (1, 3, 2), but even MKL (2024.0) erroneously (?) tells me piv = (1, 0, 3).

What am I doing wrong?
Thank you very much!

----------------------------- Code -----------------------------

#include <iostream>
#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);

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:Cholesky Decomposition for Positive Semidefinite Matrices (LAPACKE_dpstrf)

Gennady_F_Intel — Mon, 11 Mar 2024 12:13:51 GMT

MKL 2024.0 returns the following ipiv ( 1,3,2)

see below the verbose output:

MKL_VERBOSE oneMKL 2023.0 Update 2 Product build 20230613 for Intel(R) 64 architecture Intel(R) Advanced Vector Extensions 2 (Intel(R) AVX2) enabled processors, Lnx 2.30GHz lp64 intel_thread

MKL_VERBOSE DPSTRF(U,3,0x7ffd288f1480,3,0x1f3bf00,3,0x7ffd288f13b8,0x1f4d100,0) 2.58ms CNR:OFF Dyn:1 FastMM:1 TID:0 NThr:36

LAPACKE_dpstrf, info == 0

5.47723, 2.37346, 2.00832,

11, 3.37145, 2.44208,

13, 13, 0.054153,

topic Re:Cholesky Decomposition for Positive Semidefinite Matrices (LAPACKE_dpstrf) in Intel® oneAPI Math Kernel Library

Cholesky Decomposition for Positive Semidefinite Matrices (LAPACKE_dpstrf)

Re:Cholesky Decomposition for Positive Semidefinite Matrices (LAPACKE_dpstrf)