- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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; }

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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]

PARDISO:

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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/ .

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page