Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
This community is designed for sharing of public information. Please do not share Intel or third-party confidential information here.
6589 Discussions

## Pardiso in Fortran MKL giving incorrect answers

Beginner
187 Views

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

4 Replies
Black Belt
187 Views

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

Beginner
187 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?

Beginner
187 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

Black Belt
187 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.