Pardiso in Fortran MKL giving incorrect answers

[matthew_bernardnrc_g]

Hello,

I am trying to understand how the Intel MKL works and so i created a small test code that I have cobbled together from the internet to teach myself how the libraries behave so that I can implement the capability in other projects. My code compiles, runs and reports to zero errors but returns the wrong answer for my system. I have to be doing something wrong but I have no idea what that thing is. My test code is below. 

PROGRAM PardisoTest

use MKL_PARDISO

IMPLICIT NONE

   TYPE(MKL_PARDISO_HANDLE), ALLOCATABLE, DIMENSION(:) :: PT
   INTEGER, ALLOCATABLE, DIMENSION(:) :: idum
   INTEGER, ALLOCATABLE, DIMENSION(:) :: ivect
   INTEGER, ALLOCATABLE, DIMENSION(:) :: jvect
   REAL, ALLOCATABLE, DIMENSION(:) :: b
   REAL, ALLOCATABLE, DIMENSION(:) :: a
   REAL, ALLOCATABLE, DIMENSION(:) :: x
   REAL, ALLOCATABLE, DIMENSION(:) :: solution
   REAL, ALLOCATABLE, DIMENSION(:) :: error
   REAL, ALLOCATABLE, DIMENSION(:) :: ddum
   INTEGER :: maxfct, mnum, mtype, phase, nrhs, ierr, msglvl, i, nnz, nequations
   INTEGER, DIMENSION(64) :: iparm

   nrhs =1
   maxfct = 1
   mnum = 1

   nnz = 9
   nequations = 3

   do i = 1, 64
      iparm(i) = 0
   end do
!
   iparm(1) = 0 ! solver default
   error = 0 ! initialize error flag
   msglvl = 1 ! print statistical information
   mtype = 11 ! real unsymmetric --> as an example

   ALLOCATE(PT(64))
   ALLOCATE(idum(0))
   ALLOCATE(ddum(0))
   ALLOCATE(b(nEquations))
   ALLOCATE(x(nEquations))
   ALLOCATE(a(nnz))
   ALLOCATE(ivect(nequations+1))
   ALLOCATE(jvect(nnz))
   ALLOCATE(error(nEquations))

   do i = 1, 64
      PT(i)%DUMMY = 0
   end do

   a=(/1,2,3,4,5,6,7,8,9/)
   b=(/14,32,50/)
   solution = (/1.0, 2.0, 3.0/)
   jvect = (/1,2,3,1,2,3,1,2,3/)
   ivect = (/1,4,7,10/)

   phase = 11
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, ddum, ddum, ierr)
!
   phase = 22
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, ddum, ddum, ierr)
!
   phase = 33
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, b, x, ierr)

   error = solution - x

END PROGRAM PardisoTest

Any help that someone could provide would be greatly appreciated.

Matt

[mecej4]

The 3 X 3 matrix that you chose is singular -- the arithmetic mean of rows 1 and 3 equals row 2.

[matthew_bernardnrc_g]

Geez, that was an embarrassing mistake. I made a change to my system so that a(3,3,) = 10 instead of 9 and b(3) = 53. instead of 50. This fixes the linear dependence of my system and maintains the solution vector for this trivial system. Unfortunately, now my solution vector is just b = (/0.0, Nan, 0.0/). Is there anything else I need to change?

[matthew_bernardnrc_g]

I wanted to report that I got my program to work. It appears that the pardiso subroutine requires double precision vectors for the a, b, and x variables. I tested this by declaring an arbitrary type parameter and specifying different precisions for the variables. While it does surprise me that this is a requirement of the capability, it is not something that will affect its use in my other projects. I was not able to find any simple examples using the pardiso interface in fortran 90 so I am pasting my working code below. Thank you again for your help.

Matt

PROGRAM PardisoTest

use MKL_PARDISO

IMPLICIT NONE

   INTEGER, PARAMETER :: sdk = SELECTED_REAL_KIND(13, 307)
   INTEGER, PARAMETER :: sik = KIND(10000000)
!
   TYPE(MKL_PARDISO_HANDLE), ALLOCATABLE, DIMENSION(:) :: PT
   INTEGER(sik), ALLOCATABLE, DIMENSION(:) :: idum
   INTEGER(sik), ALLOCATABLE, DIMENSION(:) :: ivect
   INTEGER(sik), ALLOCATABLE, DIMENSION(:) :: jvect
   REAL(sdk), ALLOCATABLE, DIMENSION(:) :: b
   REAL(sdk), ALLOCATABLE, DIMENSION(:) :: a
   REAL(sdk), ALLOCATABLE, DIMENSION(:) :: x
   REAL(sdk), ALLOCATABLE, DIMENSION(:) :: solution
   REAL(sdk), ALLOCATABLE, DIMENSION(:) :: error
   REAL(sdk), ALLOCATABLE, DIMENSION(:) :: ddum
   INTEGER(sik) :: maxfct, mnum, mtype, phase, nrhs, ierr, msglvl, i, nnz, nequations
   INTEGER(sik), DIMENSION(64) :: iparm

   nrhs =1
   maxfct = 1
   mnum = 1

   nnz = 9
   nequations = 3

   do i = 1, 64
      iparm(i) = 0
   end do
!
   iparm(1) = 0 ! solver default
   error = 0 ! initialize error flag
   msglvl = 1 ! print statistical information
   mtype = 11 ! real unsymmetric --> as an example

   ALLOCATE(PT(64))
   ALLOCATE(idum(0))
   ALLOCATE(ddum(0))
   ALLOCATE(b(nEquations))
   ALLOCATE(x(nEquations))
   ALLOCATE(a(nnz))
   ALLOCATE(ivect(nequations+1))
   ALLOCATE(jvect(nnz))
   ALLOCATE(error(nEquations))
   ALLOCATE(solution(nEquations))

   do i = 1, 64
      PT(i)%DUMMY = 0
   end do

   a=(/1.0_sdk,2.0_sdk,3.0_sdk,4.0_sdk,5.0_sdk,6.0_sdk,7.0_sdk,8.0_sdk,10.0_sdk/)
   b=(/14.0_sdk,32.0_sdk,53.0_sdk/)
   solution = (/1.0, 2.0, 3.0/)
   jvect = (/1,2,3,1,2,3,1,2,3/)
   ivect = (/1,4,7,10/)

   phase = 11
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, ddum, ddum, ierr)
!
   phase = 22
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, ddum, ddum, ierr)
!
   phase = 33
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, b, x, ierr)

   error = solution - x

   phase = -1
   CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, ddum, ddum, ierr)

END PROGRAM PardisoTest

[mecej4]

There is no separate F95 interface to Pardiso. Secondly, to use non-default type matrices, you have to set IPARM() accordingly. For example, iparm(28) = 1 is needed if you want to use matrices whose elements are 32-bit reals. The calling sequence of Pardiso is rather complex, and there are a large number of options to understand and to pass correct values for.