topic The 3 X 3 matrix that you in Intel® oneAPI Math Kernel Library

Pardiso in Fortran MKL giving incorrect answers

matthew_bernardnrc_g — Tue, 12 Jul 2016 20:55:12 GMT

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

The 3 X 3 matrix that you

mecej4 — Tue, 12 Jul 2016 22:12:27 GMT

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

Geez, that was an

matthew_bernardnrc_g — Wed, 13 Jul 2016 12:48:42 GMT

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?

I wanted to report that I got

matthew_bernardnrc_g — Wed, 13 Jul 2016 13:58:00 GMT

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/)

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

There is no separate F95

mecej4 — Wed, 13 Jul 2016 14:17:05 GMT

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.