Pardiso test example fail

Jared_W_ — Fri, 21 Feb 2014 15:27:02 GMT

I am trying to debug my implementation of PARDISO, but unfortunately I cannot get even a simple example to work.

I created a 6x6 matrix with 16 non-zero elements and used the mkl ddnscsr routine to obtain my acsr, ia, and ja vectors to feed into the PARDISO routine. Whenever i plug things into the function call nothing happens; my code just stops running after the call to PARDISO.

Here is my code:

double *tempB;

MKL_INT lda, info, tempM, tempN;
lda = 6;

tempM = 6;

tempN = 6;

       tempB = (double*)mkl_malloc(tempM*sizeof(double),16);
       if (tempB == NULL)
       {
           cout << ">>> error allocating tempB" << endl;
           return (0);
       }

       tempB[0] = 2.0;
       tempB[1] = -6.0;
       tempB[2] = 3.0;
       tempB[3] = 1.0;
       tempB[4] = 4.0;
       tempB[5] = -1.0;

       double *tempSolution;
       tempSolution = (double*) malloc (sizeof (double)*(tempM));

       /*   MKL Sparse formatting   */
       double *adns, *acsr;
       adns = (double*)mkl_malloc(tempM*tempN*sizeof(double),16);
       if (adns == NULL)
       {
           cout << ">>> error allocating adns" << endl;
           return (0);
       }

       // Count non-zero elements
       nzElements = 0;
       for(rows = 0; rows < tempM; rows++) {
           for(cols = 0; cols < tempN; cols++) {
               if(tempA[rows][cols] != 0) {
                   adns[cols + tempN*rows] = tempA[rows][cols];
                   nzElements++;
               }
               else {
                   adns[cols + tempN*rows] = 0.0;
               }
           }
       }

       acsr = (double*)mkl_malloc(nzElements*sizeof(double),16);
       if (acsr == NULL)
       {
           cout << ">>> error allocating acsr" << endl;
           return (0);
       }

       MKL_INT *ia, *ja, *job;
       ia = (MKL_INT *)mkl_malloc((tempM + 1)*sizeof(MKL_INT),16);
       if (ia == NULL)
       {
           cout << ">>> error allocating ia" << endl;
           return (0);
       }
       ja = (MKL_INT*)mkl_malloc(nzElements*sizeof(MKL_INT),16);
       if (ja == NULL)
       {
           cout << ">>> error allocating ja" << endl;
           return (0);
       }
       job = (MKL_INT*)mkl_malloc(6*sizeof(MKL_INT),16);
       if (job == NULL)
       {
           cout << ">>> error allocating ja" << endl;
           return (0);
       }

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

       // Converts a sparse matrix into a csr format
       mkl_ddnscsr (job, &tempM, &tempN, adns, &lda, acsr, ja, ia, &info);
       if(info != 0) {
           cout << "error with mkl_ddnscsr" << endl;
           cout << "info # " << info << endl;
       }

       //   Call pardiso solver
       MKL_INT pt[64], iparm[64];
       for(i = 0; i < 64; i++) {
           pt = 0;
           iparm = 0;
       }
       MKL_INT *perm;
       perm = (MKL_INT*)mkl_malloc(tempM*sizeof(MKL_INT),16);
       if (perm == NULL)
       {
           cout << ">>> error allocating perm" << endl;
           return (0);
       }
       MKL_INT maxfct, mnum, mtype, phase, nrhs, msglvl, error;
       maxfct = 1;
       mnum = 1;
       mtype = 11;
       nrhs = 1;
       msglvl = 1;
       iparm[26] = 1;
       iparm[34] = 1;

       //   Pardiso Direct Solver
       phase = 11;
       pardiso (pt, &maxfct, &mnum, &mtype, &phase, &tempN, acsr, ia, ja, perm, &nrhs, iparm, &msglvl, tempB, tempSolution, &error);
       if (error != 0) {
           cout << "ERROR during numerical factorization: " << error << endl;
       }

Everything is fine till my code calls the PARDISO function, and then I get nothing...

Any insight as to what might be happening would be great! Sorry if I have forgotten to mention anything.

Jared

What value of error has been

Gennady_F_Intel — Sat, 22 Feb 2014 03:25:09 GMT

What value of error has been returned by Pardiso?

I see garbage into acsr

Gennady_F_Intel — Sat, 22 Feb 2014 04:57:37 GMT

I see garbage into acsr array -- pls check conversion from dense to scr

vector > tempA

Jared_W_ — Sat, 22 Feb 2014 17:33:08 GMT

vector<vector<double> > tempA(6, vector<double>(6, 0));
   for(i = 0; i < 6; i++) {
       for(j = 0; j < 6; j++) {
           tempA = 0.0;
       }
   }
   tempA[0][0] = 1.0;
   tempA[0][2] = 5.0;
   tempA[0][3] = 7.0;
   tempA[1][1] = 12.0;
   tempA[1][4] = 3.0;
   tempA[2][2] = 4.0;
   tempA[2][3] = 2.0;
   tempA[2][5] = 1.0;
   tempA[3][0] = 7.0;
   tempA[3][3] = 1.0;
   tempA[4][1] = 4.0;
   tempA[4][4] = 8.0;
   tempA[4][5] = 5.0;
   tempA[5][0] = 9.0;
   tempA[5][4] = 4.0;
   tempA[5][5] = 3.0;

Woops, I forgot to add the main matrix tempA at the top. With this added everything turns out right for the formatting routine, but the problem is the PARDISO function. I do not get any error code that is why I am not sure how to proceed.

try to set iparm[0] = 1;

Gennady_F_Intel — Sun, 23 Feb 2014 06:53:22 GMT

try to set iparm[0] = 1;

This was perfect! Thank you

Jared_W_ — Thu, 27 Feb 2014 14:29:47 GMT

This was perfect! Thank you Gennady. Thanks to you I got the small example working which in turn made my full program work and now my finite element fluid flow program solves 17K dof in 17 s instead of 6min.

You should know you are incredible and I have been stuck on this for a really long time. I appreciate you taking the time immensely!

Have a good day.

Jared

topic try to set iparm[0] = 1; in Intel® oneAPI Math Kernel Library

Pardiso test example fail

What value of error has been

I see garbage into acsr

vector > tempA

try to set iparm[0] = 1;

This was perfect! Thank you