topic Civil engineer trying to program using MKL library Pardiso in Intel® oneAPI Math Kernel Library

Civil engineer trying to program using MKL library Pardiso

Brad_Doudican — Sat, 01 Oct 2011 01:27:02 GMT

I'm in the final stages of a finite element program, and need to solve the basic Ax=B matrix for x. I'm attempting to use MKL Pardiso, as the matrix is large and sparse. I've got the code operating up until the actual call to pardiso. Error I'm getting is "Unhandled exception at 0x000000014018b14b in doudican_current.exe: 0xC0000005: Access violation writing location 0x0000000050000163."

Any thoughts? See code below.

!***********************************************************************

SUBROUTINE DOUDICAN_SOLVE (GLK, BL, SOLUTION)

USE all

USE MKL_PARDISO

implicit none

REAL(KIND=DBL), INTENT(IN), DIMENSION(NEQ, NEQ) :: GLK

REAL(KIND=DBL), INTENT(IN), DIMENSION(NEQ) :: BL

REAL(KIND=DBL), INTENT(OUT), DIMENSION(NEQ) :: SOLUTION

REAL(KIND=DBL), ALLOCATABLE, DIMENSION(:) :: Astar

INTEGER, ALLOCATABLE, DIMENSION(:) :: JA

INTEGER, DIMENSION(NEQ) :: IA

INTEGER, DIMENSION(64) :: PT

INTEGER, DIMENSION(64) :: IPARM

INTEGER, DIMENSION(NEQ) :: PERM

INTEGER :: MAXFCT, MNUM, MTYPE, PHASE, NRHS, MSGLVL, ERROR, M, N, &

COUNT, COUNT1, COUNT2, COUNT2OLD, I, J

! MODIFY GLK TO BECOME Astar

COUNT=0

DO I=1,NEQ

DO J=1,NEQ

IF (GLK(I,J)/=0) THEN

COUNT=COUNT+1

END IF

END DO

ALLOCATE ( Astar(COUNT), JA(COUNT) )

COUNT1=0

COUNT2=0

COUNT2OLD=0

DO I=1,NEQ

COUNT2=COUNT2+1

DO J=1,NEQ

IF (GLK(I,J)/=0) THEN

COUNT1=COUNT1+1

Astar(COUNT1)=GLK(I,J)

JA(COUNT1)=J

