Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Bug report: Pardiso get wrong result when solving complex symmetric matrix

Zhanghong_T_
Novice
714 Views

Dear all,

I use the following code to solve complex symmetric matrix by Pardiso and get a wrong result, could anyone help me to take a look at it?

 

Thanks

 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "mkl_pardiso.h"
#include "mkl_types.h"
#include "mkl_spblas.h"
#include "mkl.h"
double residual_csr(int sym, int n, int *ia, int *ja, MKL_Complex16 *a, MKL_Complex16 *x, MKL_Complex16 *b);
void PardisoSolver(int sym, int n, int coocsr, int nz, int *rowind, int *colind, MKL_Complex16 *a, int nrhs, MKL_Complex16 *x, MKL_Complex16 *b);

MKL_INT main( void ) {
	/* Matrix data. */
  MKL_INT n = 8;
	MKL_Complex16 *x = new MKL_Complex16, *b = new MKL_Complex16;
	for (int i = 0; i < n; i++)
    {
        b.real = 1.0;
        b.imag = 0.0;
    }
	MKL_INT row[16] = { 1, 1, 1, 1, 2,2,2,3,3,4,4,6,6,7,7,8};
	MKL_INT col[16] =
  { 1,    3,       6, 7,
       2, 3,    5,
          3,             8,
             4,       7,
       
                    6,    8,
                       7, 8,
                          8
  };

	MKL_Complex16 a[16];
	  a[0].real = 7.0; a[0].imag = 1.0; a[1].real = 1.0;  a[1].imag = 1.0; a[2].real = 2.0; a[2].imag = 1.0; a[3].real = 7.0; a[3].imag = 1.0;
                                      a[4].real = -4.0; a[4].imag = 0.0; a[5].real = 8.0; a[5].imag = 1.0; a[6].real = 2.0; a[6].imag = 1.0;
                                                                                          a[7].real = 7.0; a[7].imag = 0.0; a[8].real = 9.0; a[8].imag = 1.0;
                                                                                                a[9].real = -4.0; a[9].imag = 1.0; a[10].real = 7.0; a[10].imag = 1.0;

		                                                                                                  a[11].real = 1.0; a[11].imag = 1.0;  a[12].real = 11.0; a[12].imag = 1.0;
		                                                                                                      a[13].real = 5.0; a[13].imag = 0.0; a[14].real = 5.0; a[14].imag = 1.0;
		                                                                                                                                          a[15].real = 50.0; a[15].imag = 50.0;
	int nz = 16;
	PardisoSolver(0, n, 0, nz, row, col, a, 1, x, b);
	return 0;
}

void PardisoSolver(int sym, int n, int coocsr, int nz, int *rowind, int *colind, MKL_Complex16 *a, int nrhs, MKL_Complex16 *x, MKL_Complex16 *b)
{
	int *ia, *ja;
	MKL_Complex16 *acsr;
	if (coocsr == 0)
	{
		//// coordinate format to CSR format
		int job[6];
		job[0] = 2; /// if job(1)=2, the matrix in the coordinate format is converted to the CSR format, 
		/// and the column indices in CSR representation are sorted in the increasing order within each row.
		job[1] = 1; /// if job(2)=1, one-based indexing for the matrix in CSR format is used.
		job[2] = 1; /// if job(3)=1, one-based indexing for the matrix in coordinate format is used.
		job[3] = 0;
		job[4] = nz; /// job(5) = nzmax - maximum number of the non - zero elements allowed if job(1) = 0.
		job[5] = 0;  /// If job(6) = 0, all arrays acsr, ja, ia are filled in for the output storage.	
		ia = new int[n + 1];
		ja = new int[nz];
		acsr = new MKL_Complex16[nz];
		int info;
		mkl_zcsrcoo(job, &n, acsr, ja, ia, &nz, a, rowind, colind, &info);
	}
	else
	{
		ia = rowind;
		ja = colind;
		acsr = a;
	}
	MKL_INT mtype;
	if (sym == 1)
		mtype = 6;       /* complex symmetric matrix */
	else
		mtype = 13;       /* complex symmetric matrix */

	/* RHS and solution vectors. */
	double res;

	/* 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;

	MKL_Complex16 ddum;       /* Double dummy */
	MKL_INT idum;         /* Integer dummy. */
	/* -------------------------------------------------------------------- */
	/* .. Setup Pardiso control parameters. */
	/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++)
	{
		iparm = 0;
	}
	iparm[0] = 1;         /* No solver default */
	//iparm[1] = 2;         /* Fill-in reordering from METIS */
	iparm[1] = 3;         ///!parallel(OpenMP) version of the nested dissection algorithm
	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 */
	iparm[20] = 1;        /// !Apply 1x1 and 2x2 Bunch and Kaufman pivoting during the factorization process
	iparm[23] = 1;        /// !PARDISO uses new two - level factorization algorithm

	maxfct = 1;           /* Maximum number of numerical factorizations.  */
	mnum = 1;         /* Which factorization to use. */
