Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
6436 Discussions

Issue with Intel OneAPI MKL PARDISO

Jeremie_V_
Beginner
727 Views

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

```

Labels (1)
0 Kudos
10 Replies
ShivaniK_Intel
Moderator
704 Views

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


Gennady_F_Intel
Moderator
674 Views

Jeremie,

is that win or lin os?


Gennady_F_Intel
Moderator
670 Views

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.



Jeremie_V_
Beginner
646 Views

Thank you for the update.

 

For ifort 2021: Fedora 31 and 33

For ifort 2017: Red Hat

 

Yours sincerely,

 

Jeremie

Jeremie_V_
Beginner
483 Views

Hello,

 

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

 

Jeremie

Kirill_V_Intel
Employee
468 Views

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

or

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

Jeremie_V_
Beginner
444 Views

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

Gennady_F_Intel
Moderator
306 Views

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



Jeremie_V_
Beginner
274 Views

Great! Thank you!

Gennady_F_Intel
Moderator
257 Views

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. 



Reply