Community
cancel
Showing results for
Did you mean:
Beginner
72 Views

## PARDISO for complex symmetric matrix

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

2 Replies
Employee
72 Views

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

Beginner
72 Views

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