Pardiso Returning Wrong Results

[wagner_b_1]

Hello dear friends,

I am trying to solve a problem using the finite volume method and I am having some troubles configurating  the Pardiso to solve my linear system.

The system that I am trying to solve has 12 elements, and it is given by:

 ia = (/  1, 10, 19, 28, 40, 52, 64, 76, 88, 100, 109, 118, 127 /)

 jac = (/  1,   2,   3,   1,   2,   3,   1,   2,   3,   4,   5,   6,   4,   5,   6, &
           4,   5,   6,   7,   8,   9,   7,   8,   9,   7,   8,   9,   1,   2,   3, &
           1,   2,   3,   1,   2,   3,   4,   5,   6,   4,   5,   6,   4,   5,   6, &
           7,   8,   9,   7,   8,   9,   7,   8,   9,   10,  11,  12,   10,  11,  12, &
           10,  11,  12,   1,   2,   3,   1,   2,   3,   1,   2,   3,   4,   5,   6, &
           4,   5,   6,   4,   5,   6,   7,   8,   9,   7,   8,   9,   7,   8,   9, &
           10,  11,  12,   10,  11,  12,   10,  11,  12,   4,   5,   6,   4,   5,   6, &
           4,   5,   6,   7,   8,   9,   7,   8,   9,   7,   8,   9,   10,  11,  12, &
           10,  11,  12,   10,  11,  12 /)

 a = (/  4.19829745d-09,  -4.66477495d-06,  -1.04957436d-07,  -7.85398163d-10, &
         0.00000000d+00,   3.36766430d-05,   3.92499666d-05,   7.14315468d-18, &
        -1.37394825d-03,   0.00000000d+00,   0.00000000d+00,   1.04957436d-07, &
         7.85398163d-10,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   1.37394825d-03,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,  -1.04957436d-07,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,  -1.37394825d-03, &
         4.19829745d-09,  -4.66477495d-06,   6.29744618d-08,  -7.85398163d-10, &
         0.00000000d+00,   2.02059858d-05,   3.92499666d-05,   7.14315468d-18, &
         8.24368947d-04,   0.00000000d+00,   0.00000000d+00,  -0.00000000d+00, &
         7.85398163d-10,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,  -6.29744618d-08,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
        -8.24368947d-04,   4.19829745d-09,  -9.32954989d-07,   2.09914873d-08, &
        -7.85398163d-10,   0.00000000d+00,   6.73532860d-06,   3.92499666d-05, &
         1.78578867d-18,   2.74789649d-04,   0.00000000d+00,   0.00000000d+00, &
        -0.00000000d+00,   7.85398163d-10,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
         0.00000000d+00,   0.00000000d+00,  -2.09914873d-08,  -7.85398163d-10, &
         0.00000000d+00,   0.00000000d+00,   0.00000000d+00,   0.00000000d+00, &
        -2.74789649d-04,   4.19829745d-09,  -9.32954989d-07,   2.09914873d-08, &
         7.85398163d-10,   0.00000000d+00,   6.73532860d-06,   3.92499666d-05, &
         1.78578867d-18,   2.74789649d-04 /)

 b  = (/   0.00000000d+00,  0.00000000d+00,  0.00000000d+00,  0.00000000d+00, &
          -0.00031416d+00,  0.00000000d+00,  0.00000000d+00,  0.00000000d+00, &
           0.00000000d+00,  0.00000000d+00,  0.00000000d+00,  0.00000000d+00  /)

I have analyzed and solved this system with Matlab and the conditioning number is about 1E3, in others words, it is not too high that cannot be solvedd using direct solvers.

I am setting the Pardiso with this parameters:

maxfct=1
        mnum=1
        mtype=11 ! real and nonsymmetric
        msglvl=1 ! NOT prints statistical information to the screen.
        perm = 1

        iparm=0
        iparm(1) = 1 ! no solver default
        iparm(2) = 2 ! fill-in reordering from METIS
        iparm(3) = 1 ! numbers of processors
        iparm(4) = 0 ! no iterative-direct algorithm
        iparm(5) = 0 ! no user fill-in reducing permutation
        iparm(6) = 0 ! =0 solution on the first n compoments of x
        iparm(7) = 0 ! not in use
        iparm(8) = 9 ! numbers of iterative refinement steps
        iparm(9) = 0 ! not in use
        iparm(10) = 13 ! perturbe the pivot elements with 1E-13
        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
        iparm(12) = 0 ! not in use
        iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate ccuracy
        iparm(14) = 0 ! Output: number of perturbed pivots
        iparm(15) = 0 ! not in use
        iparm(16) = 0 ! not in use
        iparm(17) = 0 ! not in use
        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
        iparm(19) = -1 ! Output: Mflops for LU factorization
        iparm(20) = 0 ! Output: Numbers of CG Iterations
        iparm(27) = 1

        phase=13 ! Analysis, numerical factorization, solve, iterative refinement
        call pardiso(pt, maxfct, mnum, mtype, phase, m, a, ia, jac, perm, 1, iparm, msglvl, b, x, error)

        phase=-1
        call pardiso(pt, maxfct, mnum, mtype, phase, m, A, rown, col_n, perm, 1, iparm, msglvl, b, x, error)
        call mkl_free_buffers

The thing is that when I compare the Pardiso result with the Matlab one, the results does not match because Pardiso is returning strange values.

Please Help me, I cannot continue my reserach without solving this problem.

Thanks,

Wagner Barros

[mecej4 · Accepted solution]

I think that your matrix data are wrong. For instance, in row-1 you should have nine entries. The column indices of these entries should be distinct and in ascending order. What you have, however, are the values [1,   2,   3,   1,   2,   3,   1,   2,   3]. Similar errors exist in the data for subsequent rows.

You may find it helpful to set iparm(27) = 1 so that Pardiso can run some checks on the data.

Your matrix is not sparse (126 out of 144 entries are not zero). Have you considered using a dense matrix routine such as GESV?

[wagner_b_1]

Thanks a lot Mecej4,

Our matrix entries are really wrong, it is because we are using a routine to change the triplets entries to compressed row format, and this routine was not working very well. 

 About these matrix, it is a model created with only four discrete cells, used to configurate the solver. The real model will have more than 1.000 cells generating a real sparse matrix with more than 10.000 rows.

I would like to thank you for your patience, and inform that I finally made the Pardiso work.

Thanks, Wagner Barros.

[mecej4]

MKL contains a routine for converting between the COO and CSR representations: https://software.intel.com/en-us/node/468628 .