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