IF (COUNT2OLD

IA(COUNT2)=COUNT1

COUNT2OLD=COUNT2

END IF

END DO

! DEFINE PARAMETERS

DO I=1,64

PT(I)=0

IPARM(I)=0

END DO

MAXFCT =1

MNUM = 1

MTYPE = 2

PHASE = 11

NRHS=1

MSGLVL = 0

! CALL THE SPARSE MATRIX SOLVER PARDISO

CALL PARDISO_ (PT, MAXFCT, MNUM, MTYPE, PHASE, NEQ, Astar, IA, JA, PERM, IPARM, MSGLVL, BL, SOLUTION, ERROR)

RETURN

END SUBROUTINE DOUDICAN_SOLVE

!***********************************************************************

Civil engineer trying to program using MKL library Pardiso

mecej4 — Sat, 01 Oct 2011 01:40:53 GMT

Recompile with /traceback and run again to correlate the IP with source file name and line number. We don't even know in which routine the crash occurred from the information given.

Show the declarations of the arguments to DOUDICAN_SOLVE that were made in the routine that called it.

Civil engineer trying to program using MKL library Pardiso

Brad_Doudican — Sat, 01 Oct 2011 02:27:06 GMT

You can see how little computer science experience I have. Please tell me if this is helpful or what you were looking for.

LRC= 704

forrtl: severe (157): Program Exception - access violation

Image PC Routine Line Source

doudican_current. 000000014018B14B Unknown Unknown Unknown

doudican_current. 000000014018AF31 Unknown Unknown Unknown

doudican_current. 0000000140093B6F DOUDICAN_SOLVE 4599 doudican.f90

doudican_current. 0000000140007F95 MAIN__ 404 doudican.f90

doudican_current. 000000014051C09C Unknown Unknown Unknown

doudican_current. 0000000140119BE7 Unknown Unknown Unknown

doudican_current. 0000000140119AEE Unknown Unknown Unknown

kernel32.dll 00000000773D652D Unknown Unknown Unknown

ntdll.dll 0000000077ACC521 Unknown Unknown Unknown

Line 4599 is the Call to Pardiso

Line 404 is the Call from the main program to Doudican_Solve

Declarations to the argumenst to Doudican_solve

REAL(KIND=DBL), ALLOCATABLE, DIMENSION(:) :: BL,SOLUTION

REAL(KIND=DBL), ALLOCATABLE, DIMENSION(:,:) :: GLK

with KIND=DBL defined in a module as follows:

MODULE all

IMPLICIT NONE

SAVE

!Determine processor KIND number for single and double precision

INTEGER, PARAMETER :: SGL = SELECTED_REAL_KIND(p=6,r=37) !p=decimal digits of precision

INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(p=13,r=200) !r=maximum range > 10^r

Civil engineer trying to program using MKL library Pardiso

mecej4 — Sat, 01 Oct 2011 03:02:00 GMT

The traceback establishes that the crash was inside MKL/Pardiso. Usually, this is caused not by an error in Pardiso but by the user passing erroneous values of the index vectors IA, JA. Pardiso does check the input arguments, but it is possible to get erroneous input past those checks. The C0000005 trap is a second line of defence against such errors.

You have to examine whether the arguments to Pardiso satisfy the specifications, either using a debugger or by printing the contents just before the call.

I note at least one error: in the packed column storage convention, the start-of-column index array IA has dimension (NEQ+1). You have not filled in a value for the last element of the array. Here is a simpler version of the code for full matrix to packed column conversion:

[fortran]      KOUNT=COUNT(GLK /= 0d0)
!
      ALLOCATE ( Astar(KOUNT), JA(KOUNT) )
!
      KOUNT=0
      IA(1)=1
      DO J=1,NEQ        ! process columns 1, 2, ...
         IA(J+1)=IA(J)  ! cumulate
         DO I=1,NEQ
              IF (GLK(I,J)/=0) THEN
                  IA(J+1)=IA(J+1)+1
                  KOUNT=KOUNT+1
                  Astar(KOUNT)=GLK(I,J)
                  JA(KOUNT)=I
              END IF
          END DO
      END DO
[/fortran]

How does subroutine Doudican_Solv know the value of NEQ?

Civil engineer trying to program using MKL library Pardiso

Gennady_F_Intel — Sun, 02 Oct 2011 15:46:12 GMT

also, first of all in such cases, we'd recommend to check the input matrix representation by settingparm(27)=1

see the description of it into manual.

--Gennady

Civil engineer trying to program using MKL library Pardiso

Brad_Doudican — Wed, 05 Oct 2011 14:30:24 GMT

Thanks mecej4 - good catch and fix on IA. I'm still debugging (error still hits), but to answer your question NEQ is passed in from Module "all". It is the dimension of the original sparse matrix (NEQ x NEQ).

Civil engineer trying to program using MKL library Pardiso

Brad_Doudican — Wed, 05 Oct 2011 14:33:56 GMT

Hi Gennady - thanks for the insight. Can I just enter the parameter definition as follows, or do I have to spec all 64 components of IPARM?

! DEFINE PARAMETERS

DO I=1,64

PT(I)=0

IPARM(I)=0

END DO

IPARM(27)= 1

MAXFCT =1

MNUM = 1

MTYPE = 2

PHASE = 11

NRHS = 1

MSGLVL = 0

Civil engineer trying to program using MKL library Pardiso

Brad_Doudican — Wed, 05 Oct 2011 14:37:28 GMT

Thanks for all your help. After some more work, I've gotten it this far. The code below bonks out, giving a "Program Exception - Access Violation" at the Call to Pardiso. Any thoughts on why that might be?

SUBROUTINE DOUDICAN_SOLVE (GLK, BL, SOLUTION)

! THIS SUBROUTINE

USE all

USE MKL_PARDISO

implicit none

REAL(KIND=DBL), INTENT(IN), DIMENSION(NEQ, NEQ) :: GLK

REAL(KIND=DBL), INTENT(IN), DIMENSION(NEQ) :: BL

REAL(KIND=DBL), INTENT(OUT), DIMENSION(NEQ) :: SOLUTION

REAL(KIND=DBL), ALLOCATABLE, DIMENSION(:) :: Astar

INTEGER, ALLOCATABLE, DIMENSION(:) :: JA

INTEGER, DIMENSION(NEQ+1) :: IA

!INTEGER, DIMENSION(64) :: PT

TYPE(MKL_PARDISO_HANDLE), ALLOCATABLE :: PT(:)

INTEGER, ALLOCATABLE :: IPARM(:)

INTEGER, DIMENSION(NEQ) :: PERM

INTEGER :: MAXFCT, MNUM, MTYPE, PHASE, NRHS, MSGLVL, ERROR, M, N, &

COUNT, COUNT1, COUNT2, COUNT2OLD, I, J, KOUNT, DUMMY

! MODIFY GLK TO BECOME Astar

KOUNT=COUNT(GLK/=0d0)

ALLOCATE ( Astar(KOUNT), JA(KOUNT) )

KOUNT = 0

IA(1)=1

DO J=1,NEQ

IA(J+1)=IA(J)

DO I=1,NEQ

IF (GLK(I,J)/=0) THEN

IA(J+1)=IA(J+1)+1

KOUNT=KOUNT+1

Astar(KOUNT)=GLK(I,J)

JA(KOUNT)=I

END IF

END DO

! DEFINE PARAMETERS

ALLOCATE ( PT (64) )

ALLOCATE ( IPARM (64) )

DO I=1,64

IPARM(I)=0

PT(I)%DUMMY=0

END DO

IPARM(1) = 0 ! solver default

IPARM(27) = 1

MTYPE = 2

MAXFCT =1

MNUM = 1

NRHS = 1

MSGLVL = 0

ERROR = 0

! CALL THE SPARSE MATRIX SOLVER PARDISO

PHASE = 13

CALL PARDISO_ (PT, MAXFCT, MNUM, MTYPE, PHASE, NEQ, Astar, IA, JA, PERM, NRHS, IPARM, MSGLVL, BL, SOLUTION, ERROR)

RETURN

END SUBROUTINE DOUDICAN_SOLVE