Odd behavior in Pardiso Back Substitution with MKL 2025

[SChakraborty]

Hello,

I have a strange problem with MKL 2025  Pardiso Solve. 

I have simple code (attached) which is essentially loading a sparse matrix,  running a two-step factorization and then running back substitution for an input vector.  In short, the code looks like

A.load(File)  //  Square sparse matrix in Yale format

B =  Simple vector  internally crated // Size == A.rowSize()

X = 0  // Solution initialized to 0

// Phase 11

A.symbolic_factorization(....)

// Phase 22 

A.numerical_factorization(....) 

//Phase 33 

A.BackSubstitute(B, X)  //  Expecting   X to have the solution of AX = B 

I am attaching the matrix ( named Good and Bad).  

Using MKL 2021.3, for both Good and Bad matrix I am getting  AX = B 

However

Using MKL 2025.0 ,  I am getting the expected results for the Good case,  but it is not working for the Bad case.  In fact strangely,  I am getting X ==0  after the backsubstibution call  as if A is behaving like an identity matrix.  I am not getting any errors in Phases 11, 22 or 33.  

Clearly it is a case specific problem  as the Good matrix is showing ok ( also we have tested this for many different matrices and everything works fine except this one case).  Any help will be appreciated.

I have included the C++ file,  along with the input matrices ( Yale*.txt)  and the B's and X's for the good and bad cases. 

Thanks

Swagato

PS:  For some reason it is not attaching my C++ code even if I am calling it a text file.  So  pasting the code here ( apologies for the indentation getting messed up) :

#include <stdlib.h>

#include <stddef.h>

#include <iostream>

#include <string>

#include <fstream>

#include <mkl.h>

#include <mkl_pardiso.h>

#include <mkl_spblas.h>

#include <vector>

typedef struct { double r, i; } doublecomplex;

void loadIaJaAFile(const char* fileName, int& n, int& nnz, int** ia, int** ja, doublecomplex** a)

{

std::ifstream myfile;

myfile.open(fileName);

std::string text;

char* pch;

if (myfile.is_open())

{

// read number of equations

getline(myfile, text);

n = atoi(text.c_str());

// read number of nonzeros

getline(myfile, text);

nnz = atoi(text.c_str());

// allocate memory for csr matrix

*ia = (int*)malloc((n + 1) * sizeof(int));

*ja = (int*)malloc(nnz * sizeof(int));

*a = (doublecomplex*)malloc(nnz * sizeof(doublecomplex));

// read and convert input data

for (int i = 0; i <= n; i++)

{

getline(myfile, text);

(*ia)[i] = atoi(text.c_str()) - 1;

}

for (int i = 0; i < nnz; i++)

{

getline(myfile, text);

(*ja)[i] = atoi(text.c_str()) - 1;

}

for (int i = 0; i < nnz; i++)

{

getline(myfile, text);

pch = strtok((char*)text.c_str(), " ,");

if (pch[0] == '(')

pch[0] = ' ';

(*a)[i].r = atof(pch);

pch = strtok(NULL, " ");

if (pch[(strlen(pch) - 1)] == ')')

pch[(strlen(pch) - 1)] = ' ';

(*a)[i].i = atof(pch);

}

}

else

{

std::cout << "Unable to open file" << fileName;

exit(0);

}

std::cout << "# of equations " << n << ", # of  nonzeros " << (*ia)[n] - 1 << std::endl;

}

void freeMem(int** ia, int** ja, doublecomplex** a, doublecomplex** b, doublecomplex** x)

{

free(*ia);

free(*ja);

free(*a);

free(*b);

free(*x);

}

int testMKLPardisoDirect(const char* iajaaFile, const char* inputFile, const char* outputFile)

