Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

Issue with Intel OneAPI MKL PARDISO

Jeremie_V_
Beginner
2,114 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
1 Solution
Gennady_F_Intel
Moderator
1,693 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



View solution in original post

0 Kudos
10 Replies
ShivaniK_Intel
Moderator
2,091 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
2,061 Views

Jeremie,

is that win or lin os?


0 Kudos
Gennady_F_Intel
Moderator
2,057 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.



0 Kudos
Jeremie_V_
Beginner
2,033 Views

Thank you for the update.

 

For ifort 2021: Fedora 31 and 33

For ifort 2017: Red Hat

 

Yours sincerely,

 

Jeremie

0 Kudos
Jeremie_V_
Beginner
1,870 Views

Hello,

 

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

 

Jeremie

0 Kudos
Kirill_V_Intel
Employee
1,855 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

0 Kudos
Jeremie_V_
Beginner
1,831 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

0 Kudos
Gennady_F_Intel
Moderator
1,694 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



0 Kudos
Jeremie_V_
Beginner
1,661 Views

Great! Thank you!

0 Kudos
Gennady_F_Intel
Moderator
1,644 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. 



0 Kudos
Reply