topic Can someone tell why pardiso in Intel® oneAPI Math Kernel Library

pardiso_getdiag problem

Oleg_S_ — Fri, 03 Mar 2017 22:43:01 GMT

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.

Your matrix is fine, and I

mecej4 — Sat, 04 Mar 2017 02:04:46 GMT

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.

*msglvl = 0; // Не выводить

Oleg_S_ — Sat, 04 Mar 2017 08:41:57 GMT

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.

Since you posted only a part

mecej4 — Sat, 04 Mar 2017 11:55:29 GMT

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.

I'm running visual studio on

Oleg_S_ — Sat, 04 Mar 2017 16:07:09 GMT

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

Can someone tell why pardiso

Oleg_S_ — Sun, 05 Mar 2017 11:47:03 GMT

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

I had not used pardiso

mecej4 — Sun, 05 Mar 2017 13:47:00 GMT

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.