#ifdef _DEBUG
	msglvl = 1;           /* Print statistical information  */
#else
	msglvl = 0;
#endif
	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 = 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, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

	if (error != 0)
	{
#ifdef _DEBUG
		printf("\nERROR during symbolic factorization: %d", error);
#endif
		exit(1);
	}
#ifdef _DEBUG
	printf("\nReordering completed ... ");
	printf("\nNumber of nonzeros in factors  = %d", iparm[17]);
	printf("\nNumber of factorization MFLOPS = %d", iparm[18]);
#endif

	/* -------------------------------------------------------------------- */
	/* ..  Numerical factorization.                                         */
	/* -------------------------------------------------------------------- */
	phase = 22;

	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

	if (error != 0)
	{
#ifdef _DEBUG
		printf("\nERROR during numerical factorization: %d", error);
#endif
		exit(2);
	}
#ifdef _DEBUG
	printf("\nFactorization completed ...\n ");
#endif

	/* -------------------------------------------------------------------- */
	/* ..  Back substitution and iterative refinement.                      */
	/* -------------------------------------------------------------------- */
	phase = 33;

#ifdef _DEBUG
	printf("\n\nSolving system...\n");
#endif
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);

	if (error != 0)
	{
#ifdef _DEBUG
		printf("\nERROR during solution: %d", error);
#endif
		exit(3);
	}
	// Compute residual
#ifdef _DEBUG
	res = residual_csr(sym, n, ia, ja, acsr, x, b);
	printf("\nRelative residual = %e", res);
#endif

	/* -------------------------------------------------------------------- */
	/* ..  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);
	if (coocsr == 0)
	{
		delete[] ia;
		delete[] ja;
		delete[] acsr;
	}
}

double residual_csr(int sym, int n, int *ia, int *ja, MKL_Complex16 *a, MKL_Complex16 *x, MKL_Complex16 *b)
{
	MKL_Complex16 *r = new MKL_Complex16;
	if (sym == 1)
		mkl_zcsrsymv("u", &n, a, ia, ja, x, r);
	else
		mkl_zcsrgemv("n", &n, a, ia, ja, x, r);
	double res = 0.0;
	double res0 = 0.0;
	for (int i = 0; i < n; i++)
	{
		r.real = r.real - b.real;
		r.imag = r.imag - b.imag;
	}
	for (int i = 0; i < n; i++)
	{
		res += r.real*r.real + r.imag*r.imag;
		res0 += b.real*b.real + b.imag*b.imag;
	}
	res = cblas_dznrm2(n, r, 1);
	res0 = cblas_dznrm2(n, b, 1);
	delete[]r;
	return sqrt(res) / sqrt(res0);
}

 

0 Kudos
8 Replies
Zhanghong_T_
Novice
714 Views

Dear all,

I have updated the code to verify the bug: The program read a matrix from the attached file, it is a symmetric matrix and only upper triangular part is inputted. Then the PARDISO is called with mtype=6. The residual error is very large. And then, the lower triangular part is filled and the PARDISO is called again with mtype=13. At this time the residual error is acceptable.

Could anyone help me to take a look at it?

Thanks

Zhanghong Tang

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "mkl_pardiso.h"
#include "mkl_types.h"
#include "mkl_spblas.h"
#include "mkl.h"
double residual_csr(int sym, int n, int *ia, int *ja, MKL_Complex16 *a, MKL_Complex16 *x, MKL_Complex16 *b);
void PardisoSolver(int sym, int n, int coocsr, int nz, int *rowind, int *colind, MKL_Complex16 *a, int nrhs, MKL_Complex16 *x, MKL_Complex16 *b);

MKL_INT main( void ) {
	/* Matrix data. */
  MKL_INT n;
	MKL_Complex16 *x, *b;
	int nz;
	FILE *fp = fopen("mat.txt", "rt");
	fscanf(fp, "%d %d\n", &n, &nz);
	x = new MKL_Complex16;
	b = new MKL_Complex16;
	MKL_INT *row = new MKL_INT[nz];
	MKL_INT *col = new MKL_INT[nz];
	MKL_Complex16 *a = new MKL_Complex16[nz];
	for (int i = 0; i < n; i++)
	{
		b.real = 1.0;
		b.imag = 0.0;
	}
	for (int i = 0; i < nz; i++)
	{
		fscanf(fp, "%d %d %Lf %Lf\n", &row, &col, &a.real, &a.imag);
	}
	fclose(fp);
	PardisoSolver(1, n, 0, nz, row, col, a, 1, x, b);

	/// test unsymmetric matrix
	MKL_INT *row2 = new MKL_INT[2*nz-n];
	MKL_INT *col2 = new MKL_INT[2*nz-n];
	MKL_Complex16 *a2 = new MKL_Complex16[2*nz-n];
	int k = 0;
	for (int i = 0; i < nz; i++)
	{
		row2 = row;
		col2 = col;
		a2.real = a.real;
		a2.imag = a.imag;
		if (row != col)
		{
			k++;
			row2 = col;
			col2 = row;
			a2.real = a.real;
			a2.imag = a.imag;
		}
		k++;
	}
	PardisoSolver(0, n, 0, nz, row, col, a, 1, x, b);
	delete[]x;
	delete[]b;
	delete[]row;
	delete[]col;
	delete[]a;
	delete[]row2;
	delete[]col2;
	delete[]a2;
	return 0;
}

