pardiso_getdiag problem

[Oleg_S_]

Hello everyone!

I am using pardiso to solve nonsymmetric system, but i don't get result that i expect. Today i tried pardiso_getdiag to understand what is going on with diagonal after factorization and i noticed that pardiso thinks that i have zero elements on diagonal which is not true. And if size of the matrix rises i get more zeros. Matrix_checkers returns MKL_SPARSE_CHECKER_SUCCESS. How can i solve this problem? I attached my matrix in csr format. Indexing starts from zero.

Thank you.

[mecej4]

Your matrix is fine, and I was able to obtain a solution of A x = b with b = 1 (all elements). I used the pardiso_unsym_f.f example code from the MKL distribution, modified to read your data files.

Perhaps if you post your code one could spot problems with the Pardiso calls and arguments.

[Oleg_S_]

Thanks for quick response. Here is my code:

	*msglvl = 0;		//	Не выводить полезную информацию
	*mtype = 11;		//	Не симметричная реальная матрица
	for (int i = 0; i < 64; i++)
	{
		iparm = 0;
	}
	iparm[0] = 1;         /* No solver default */
	iparm[1] = 0;         /* Fill-in reordering from METIS */
	iparm[2] = 2;
	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] = 0;         /* Max numbers of iterative refinement steps */
	iparm[8] = 0;         /* Not in use */
	iparm[9] = 13;        /* Perturb the pivot elements with 1E-13 */
	iparm[10] = 0;        /* 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[26] = 1;
	iparm[34] = 1;
	iparm[55] = 1;

	MKL_INT nrhs = 1;     /* Number of right hand sides. */
	MKL_INT maxfct, mnum, phase, error;
	/* Auxiliary variables. */


	double ddum;          /* Double dummy */
	MKL_INT idum;         /* Integer dummy. */


	maxfct = 1;           /* Maximum number of numerical factorizations. */
	mnum = 1;         /* Which factorization to use. */
	//	msglvl = 1;           /* Print statistical information in file */
	error = 0;            /* Initialize error flag */
	/* -------------------------------------------------------------------- */
	/* .. Initialize the internal solver memory pointer. This is only */
	/* necessary for the FIRST call of the PARDISO solver. */
	/* -------------------------------------------------------------------- */
	for (int i = 0; i < 64; i++)
	{
		pt = 0;
	}

	

	sparse_checker_error_values check_err_val;
	sparse_struct pt1;


	sparse_matrix_checker_init(&pt1);
	pt1.n = n;
	pt1.csr_ia = ia;
	pt1.csr_ja = ja;
	pt1.indexing = MKL_ZERO_BASED;
	pt1.matrix_structure = MKL_GENERAL_STRUCTURE;
	pt1.print_style = MKL_C_STYLE;
	pt1.message_level = MKL_PRINT;

	check_err_val = sparse_matrix_checker(&pt1);

	printf("Matrix check details: (%d, %d, %d)\n", pt1.check_result[0], pt1.check_result[1], pt1.check_result[2]);

	if (check_err_val == MKL_SPARSE_CHECKER_NONTRIANGULAR) {
		printf("Matrix check result: MKL_SPARSE_CHECKER_NONTRIANGULAR\n");
		error = 0;
	}
	else {
		if (check_err_val == MKL_SPARSE_CHECKER_SUCCESS) { printf("Matrix check result: MKL_SPARSE_CHECKER_SUCCESS\n"); }
		if (check_err_val == MKL_SPARSE_CHECKER_NON_MONOTONIC) { printf("Matrix check result: MKL_SPARSE_CHECKER_NON_MONOTONIC\n"); }
		if (check_err_val == MKL_SPARSE_CHECKER_OUT_OF_RANGE) { printf("Matrix check result: MKL_SPARSE_CHECKER_OUT_OF_RANGE\n"); }
		if (check_err_val == MKL_SPARSE_CHECKER_NONORDERED) { printf("Matrix check result: MKL_SPARSE_CHECKER_NONORDERED\n"); }
		error = 1;
	}

	
	phase = 11;
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, 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]); fflush(0);
	/* -------------------------------------------------------------------- */
	/* .. Numerical factorization. */
	/* -------------------------------------------------------------------- */

	phase = 22;
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0)
	{
		system("pause");
		printf("\nERROR during numerical factorization: %d", error);
		exit(2);
	}
	
	double *df = NULL;
	double *da = NULL;
	MKL_INT er_d = 0;

	df = (double*)mkl_malloc(n * sizeof(double), 64);

	da = (double*)mkl_malloc(n * sizeof(double), 64);



	pardiso_getdiag(pt, df, da, &mnum, &er_d);
	printf("\nFirst\n");
	for (int i = 0; i < n; i++)
	{
		printf("d_fact[%d]= %lf d_a=%lf \n", i, df, da);
	}

	phase = 33;
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, X, &error);
	if (error != 0)
	{
			printf("\nERROR during solution: %d", error);
			exit(3);
	}
	phase = -1;           /* Release internal memory. */
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&n, &ddum, ia, ja, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);

After pardiso_getdiag i see this: d_fact[41]= 4.269444 d_a=4.625000 d_fact[42]= 4.269444 d_a=4.625000 d_fact[43]= 37.925772 d_a=-0.000000 d_fact[44]= 37.924292 d_a=-0.000000 d_fact[45]= 8.519367 d_a=-0.000000 d_fact[46]= 4.189655 d_a=-0.000000 d_fact[47]= 2.036909 d_a=-0.000000 d_fact[48]= 0.412003 d_a=-0.000000 d_fact[49]= 0.402489 d_a=-0.000000

I also attached my b vector. Maybe it would help.

[mecej4]

Since you posted only a part of the code, it is difficult to comment. Similarly, when you simply say "i don't get result that i expect" we are not told what you expect.

When you report problems with code, it is usually necessary for you to give details of MKL and compiler version, OS, CPU model, etc.

[Oleg_S_]

I'm running visual studio on win 7. I have MKL 2017 Update 2. I have intel compiler 2017.1.143. My cpu model is i7-3615QM

[Oleg_S_]

Can someone tell why pardiso_getdiag thinks that i have zeros on a diagonal?

[mecej4]

I had not used pardiso_getdiag() before you brought it to my attention. The MKL document page contains the following contradictory statements:

da: Array with a dimension of n. Contains diagonal elements of the initial matrix

and

Elements of da correspond to diagonal elements of matrix L computed during phase 22

It cannot be both (initial matrix A and L factor). So, MKL group, which is it? If the latter, how is da different from df ?

Oleg: The solution of your equation set given by the dense Lapack solver DGESV agrees with the solution given by MKL-Pardiso, as well as Pardiso-5 from www.pardiso-project.org . If you maintain that the results are incorrect you should  provide supporting arguments. Until I learn more about pardiso_getdiag() I would not not depend entirely on it to decide whether a solution is correct or not.