Community
cancel
Showing results for 
Search instead for 
Did you mean: 
47 Views

Pardiso in Fortran MKL giving incorrect answers

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

0 Kudos
4 Replies
mecej4
Black Belt
47 Views

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

47 Views

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?

47 Views

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
Black Belt
47 Views

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.