PARDISO for complex symmetric matrix

Zhanghong_T_ · ‎05-24-2015

Dear all,

I used PARDISO for complex symmetric matrix. The code is as follows:

void PardisoSolver(int n, int nz, int *rowind, int *colind, MKL_Complex16 *a, int nrhs, MKL_Complex16 *x, MKL_Complex16 *b)
{
	//// 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.	
	int *ia = new int[n + 1], *ja = new int[nz];
	MKL_Complex16 *acsr = new MKL_Complex16[nz];
	int info;
	mkl_zcsrcoo(job, &n, acsr, ja, ia, &nz, a, rowind, colind, &info);


  MKL_INT mtype = 6;       /* 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. */

  msglvl = 1;           /* 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 = 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)
    {
      printf ("\nERROR during symbolic factorization: %d", error);
      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, acsr, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);

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

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

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

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

/* -------------------------------------------------------------------- */
/* ..  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);
	delete[] ia;
	delete[] ja;
	delete[] acsr;
}

However, the output shows that number of elements in U matrix is 1, and the res is very large.

Then I output the full matrix (both upper and lower triangular part) and set mtype=13 and solve the matrix again. The result is correct. Could anyone help me to take a look at the problem?

Thanks,

Zhanghong Tang

Alexander_K_Intel2 · ‎06-01-2015

Hi,

Can I ask you to send us example of matrix on which you see the issue? It could help us to reproduce and fix it.

Thanks,

Alex

Zhanghong_T_ · ‎06-04-2015

Dear Alex,

A complete example is given here:

https://software.intel.com/en-us/forums/topic/558881

I have also reported this issue and your developer has reproduced this problem.

Thanks,

Zhanghong Tang