topic Re: calling PARDISO from fortran 90 in Intel® oneAPI Math Kernel Library

calling PARDISO from fortran 90

susannebuiter — Fri, 21 Aug 2009 08:30:30 GMT

Hi all,

I have seen some related issues in the discussion forum, but could not find help for my problem, so i hope someone can point me to a good place or help me out!

I am trying to call PARDISO from fortran 90. I first rewrote the example pardiso_sym_f.f to .f90. I run into problems with the type declarations for several of the variables:

ifort -O2 -check all -L /libmkl_solver_sequential.a -Wl,--start-group -lmkl_intel -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -c pardiso_sym_f.f90
pardiso_sym_f.f90(101): error #6633: The type of the actual argument differs from the type of the dummy argument. [PT]
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
-------------------^
pardiso_sym_f.f90(101): error #7836: If the actual argument is scalar, the corresponding dummy argument shall be scalar unless the actual argument is an element of an array that is not an assumed-shape or pointer array, or a substring of such an element. [PERM]
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
-----------^
pardiso_sym_f.f90(101): error #7836: If the actual argument is scalar, the corresponding dummy argument shall be scalar unless the actual argument is an element of an array that is not an assumed-shape or pointer array, or a substring of such an element.
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
-----------^
pardiso_sym_f.f90(101): error #7836: If the actual argument is scalar, the corresponding dummy argument shall be scalar unless the actual argument is an element of an array that is not an assumed-shape or pointer array, or a substring of such an element.
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
-----------^

etc

The code follows here:

INCLUDE 'mkl_pardiso.f90' ! Include the standard pardiso "header file."
PROGRAM pardiso_sym
!
use mkl_pardiso
IMPLICIT NONE
INTEGER, PARAMETER :: dp = KIND(1.0D0)
!.. Internal solver memory pointer
INTEGER,allocatable::pt(:)
!.. All other variables
INTEGER:: maxfct, mnum, mtype, phase, n, nrhs, error, msglvl
INTEGER,allocatable::iparm(:), ia(:),ja(:)
REAL(kind=dp),allocatable::a(:),b(:),x(:)
INTEGER:: i, idum
real(kind=dp):: waltime1, waltime2, ddum
allocate(pt(64))
allocate(iparm(64),ia(9),ja(18))
allocate(a(18),b(8),x(8))
!
!.. Fill all arrays containing matrix data.
n = 8
nrhs = 1
maxfct = 1
mnum = 1
ia = (/1,5,8,10,12,15,17,18,19/)
ja = (/1,3,6,7,2,3,5,3,8,4,7,5,6,7,6,8,7,8/)
a = (/7.d0,1.d0,2.d0,7.d0,-4.d0,8.d0,2.d0,1.d0,5.d0, &
7.d0,9.d0,5.d0,1.d0,5.d0,-1.d0,5.d0,11.d0,5.d0/)
!..
!.. Set up PARDISO control parameter
!..
iparm = 0
iparm(1) = 1 ! no solver default
iparm(2) = 2 ! fill-in reordering from METIS
iparm(3) = 1 ! numbers of processors
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n compoments of x
iparm(7) = 0 ! not in use
iparm(8) = 9 ! numbers of iterative refinement steps
iparm(9) = 0 ! not in use
iparm(10) = 13 ! perturbe the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(12) = 0 ! not in use
iparm(13) = 0 ! maximum weighted matching algorithm is switched-off
!(default for symmetric).
!Try iparm(13) = 1 in case of inappropriate accuracy
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(15) = 0 ! not in use
iparm(16) = 0 ! not in use
iparm(17) = 0 ! not in use
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
!
error = 0 ! initialize error flag
msglvl = 1 ! print statistical information
mtype = -2 ! symmetric, indefinite
!
!.. Initialize the internal solver memory pointer. This is only
! necessary for the FIRST call of the PARDISO solver.
pt = 0
!
!.. Reordering and Symbolic Factorization, This step also allocates
! all memory that is necessary for the factorization
phase = 11 ! only reordering and symbolic factorization
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
idum, nrhs, iparm, msglvl, ddum, ddum, error)
WRITE(*,*) 'Reordering completed ... '
IF (error .NE. 0) THEN
WRITE(*,*) 'The following ERROR was detected: ', error
STOP
END IF
WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
!
!.. Factorization.
phase = 22 ! only factorization
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
idum, nrhs, iparm, msglvl, ddum, ddum, error)
WRITE(*,*) 'Factorization completed ... '
IF (error .NE. 0) THEN
WRITE(*,*) 'The following ERROR was detected: ', error
STOP
ENDIF
!
!.. Back substitution and iterative refinement
iparm(8) = 2 ! max numbers of iterative refinement steps
phase = 33 ! only factorization
do i = 1, n
b(i) = 1.d0
enddo
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
idum, nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Solve completed ... '
WRITE(*,*) 'The solution of the system is '
DO i = 1, n
WRITE(*,*) ' x(',i,') = ', x(i)
ENDDO
!
!.. Termination and release of memory
phase = -1 ! release internal memory
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, &
idum,idum, nrhs, iparm, msglvl, ddum, ddum, error)
!
deallocate(iparm,ia,ja,a,b,x,pt)
END PROGRAM pardiso_sym

