Pardiso test example fail

[Jared_W_]

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

[Gennady_F_Intel]

What  value of error has been returned by Pardiso? 

[Gennady_F_Intel]

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

[Jared_W_]

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.

[Gennady_F_Intel]

try to set iparm[0] = 1; 

[Jared_W_]

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