{

int    n, nnz, nrhs = 1;

int* ia, * ja;

doublecomplex* a;

std::cout << "Loading matrix data..." << std::endl;

loadIaJaAFile(iajaaFile, n, nnz, &ia, &ja, &a);

//loadIaJaAFile("small.iajaa", n, nnz, &ia, &ja, &a);

doublecomplex* b = (doublecomplex*)malloc(n * sizeof(doublecomplex));

doublecomplex* x = (doublecomplex*)malloc(n * sizeof(doublecomplex));

std::cout << "Setting up RHS..." << std::endl;

for (int i = 0; i < n * nrhs; i++)

{

b[i].r = (double)i / (double)n;  b[i].i = (double)i / (double)n;

x[i].r = 0.0;  x[i].i = 0.0;

}

std::cout << "Finished setting up RHS..." << std::endl;

mkl_set_num_threads(1);

/* Matrix data. */

MKL_INT mtype = 13;       /* Complex symmetric matrix */

/* RHS and solution vectors. */

/* Internal solver memory pointer pt, */

/* 32-bit: int pt[64]; 64-bit: long int pt[64] */

/* or void *pt[64] should be OK on both architectures */

void* pt[64];

/* Pardiso control parameters. */

MKL_INT iparm[64];

MKL_INT maxfct, mnum, phase, error, msglvl;

/* Auxiliary variables. */

MKL_INT i;

double ddum;          /* Double dummy */

std::vector<MKL_INT> irowSel;

/* -------------------------------------*/

/* .. Setup Pardiso control parameters. */

/* -------------------------------------*/

for (i = 0; i < 64; i++)

iparm[i] = 0;

iparm[0] = 1;         /* No solver default */

iparm[1] = 2;         /* Fill-in reordering from METIS */

iparm[5] = 1; //New

iparm[7] = 0;// 2;         /* Max numbers of iterative refinement steps */

iparm[9] = 13;        /* Perturb the pivot elements with 1E-13 */

iparm[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */

iparm[12] = 1;//

iparm[17] = -1;       /* Output: Number of nonzeros in the factor LU */

iparm[18] = -1;       /* Output: Mflops for LU factorization */

iparm[24] = 1; // 

iparm[34] = 1;        /* PARDISO use C-style indexing for ia and ja arrays */

maxfct = 1;           /* Maximum number of numerical factorizations. */

mnum = 1;         /* Which factorization to use. */

msglvl = 1;           /* Print statistical information in file */

error = 0;            /* Initialize error flag */

/* ----------------------------------------------------------------*/

/* .. Initialize the internal solver memory pointer. This is only  */

/*   necessary for the FIRST call of the PARDISO solver.           */

/* ----------------------------------------------------------------*/

for (i = 0; i < 64; i++)

pt[i] = 0;

std::vector<int> rowSelVec;

// Selective output rows

irowSel.resize(1, 0);

/* --------------------------------------------------------------------*/

/* .. Reordering and Symbolic Factorization. This step also allocates  */

/*    all memory that is necessary for the factorization.              */

/* --------------------------------------------------------------------*/

phase = 11;

PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &irowSel[0], &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

if (error != 0)

{

printf("\nERROR during symbolic factorization: %d", error);

freeMem(&ia, &ja, &a, &b, &x);

exit(1);

}

printf("\nReordering completed ... ");

printf("\nNumber of nonzeros in factors = %d", iparm[17]);

printf("\nNumber of factorization MFLOPS = %d", iparm[18]);

/* ----------------------------*/

/* .. Numerical factorization. */

/* ----------------------------*/

phase = 22;

PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &irowSel[0], &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

if (error != 0)

{

printf("\nERROR during numerical factorization: %d", error);

freeMem(&ia, &ja, &a, &b, &x);

exit(2);

}

printf("\nFactorization completed ... ");

/* -----------------------------------------------*/

/* .. Back substitution and iterative refinement. */

/* -----------------------------------------------*/

phase = 33;

iparm[5] = 0;   /* Write solution into 0: x, 1: b*/

PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &irowSel[0], &nrhs, iparm, &msglvl, b, x, &error);

if (error != 0)

{

printf("\nERROR during solution: %d", error);

freeMem(&ia, &ja, &a, &b, &x);

exit(3);

}

printf("\nSolve completed ... ");

std::ofstream inFile(inputFile);

for (i = 0; i < n; i++)

inFile << b[i].r << " " << b[i].i << std::endl;

inFile.close();

std::ofstream outFile(outputFile);

for (i = 0; i < n; i++)

outFile << x[i].r << " " << x[i].i << std::endl;

outFile.close();

/* --------------------------------------*/

/* .. Termination and release of memory. */

/* --------------------------------------*/

phase = -1;           /* Release internal memory. */

PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, &ddum, ia, ja, &irowSel[0], &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

freeMem(&ia, &ja, &a, &b, &x);

return 0;

}

int main(void)

{

std::string goodMatrix = "Yale_Good_490.txt";

std::string goodB = "Good_B.txt";

std::string goodX = "Good_X.txt";

testMKLPardisoDirect(goodMatrix.c_str(), goodB.c_str(), goodX.c_str());

std::string badMatrix = "Yale_Bad_5289.txt";

std::string badB = "Bad_B.txt";

std::string badX = "Bad_X.txt";

testMKLPardisoDirect(badMatrix.c_str(), badB.c_str(), badX.c_str());

return 0;

}

[Fengrui]

Hi,

Thank you for posting in the forum. All-zero output for some matrices with iparm[7]=0 is a known issue for oneMKL 2025.0. It has been fixed and the fix will be included in the upcoming 2025.1 release. 

I tried the Bad_matrix case with the fixed version and the results are expected.

Btw, starting from oneMKL 2025.0 iparm[7]=0 will turn off the iterative refinement completely, while in prior versions iparm[7]=0 performs 2 iterative refinement steps.

Thanks,

Fengrui

[SChakraborty]

@Fengrui :  Thanks for the quick response.   What is the ETA of 2025.1  

To make the 2025.0 behave similar to previous versions, should we change iparm[7] = 2?  

[Fengrui]

It is planned for late March timeframe. I will update here when it is available.

Yes, setting iparm[7]=2 will bring the default behavior (2 iterative refinement steps) of previous versions back. Please give it a try. I just tried on my side, and the results seem to be expected (not trivial solutions).

[SChakraborty]

Hello,

  Is there a list of other known issues with MKL 2025.0  Pardiso ?   

  I ran an experiment with iparm[7]=2 and in my isolated matrix solve environment, that is indeed producing similar ( not identical) results to what I get with  2021.  However,  when I am running our entire solution ( 3D EM Solver, with many sparse matrix operations)  I am getting bad result with the 2025 version with ipard[7]=2 and that gets fixed with 2021 version.  

Before we dive deeper to isolate what is going wrong at a single matrix operation level,  wanted to check if there are other known issues that might be impacting us. 

Thanks

[c_sim]

Hi,

Thank you for posting the issue and testing iparm[7]=2 setting.

An associated problem might happen if iparm[20] is set to 2 or 3 for symmetric indefinite matrices (mtype=-2,-4 or 6). Do you set them?

We are not aware of any other new PARDISO issue in oneMKL 2025.0.

Could you maybe try to isolate the matrix and share it with us.

Thank you,

Chris