Thanks!

Re: calling PARDISO from fortran 90

ArturGuzik — Sat, 22 Aug 2009 04:21:47 GMT

Quoting - susannebuiter

call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &

The code follows here:

!.. Internal solver memory pointer
INTEGER,allocatable::pt(:)

Compiler is pretty clear on what it doesn't like. You declared pt as allocatable array of integers (which platform are you on? You need observe as it is integer*4 or integer*8), while (after checking the subroutine interface, defined in Include file) it expects TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*).
You need to fix this (I didn't check more than that, so other problems may be down the road) as well as argument corresponding to perm (expected to be array) while you're passing idum (scalar, integer).

A.

Note: if you USE mkl_pardiso then you don't need include.

Re: calling PARDISO from fortran 90

Gennady_F_Intel — Mon, 24 Aug 2009 06:25:50 GMT

Susanne,

It seems we need to add the example about calling pardiso from fortran 90.

In general , your program should be looks like ( this is just example)

PROGRAM PardisoTest

use MKL_PARDISO

IMPLICIT NONE

TYPE(MKL_PARDISO_HANDLE), pointer :: PT(:) => null()
integer, pointer :: idum(:) => null(), ivect(:) => null(), jvect(:) => null()
real(KIND=8), pointer :: ddum(:) => null(), a(:) => null()
INTEGER maxfct, mnum, mtype, phase, nrhs, error, msglvl, i, nnz, nequations
INTEGER iparm(64)

DATA nrhs /1/, maxfct /1/, mnum /1/

nnz = XXX
nequations = XXX

do i = 1, 64
iparm(i) = 0
end do
iparm(1) = 0 ! solver default
error = 0 ! initialize error flag
msglvl = 1 ! print statistical information
mtype = 11 ! real unsymmetric --> as an example

allocate( PT(64), idum(0), ddum(0), a(nnz), ivect(nequations+1), jvect(nnz) )

do i = 1, 64
PT(i)%DUMMY = 0
end do

open (unit = 138, file = " ..... ", .... )
do i = 1, nnz
read (138) a(i)
end do

do i = 1, nnz
read (138) jvect(i)
end do

do i = 1, (nequations+1)
read (138) ivect(i)
end do
close(138)

phase = 11

CALL pardiso (PT, maxfct, mnum, mtype, phase, nequations, a, ivect, jvect, idum, nrhs, iparm, msglvl, ddum, ddum, error)

! check the error

END PROGRAM PardisoTest

Regards, Gennady

Re: calling PARDISO from fortran 90

Gennady_F_Intel — Mon, 24 Aug 2009 06:33:45 GMT

added the MVSC 2005 project with the example described above. Probably it will help you. Please pay attention the project contains F90 header file "mkl_pardiso.f90" included into the project.
--Gennady

Re: calling PARDISO from fortran 90

vchan1186 — Fri, 23 Oct 2009 21:03:09 GMT

Quoting - Gennady Fedorov (Intel)

Hi Gennady,

I'm trying to use the MKL pardiso solver. I am able to compile the program, however when I run the program I run into an "abort trap" message.

I downloaded the sample code from your code and the only changes I made was changing the real(kind=8) pointers to complex and instead of reading in a file for a,ivect, and jvect, I created them with another program.

Do you have any suggestions as to why I'm running into such a problem. Are there certain test I can do to pinpoint the problem?

Re: calling PARDISO from fortran 90

Gennady_F_Intel — Sat, 24 Oct 2009 05:45:40 GMT

Quoting - vchan1186

Quoting - Gennady Fedorov (Intel)

added the MVSC 2005 project with the example described above. Probably it will help you. Please pay attention the project contains F90 header file "mkl_pardiso.f90" included into the project.
--Gennady

Hi Gennady,

I'm trying to use the MKL pardiso solver. I am able to compile the program, however when I run the program I run into an "abort trap" message.

I downloaded the sample code from your code and the only changes I made was changing the real(kind=8) pointers to complex and instead of reading in a file for a,ivect, and jvect, I created them with another program.

Do you have any suggestions as to why I'm running into such a problem. Are there certain test I can do to pinpoint the problem?

No idea at this moment what's going wrong with the code.
Can you enable matrix checker? iparm(27)=1.
What phase of calculation do you see the crash?
--Gennady