Inaccurate results while using PARDISO

Bounrajbanditt_K_ · ‎02-04-2018

Hello,

I am trying to utilize the PARDISO sparse matrix solver to test out a solution with my sparse matrix. The goal of this experiment was to see if I can replicate the same solution that I had while using MATLAB's backslash operator. The results of PARDISO and the solution using the backslash operator were no where near the same.

In order to do this, I exported the sparse matrix arrays from MATLAB into binary files in COO format. I re-imported the sparse matrix arrays into my C program in which I then converted it into CSR format for use with PARDISO using mkl_?csrcoo. I then set up PARDISO to solve it but ended up with inaccurate results. Attached is my code, the matrix bin files, a makefile (just to see what my linking and compiling parameters were) and the original MATLAB sparse matrix in .mat format

My sparse matrix is set up as a nonsymmetric real matrix. I wonder if it's possible to achieve matching results if I set up the right iparm parameters or if the problem lies else where in my program. Either way, I wish to gain insight into my mistake.

The matrix is of size 122860x122860.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>

#include "mkl.h"
#include "mkl_pardiso.h"
#include "mkl_types.h"
#include <time.h>

int main()
{
#define N 122860
#define NNZ_G 491063
#define INFO 0
#define P 64
	//***************************************************************************************************
	//* Import sparse G matrix data in COO sparse matrix format
	//***************************************************************************************************

	FILE *ptr;
	char pathG1[] = "./Gmat122860x122860/Gvalues122860x122860.bin";
	char pathG2[] = "./Gmat122860x122860/Grows122860x122860.bin";
	char pathG3[] = "./Gmat122860x122860/Gcols122860x122860.bin";

	double *gval = (double*) calloc(NNZ_G, sizeof(double));
	MKL_INT *grow = (MKL_INT*) calloc(NNZ_G, sizeof(MKL_INT));
	MKL_INT *gcol = (MKL_INT*) calloc(NNZ_G, sizeof(MKL_INT));

	double *gcsr = (double*) calloc(NNZ_G, sizeof(double));
	MKL_INT *jg = (MKL_INT*) calloc(NNZ_G, sizeof(MKL_INT));
	MKL_INT *ig = (MKL_INT*) calloc(N + 1, sizeof(MKL_INT));

	ptr = fopen(pathG1, "rb");
	fread(gval, sizeof(double), NNZ_G, ptr);
	fclose(ptr);

	ptr = fopen(pathG2, "rb");
	fread(grow, sizeof(MKL_INT), NNZ_G, ptr);
	fclose(ptr);

	ptr = fopen(pathG3, "rb");
	fread(gcol, sizeof(MKL_INT), NNZ_G, ptr);
	fclose(ptr);
	
	//***************************************************************************************************
	//* Convert sparse matrix format of G from COO to CSR - base 1 indexing
	//***************************************************************************************************

	MKL_INT job[6];
	MKL_INT n = N, nnz_g = NNZ_G;

	job[0] = 2;
	job[1] = 1;
	job[2] = 1;
	job[3] = 0;
	job[4] = 0;
	job[5] = 0;

	MKL_INT info = INFO;
	mkl_dcsrcoo(job, &n, gcsr, jg, ig, &nnz_g, gval, grow, gcol, &info);

	//***************************************************************************************************
	// Set up for PARDISO and loop
	//***************************************************************************************************

	void    *pt;
	MKL_INT iparm
;
	MKL_INT i, maxfct, mnum, mtype, idum;
	MKL_INT msglvl, error, phase;
	MKL_INT nrhs = 1;

	mtype = 11;
	maxfct = 1;        //  Maximum number of numerical factorizations. 
	mnum = 1;        //  Which factorization to use. 
	msglvl = 0;        // Print statistical information in file 
	error = 0;        //  Initialize error flag 
	
	pardisoinit(pt, &mtype, iparm);

	iparm[34] = 0; // One based indexing

	char transa = 'n';
	
	double *b0 = (double*) calloc(N, sizeof(double));
	b0[N-1] = -0.8; 
	 
	double *x0 = (double*) calloc(N, sizeof(double));  
	
	phase = 13; 
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &n, gcsr, ig, jg, &idum, &nrhs, iparm, &msglvl, b0, x0, &error);
       
	return 0; 
}

mecej4 · ‎02-05-2018

I suggest that you start with a much smaller matrix and work the bugs out.

There are many pieces that could have gone wrong: (i) exporting the data from Matlab to a file; (ii) importing the data into a C program that calls Pardiso, (iii) calling Pardiso with the correct arguments, and (iv) interpreting the results.

Questions to consider: what is the condition number of the matrix? How accurate do you expect the solution to be? Which results from Matlab did you use in deciding that the Pardiso results did not agree?

Bounrajbanditt_K_ · ‎02-05-2018

Thank you for your reply. I have tested my program on a smaller matrix that is in the MKL 2017 C developer's reference pdf under the sparse matrix format appendix (they all seem to use the same sparse matrix). I exported and reimported the matrix into my program to test with PARDISO and I also tested them with Matlab. This is the matlab script I used. I found out that the results were matching between the c program and the matlab script, both giving the solution:

x = [-0.099290780141844;-0.039716312056738;-0.019858156028369;-0.051063829787234;0.096453900709220]

Script:

aval = [1, -1, -3, -2, 5, 4, 6, 4, -4, 2, 7, 8, -5];

arow = [1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5];

acol = [1, 2, 3, 1, 2, 3, 4, 5, 1, 3, 4, 2, 5];

A = sparse(arow, acol, aval, 5, 5);

b0 = [0;0;0;0;-0.8];

x = A\b0;

I reckon I can knock out (i) and (ii). I also found out using Matlab's condest() function that my original matrix (called 'G') has a condition number of 4.8796e+07 meaning it is largely ill-conditioned. I imagine it might be in setting up PARDISO. In my original program b was an all zeros vector with the last value set to -0.8. I read the iparm array in the manual as well as the individual settings for each index and its settings yet I cannot understand how I can set it so that my matrix may be properly preconditioned. I kept experimenting and fiddling around but still got wild results. I wish to know what I can do to set PARDISO to have matching results to what Matlab has when doing G\b. Here are the first 10 solution values of both Matlab and PARDISO for my original program:

MATLAB results :

2.65767860454835e-16

6.93755753591888e-16

6.69682277981033e-16

-7.80650219198111e-16

4.26942724439305e-16

-1.94635032810143e-15

2.42380892307184e-16

-1.75193726393755e-16

1.23125245192626e-15

1.49153731417372e-15

PARDISO:

vec[1] = 6.828214E+00

vec[2] = 1.418508E+02

vec[3] = 1.958590E+02

vec[4] = 7.672455E+02

vec[5] = -1.905387E+03

vec[6] = 4.334156E+00

vec[7] = 7.674021E-01

vec[8] = 2.891470E-01

vec[9] = -7.786293E-01

vec[10] = -1.586327E+00

mecej4 · ‎02-06-2018

Further investigation showed that your matrix is somewhat peculiar. As you noted, it is not well-conditioned, and you may have only 7-8 correct digits in the solution. Furthermore, of the 122,860 elements of the solution x, only 12,324 differ from 0 by more than 10^-10. Furthermore, all of these non-zero solution elements equal the only nonzero value in b, namely, 0.8! Is this what you expect?

I think that Pardiso ought to do better than it does on this problem, but I do not know if that is an issue that is worth pursuing. You may try more steps of iterative refinement to see if that cures the problem.

As you may know, Matlab uses SuperLU for A\b operations on sparse matrices, and you may use with your code the freely downloadable SuperLU solver from http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ .