topic Issue with Intel OneAPI MKL PARDISO in Intel® oneAPI Math Kernel Library

Issue with Intel OneAPI MKL PARDISO

Jeremie_V_ — Thu, 03 Jun 2021 17:05:25 GMT

Dear,

Using ifort 2021.2.0, MKL Pardiso provides a wrong answer when a permutation vector is provided, while it provides a correct answer when no permutation vector is provided.

The correct answer is computed in both cases with Intel ifort 17.0.0

I checked the API, but it doesn't seem that it changed since ifort 17.0.0

I wonder if the problem is in my code, or in the library.

Attached is a small program that reproduces the issue.

In advance thank you.

Jeremie

```fortran

include 'mkl_pardiso.f90'

program test
use iso_fortran_env, only: int32, real64, output_unit
use mkl_pardiso, only: MKL_PARDISO_HANDLE
implicit none
integer(int32), allocatable :: crs_i(:), crs_j(:)
integer(int32), allocatable :: perm(:)
integer(int32) :: getdim
real(real64), allocatable :: crs_a(:)
real(real64), allocatable :: rhs(:), x(:), x_true(:)

logical :: lperm = .true.

!Pardiso variables
integer(kind=int32) :: i, error
integer(kind=int32)::nrhs
integer(kind=int32)::idum(1)
real(kind=real64)::ddum(1)
integer(kind=int32)::maxfct
integer(kind=int32)::mnum
integer(kind=int32)::mtype
integer(kind=int32)::msglvl
integer(kind=int32)::phase
integer(kind=int32)::solver
integer(kind=int32),allocatable::iparm(:)
type(MKL_PARDISO_HANDLE),allocatable::pt(:)

crs_i = [1, 4, 6, 7, 8, 9]
crs_j = [1 , 2 , 3 , 2 , 3 , 3 , 4 , 5]
crs_a = [10._real64, 1._real64, 1.2_real64, 10._real64, 0.3_real64, 13.4_real64, 0._real64, 1._real64]
rhs = [1._real64, 2._real64, 3._real64, 0._real64, 5._real64]

x_true = [5.543360619502753d-002, 0.188015426594107_real64, 0.214707092879682_real64&
,0._real64, 5._real64]

allocate(x, mold = rhs)

perm = [3, 2, 1, 5, 4]

getdim = size(crs_i) - 1

call printsquare(crs_i, crs_j, crs_a)
! print*,getdim
! print*,crs_i
! print*,crs_j
! print*,crs_a
! print*,rhs

maxfct = 1
mnum = 1
mtype = -2
solver = 0
msglvl = 1
nrhs = 1
idum = 0
ddum = 0

allocate(iparm(64))
iparm = 0

allocate(pt(64))
do i=1,64
pt(i)%dummy = 0
enddo

call pardisoinit(pt,mtype,iparm)

!Ordering and factorization
phase = 12
iparm(2) = 3
iparm(19) = -1
iparm(24) = 1
iparm(27) = 1
iparm(28) = 0
write(*,'(a)')' Start ordering and factorization'

if(lperm)then
!provide a permutation vector
iparm(5) = 1
iparm(31) = 0
iparm(36) = 0
call pardiso(pt,maxfct,mnum,mtype,phase,&
getdim,crs_a,crs_i,crs_j,&
perm,nrhs,iparm,msglvl,ddum,ddum,error)
else
!don't provide a permutation vector
call pardiso(pt,maxfct,mnum,mtype,phase,&
getdim,crs_a,crs_i,crs_j,&
idum,nrhs,iparm,msglvl,ddum,ddum,error)
endif
call checkpardiso(phase,error)

iparm(27)=0 !disable Pardiso checker

!Solving
phase=33
call pardiso(pt,maxfct,mnum,mtype,phase,&
getdim,crs_a,crs_i,crs_j,&
idum,nrhs,iparm,msglvl,rhs,x,error)
call checkpardiso(phase,error)

write(*,'(a,*(f10.4,x))')'x = ',x

if(sum(abs(x - x_true))/size(x) > 1.d-6)then
write(*,'(a,*(f10.4,x))')'x_true = ',x_true
print*,'Rel abs error = ',sum(abs(x - x_true))/size(x)
error stop 'WRONG ANSWER FROM PARDISO'
else
print*,'CORRECT ANSWER'
endif

contains

subroutine checkpardiso(phase,error,un)
integer(kind=int32),intent(in)::phase,error
integer(kind=int32),intent(in),optional::un

integer(kind=int32)::unlog

unlog=output_unit
if(present(un))unlog=un

if(error.ne.0)then
write(unlog,'(2(a,i0))')' The following error for phase ',phase,' was detected: ',error
error stop
endif

end subroutine

subroutine printsquare(ia, ja, a)
integer(int32), intent(in) :: ia(:), ja(:)
real(real64), intent(in) :: a(:)

integer :: i, j
real(real64), allocatable :: dense(:,:)

allocate(dense(size(ia)-1,size(ia)-1), source = 0._real64)

do i = 1, size(ia)-1
do j = ia(i), ia(i+1) - 1
dense(i, ja(j)) = a(j)
dense(ja(j), i) = a(j)
enddo
enddo

do i = 1, size(ia)-1
write(*,'(*(f10.4,x))')dense(:,i)
enddo

end subroutine

end program test

```

Re:Issue with Intel OneAPI MKL PARDISO

ShivaniK_Intel — Fri, 04 Jun 2021 11:05:54 GMT

Hi,

Thanks for reaching out to us.

We are able to reproduce the issue at our end. We are working on it and will get back to you soon.

Thanks & Regards

Shivani

Re:Issue with Intel OneAPI MKL PARDISO

Gennady_F_Intel — Mon, 07 Jun 2021 02:11:53 GMT

Jeremie,

is that win or lin os?

Re:Issue with Intel OneAPI MKL PARDISO

Gennady_F_Intel — Mon, 07 Jun 2021 03:11:31 GMT

Jeremy,

I checked versions 2018,2019, 2020 and the latest 2021 - everywhere the problem is reproduce. The latest version which works correct is 2017u5. We will investigate the cause of the issue and keep this thread inform with the status.

thanks for the report.

Re: Re:Issue with Intel OneAPI MKL PARDISO

Jeremie_V_ — Mon, 07 Jun 2021 08:07:39 GMT

Thank you for the update.

For ifort 2021: Fedora 31 and 33

For ifort 2017: Red Hat

Yours sincerely,

Jeremie

Re: Re:Issue with Intel OneAPI MKL PARDISO

Jeremie_V_ — Mon, 05 Jul 2021 10:18:29 GMT

Hello,

Thank you. Meanwhile, are there any ways to avoid this issue? Maybe with a different combination of options?

Jeremie

Re: Re:Issue with Intel OneAPI MKL PARDISO

Kirill_V_Intel — Mon, 05 Jul 2021 19:48:54 GMT

Hi Jeremie!

There are two temporary workarounds before we release the new version with the fix:

1) set iparm(24)=0 (classic factorization) or 10 (improved two-level factorization).

2) if your system is not small, use more than one thread. if your system is small, try to use more than one thread and set MKL_DYNAMIC to false.

I've quickly checked that both ideas worked.

Best,
Kirill

Re: Re:Issue with Intel OneAPI MKL PARDISO

Jeremie_V_ — Wed, 07 Jul 2021 14:00:44 GMT

Dear,

Thank you for your answer.

My usual systems are not small in terms of equations, but highly sparse (that is a few million equations).

I will test both options, and come back if needed!

Thanks.

Yours sincerely,

Jeremie

Re:Issue with Intel OneAPI MKL PARDISO

Gennady_F_Intel — Sat, 02 Oct 2021 04:53:35 GMT

Jeremie,

the issue is fixed into the latest MKL 2021.4 which is available for download.

here is the result I see when linked against this update:

$ ./a.out

10.0000 1.0000 1.2000 0.0000 0.0000

1.0000 10.0000 0.3000 0.0000 0.0000

1.2000 0.3000 13.4000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 1.0000

Start ordering and factorization

Memory allocated on phase 11 -0.0005 MB

Number of non-zeros in L 11

Number of non-zeros in U 1

Memory allocated on phase 22 -0.0001 MB

Percentage of computed non-zeros for LL^T factorization

81 % 90 % 100 %

x = 0.0554 0.1880 0.2147 0.0000 5.0000

CORRECT ANSWER

$ echo echo $MKLROOT

echo /opt/intel/oneapi/mkl/2021.4.0

Re: Issue with Intel OneAPI MKL PARDISO

Jeremie_V_ — Mon, 04 Oct 2021 20:10:14 GMT

Great! Thank you!

Re:Issue with Intel OneAPI MKL PARDISO

Gennady_F_Intel — Tue, 05 Oct 2021 07:34:41 GMT

This issue has been resolved and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.