void PardisoSolver(int sym, int n, int coocsr, int nz, int *rowind, int *colind, MKL_Complex16 *a, int nrhs, MKL_Complex16 *x, MKL_Complex16 *b)
{
	int *ia, *ja;
	MKL_Complex16 *acsr;
	if (coocsr == 0)
	{
		//// coordinate format to CSR format
		int job[6];
		job[0] = 2; /// if job(1)=2, the matrix in the coordinate format is converted to the CSR format, 
		/// and the column indices in CSR representation are sorted in the increasing order within each row.
		job[1] = 1; /// if job(2)=1, one-based indexing for the matrix in CSR format is used.
		job[2] = 1; /// if job(3)=1, one-based indexing for the matrix in coordinate format is used.
		job[3] = 0;
		job[4] = nz; /// job(5) = nzmax - maximum number of the non - zero elements allowed if job(1) = 0.
		job[5] = 0;  /// If job(6) = 0, all arrays acsr, ja, ia are filled in for the output storage.	
		ia = new int[n + 1];
		ja = new int[nz];
		acsr = new MKL_Complex16[nz];
		int info;
		mkl_zcsrcoo(job, &n, acsr, ja, ia, &nz, a, rowind, colind, &info);
	}
	else
	{
		ia = rowind;
		ja = colind;
		acsr = a;
	}

	//for (int i = 0; i < nz; i++)printf("%d,%e,%e\n", i,acsr.real, acsr.imag);
	MKL_INT mtype;
	if (sym == 1)
		mtype = 6;       /* complex symmetric matrix */
	else
		mtype = 13;       /* complex symmetric matrix */

	/* RHS and solution vectors. */
	double res;

	/* 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;

	MKL_Complex16 ddum;       /* Double dummy */
	MKL_INT idum;         /* Integer dummy. */
	/* -------------------------------------------------------------------- */
	/* .. Setup Pardiso control parameters. */
	/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++)
	{
		iparm = 0;
	}
	iparm[0] = 1;         /* No solver default */
	//iparm[1] = 2;         /* Fill-in reordering from METIS */
	iparm[1] = 3;         ///!parallel(OpenMP) version of the nested dissection algorithm
	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 */
	iparm[20] = 1;        /// !Apply 1x1 and 2x2 Bunch and Kaufman pivoting during the factorization process
	iparm[23] = 1;        /// !PARDISO uses new two - level factorization algorithm

	maxfct = 1;           /* Maximum number of numerical factorizations.  */
	mnum = 1;         /* Which factorization to use. */
#ifdef _DEBUG
	msglvl = 1;           /* Print statistical information  */
#else
	msglvl = 0;
#endif
	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 = 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, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

	if (error != 0)
	{
#ifdef _DEBUG
		printf("\nERROR during symbolic factorization: %d", error);
#endif
		exit(1);
	}
