I am attempting to analyze a program in Intel Inspector, but have experienced that the program crashes when PARDISO-functions are called.
The problem can reproduced by compiling and running one of the test examples from the oneAPI installation, e.g:
C:\Program Files (x86)\Intel\oneAPI\mkl\2023.0.0\examples\examples_core_c.zip\c\sparse_directsolvers\source\pardiso_unsym.c
Code example is in the below.
I build the program using default release settings in visual studio with oneAPI c++ (2023). Linking and include of MKL lib is archived with "Use oneMKL" set to "Parallel". Furthermore static linking of runtimes is employed (/MT).
With such settings, the example program runs without any problems when called normally from Visual Studio or from cmd.
However, when the program is run from Intel Inspector (v. 2023), it crashes.
Further debugging indicates that the crash occurs in first call of PARDISO(...) at line 129 and that it is due to an access violation within this call.
The crash is observed at least for the leak detection configuration runs (with default settings) in Inspector, but also seems to occur in other configurations. Using Inspector for other programs not calling PARDISO routines does not seem to cause any problems.
Is this a problem related to Intel Inspector, or is it perhaps resolvable by changing certain settings?
/*******************************************************************************
* Copyright 2004-2022 Intel Corporation.
*
* This software and the related documents are Intel copyrighted materials, and
* your use of them is governed by the express license under which they were
* provided to you (License). Unless the License provides otherwise, you may not
* use, modify, copy, publish, distribute, disclose or transmit this software or
* the related documents without Intel's prior written permission.
*
* This software and the related documents are provided as is, with no express
* or implied warranties, other than those that are expressly stated in the
* License.
*******************************************************************************/
/*
* Content : Intel(R) oneAPI Math Kernel Library (oneMKL) PARDISO C example
*
********************************************************************************
*/
/* -------------------------------------------------------------------- */
/* Example program to show the use of the "PARDISO" routine */
/* on unsymmetric linear systems */
/* -------------------------------------------------------------------- */
/* This program can be downloaded from the following site: */
/* www.pardiso-project.org */
/* */
/* (C) Olaf Schenk, Department of Computer Science, */
/* University of Basel, Switzerland. */
/* Email: olaf.schenk@unibas.ch */
/* -------------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mkl_pardiso.h"
#include "mkl_types.h"
#include "mkl_spblas.h"
// Define the format to printf MKL_INT values
#if !defined(MKL_ILP64)
#define IFORMAT "%i"
#else
#define IFORMAT "%lli"
#endif
int main (void)
{
/* Matrix data. */
MKL_INT n = 5;
MKL_INT ia[6] = { 1, 4, 6, 9, 12, 14};
MKL_INT ja[13] =
{ 1, 2, 4,
1, 2,
3, 4, 5,
1, 3, 4,
2, 5
};
double a[13] =
{ 1.0,-1.0, -3.0,
-2.0, 5.0,
4.0, 6.0, 4.0,
-4.0, 2.0, 7.0,
8.0, -5.0
};
MKL_INT mtype = 11; /* Real unsymmetric matrix */
// Descriptor of main sparse matrix properties
struct matrix_descr descrA;
// Structure with sparse matrix stored in CSR format
sparse_matrix_t csrA;
sparse_operation_t transA;
/* RHS and solution vectors. */
double b[5], x[5], bs[5], res, res0;
MKL_INT nrhs = 1; /* Number of right hand sides. */
/* 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, j;
double ddum; /* Double dummy */
MKL_INT idum; /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. 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[3] = 0; /* No iterative-direct algorithm */
iparm[4] = 0; /* No user fill-in reducing permutation */
iparm[5] = 0; /* Write solution into x */
iparm[6] = 0; /* Not in use */
iparm[7] = 2; /* Max numbers of iterative refinement steps */
iparm[8] = 0; /* Not in use */
iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */
iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */
iparm[11] = 0; /* Conjugate transposed/transpose solve */
iparm[12] = 1; /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
iparm[13] = 0; /* Output: Number of perturbed pivots */
iparm[14] = 0; /* Not in use */
iparm[15] = 0; /* Not in use */
iparm[16] = 0; /* Not in use */
iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
iparm[18] = -1; /* Output: Mflops for LU factorization */
iparm[19] = 0; /* Output: Numbers of CG Iterations */
maxfct = 1; /* Maximum number of numerical factorizations. */
mnum = 1; /* Which factorization to use. */
msglvl = 0; /* Print statistical information */
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;
}
/* -------------------------------------------------------------------- */
/* .. 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, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if ( error != 0 )
{
printf ("\nERROR during symbolic factorization: " IFORMAT, error);
exit (1);
}
printf ("\nReordering completed ... ");
printf ("\nNumber of nonzeros in factors = " IFORMAT, iparm[17]);
printf ("\nNumber of factorization MFLOPS = " IFORMAT, iparm[18]);
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
phase = 22;
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if ( error != 0 )
{
printf ("\nERROR during numerical factorization: " IFORMAT, error);
exit (2);
}
printf ("\nFactorization completed ... ");
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
phase = 33;
descrA.type = SPARSE_MATRIX_TYPE_GENERAL;
descrA.mode = SPARSE_FILL_MODE_UPPER;
descrA.diag = SPARSE_DIAG_NON_UNIT;
mkl_sparse_d_create_csr ( &csrA, SPARSE_INDEX_BASE_ONE, n, n, ia, ia+1, ja, a );
/* Set right hand side to one. */
for ( i = 0; i < n; i++ )
{
b[i] = 1;
}
// Loop over 3 solving steps: Ax=b, AHx=b and ATx=b
for ( i = 0; i < 3; i++ )
{
iparm[11] = i; /* Conjugate transposed/transpose solve */
if (i == 0)
transA = SPARSE_OPERATION_NON_TRANSPOSE;
else if (i == 1)
transA = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;
else
transA = SPARSE_OPERATION_TRANSPOSE;
printf ("\n\nSolving system with iparm[11] = " IFORMAT " ...\n", iparm[11]);
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
if ( error != 0 )
{
printf ("\nERROR during solution: " IFORMAT, error);
exit (3);
}
printf ("\nThe solution of the system is: ");
for ( j = 0; j < n; j++ )
{
printf ("\n x [" IFORMAT "] = % f", j, x[j]);
}
printf ("\n");
// Compute residual
mkl_sparse_d_mv( transA, 1.0, csrA, descrA, x, 0.0, bs);
res = 0.0;
res0 = 0.0;
for ( j = 1; j <= n; j++ )
{
res += (bs[j - 1] - b[j - 1]) * (bs[j - 1] - b[j - 1]);
res0 += b[j - 1] * b[j - 1];
}
res = sqrt (res) / sqrt (res0);
printf ("\nRelative residual = %e", res);
// Check residual
if ( res > 1e-10 )
{
printf ("Error: residual is too high!\n");
exit (10 + i);
}
}
mkl_sparse_destroy(csrA);
/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
phase = -1; /* Release internal memory. */
PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
&n, &ddum, ia, ja, &idum, &nrhs,
iparm, &msglvl, &ddum, &ddum, &error);
return 0;
}
