- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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);
}
Link Copied
8 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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;
}
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
I've reproduce behavior that you see, look's like an issue.
Thanks,
Alex
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear Alex,
Thank you very much for your kindly reply. I am waiting for this problem to be solved.
Thanks
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page