PARDISO for complex symmetric matrix

[Zhanghong_T_]

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

[Zhanghong_T_]

Dear all,

Could anyone help me to take a look at this part of code, or could anyone provide an example code to solve complex symmetric problem? The input matrix is coordinate format, only the upper triangular part is inputted.

Thanks