Pardiso test example fail

Pardiso test example fail

Ritratto di 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[i] = 0;
            iparm[i] = 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

6 post / 0 new
Ultimo contenuto
Per informazioni complete sulle ottimizzazioni del compilatore, consultare l'Avviso sull'ottimizzazione
Ritratto di Gennady Fedorov (Intel)

What  value of error has been returned by Pardiso? 

Ritratto di Gennady Fedorov (Intel)

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

Ritratto di Jared W.

vector<vector<double> > tempA(6, vector<double>(6, 0));
    for(i = 0; i < 6; i++) {
        for(j = 0; j < 6; j++) {
            tempA[i][j] = 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.

 

Ritratto di Gennady Fedorov (Intel)

try to set iparm[0] = 1; 

Ritratto di 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

Accedere per lasciare un commento.