Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Brad_Doudican
Beginner
126 Views

Civil engineer trying to program using MKL library Pardiso

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
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 IF
END DO
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
!
!***********************************************************************
0 Kudos
7 Replies
mecej4
Black Belt
126 Views

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.
Brad_Doudican
Beginner
126 Views

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
mecej4
Black Belt
126 Views

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?
Gennady_F_Intel
Moderator
126 Views

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
Brad_Doudican
Beginner
126 Views

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).
Brad_Doudican
Beginner
126 Views

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
Brad_Doudican
Beginner
126 Views

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
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
Reply