#ifdef _DEBUG
	printf("\nReordering completed ... ");
	printf("\nNumber of nonzeros in factors  = %d", iparm[17]);
	printf("\nNumber of factorization MFLOPS = %d", iparm[18]);
#endif

	/* -------------------------------------------------------------------- */
	/* ..  Numerical factorization.                                         */
	/* -------------------------------------------------------------------- */
	phase = 22;

	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

	if (error != 0)
	{
#ifdef _DEBUG
		printf("\nERROR during numerical factorization: %d", error);
#endif
		exit(2);
	}
#ifdef _DEBUG
	printf("\nFactorization completed ...\n ");
#endif

	/* -------------------------------------------------------------------- */
	/* ..  Back substitution and iterative refinement.                      */
	/* -------------------------------------------------------------------- */
	phase = 33;

#ifdef _DEBUG
	printf("\n\nSolving system...\n");
#endif
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);

	if (error != 0)
	{
#ifdef _DEBUG
		printf("\nERROR during solution: %d", error);
#endif
		exit(3);
	}
	// Compute residual
#ifdef _DEBUG
	res = residual_csr(sym, n, ia, ja, acsr, x, b);
	printf("\nRelative residual = %e\n", res);
#endif

	/* -------------------------------------------------------------------- */
	/* ..  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);
	if (coocsr == 0)
	{
		delete[] ia;
		delete[] ja;
		delete[] acsr;
	}
}

double residual_csr(int sym, int n, int *ia, int *ja, MKL_Complex16 *a, MKL_Complex16 *x, MKL_Complex16 *b)
{
	MKL_Complex16 *r = new MKL_Complex16;
	if (sym == 1)
		mkl_zcsrsymv("u", &n, a, ia, ja, x, r);
	else
		mkl_zcsrgemv("n", &n, a, ia, ja, x, r);
	for (int i = 0; i < n; i++)
	{
		r.real = r.real - b.real;
		r.imag = r.imag - b.imag;
	}
	double res, res0;
	res = cblas_dznrm2(n, r, 1);
	res0 = cblas_dznrm2(n, b, 1);
	delete[]r;
	return res / res0;
}

 

0 Kudos
Alexander_K_Intel2
714 Views

Hi,

I've reproduce behavior that you see, look's like an issue.

Thanks,

Alex

0 Kudos
Zhanghong_T_
Novice
714 Views

Dear Alex,

Thank you very much for your kindly reply. I am waiting for this problem to be solved.

Thanks

0 Kudos
Alexander_K_Intel1
714 Views
Hi Zhanghong, Your problem heavily relies upon scaling which is better done for nonsymmetrical case: elements can be permuted without preserving matrix symmetry. Turning off scaling for nonsymmetrical case (essentially reducing problem reliability improvement abilities to symmetric case) leads to same large residual. Best regards, Alexander
0 Kudos
Zhanghong_T_
Novice
714 Views

Dear Alexander,

Thank you very much for your kindly reply.

Yes, the scaling is very important for my problem. I think it is also common for most of finite element problems when the domain is complex. Do you mean that it is an unsolved problem that the scaling is much worse if we need to preserving matrix symmetry?

Thanks

0 Kudos
Alexander_K_Intel1
714 Views
Hi Zhanghong, Only diagonal elements can be permuted in symmetric matrices: it is impossible to choose arbitrary (i.e. largest) element from row/column and put it on a diagonal without losing symmetry. So, yes, scaling is naturally limited if we work with only one triangle part of the matrix. Best regards, Alexander
0 Kudos
Zhanghong_T_
Novice
714 Views

Dear Alexander,

Thank you very much for your kindly reply. I understood now. However, I think that the PARDISO should do something to make sure it find a correct solution. The worst condition is that when it found that the permutation can't be done without losing symmetry, at that time it should expand the matrix into full one and take it as asymmetric matrix.

Thanks

0 Kudos
Alexander_K_Intel2
714 Views

Hi,

The symmetric permutation is worse than nonsymmetric for almost all case (that is clear) but symmetric permutation is most cases provide quite good solution without significant memory and computational time increase (switching to nonsymmetric mean change LDL^t decomposition on LDU). And the fact is that understand poor factorization of symmetric case one can only after factorization is over - by calculating residual for example. Practically this idea similar to your example - you implement symmetric factorization, check it's residual and switch to nonsymmetric and i don't see any chance to implement it in quicker way.

Thanks,

Alex   

0 Kudos
Reply