how to use pardiso in f90 program

Brian_Murphy · ‎03-01-2018

What is the right way to "use" or "include" mkl_pardiso in an f90 program?
use mkl_pardiso generates a compile about not being to find the file in the include directory, but it can find blas95 and lapack95 without any trouble.

Brian_Murphy · ‎03-01-2018

I got use mkl_pardiso to work, but I'm not exactly sure what fixed it. I put this use statement inside a module of my own, right after Use Blas95 and Use Lapack95, and that eliminated the compiler error. I then put "use" for my own module where I want to call pardiso, and that's working. I don't understand why I couldn't put use mkl_pardiso in my program where I'm actually calling pardiso.

Julian_H_ · ‎03-01-2018

There is an example in your fortran installation folder

..\Program Files (x86)\Intel\Composer XE 2015\mkl\examples\examples_core_f.zip\solverf\source\pardiso_sym_f90.f90

Basically you include

INCLUDE 'mkl_pardiso.f90'

and then you can use

USE mkl_pardiso

Also make sure you have acitvated the mkl library in your compiler settings

mecej4 · ‎03-02-2018

The source file mkl_pardiso.f90 contains the definition of a complete module. Therefore, the INCLUDE statement that was shown in #3 should be situated outside the body of any other subprogram or main program. In fact, since there is no need to generate the module file mkl_pardiso.mod more than once, you should simply compile mkl_pardiso.f90 separately, and place the resulting module file in an accessible place. After that has been done, it should suffice to USE the module, and no INCLUDE needs to be done.

This is even more notable in the case of modules such as BLAS95 and Lapack95, which take quite a bit of time to compile. Generate the module files once, place them properly, and USE the generated modules as needed.

Brian_Murphy · ‎03-02-2018

Something has changed, but I don't know what. I can now put use mkl_pardiso in a subprogram unit, and that will compile without errors. I was getting errors before, but not now. I do not have include 'mkl_pardiso.90' anywhere. So anyhow, I'm past that problem, and on the next one regarding pardiso.

I've created a defined type for a csr format sparse matrix to use with pardiso.

 Type csrMatrix
  real(8), allocatable :: a(:)
  integer, allocatable :: ia(:), ja(:)
  integer :: nx, ne
  type(MKL_PARDISO_HANDLE)  :: pt(64)
  integer :: iparm(64)
 end type csrMatrix

Are there any obvious problems with this defined Type? In particular, are allocatable arrays ok?

I am using this with pardiso for eigenvalue calculations with real general asymmetric matrices (i.e. ARPACK). I allocate the arrays, load the data, and the calculation of eigenvalues goes fine. But afterwards when I deallocate the arrays in the csrMatrix object and call pardiso with phase=-1 to "Release all internal memory for all matrices", something goes wrong and other allocated arrays are mysteriously becoming deallocated. Anyhow, I'm hoping this is a simple problem with how I've setup the Type csrMatrix.

Brian_Murphy · ‎03-02-2018

I think I solved my lost allocation problem. When calling pardiso with phase=-1 I was passing a dummy array for iparm, and it evidently needs the actual iparm array. This change got the code running without any more trouble. But lots more testing to do.

mecej4 · ‎03-03-2018

Brian Murphy wrote:
Something has changed, but I don't know what. I can now put use mkl_pardiso in a subprogram unit, and that will compile without errors. I was getting errors before, but not now. I do not have include 'mkl_pardiso.90' anywhere.

There is no mystery here. You inserted include 'mkl_pardiso.f90' in one of your source files and compiled that source file at least once. That compilation generated the mkl_pardiso.mod file in an accessible place in the file system. During subsequent compilations of sources containing use mkl_pardiso the compiler found that module file and used the definitions in the module.

I commented on this matter in #4.

Brian_Murphy · ‎03-03-2018

I think I follow you. My Composer mkl\include\ia32 folder contains .mod files for blas and lapack, but not for pardiso. I have been compiling in Debug mode from visual studio, and I indeed see mkl_pardiso.mod in my project's \Debug folder. So that confirms I did what you described. Do you recommend I copy this file to the MKL\include\ia32 folder?

I just switched to Release mode, and after covering the bases with the include files and visual studio project settings, I got that to build successfully.

The Debug and Release versions of mkl_pardiso.mod appear to be the same.

mecej4 · ‎03-03-2018

Brian Murphy wrote:
... So that confirms I did what you described. Do you recommend I copy this file to the MKL\include\ia32 folder?

Sure, if it was generated during the course of building a 32-bit EXE, DLL or LIB. Doing so will make that file available for other projects where that module will be used.

JohnNichols · ‎03-03-2018

This is the Pardiso solver routine that mecej4 "helped" develop over the last few years. Actually he has all the credit all I did was ask the questions and then follow his answers.

It has the advantage you can pass in a square matrix and it does everything. You just need to set the lib in the properties.

I have actually been looking at in the last few weeks as we grapple with billions of vibration results from all over the world. I do cross checks on the results occasionally referring to the original Gaussian solver provided by Harrison in 1972 and coded in VOF - very old Fortran -- there are error differences of about 1 - 2% on the solution - which is what one expects, but it is interesting to compare the results from published data which compares numerical results using MATLAB to theoretical problems and says the results are within 0.5%.

My students all want to use MATLAB and I can say with a passion - it is destroying engineering and programming.

      !---------------------------------------------------------------------------------------------------------------------------
      !
      !
      !
      !---------------------------------------------------------------------------------------------------------------------------


!===============================================================================
! Copyright 2004-2016 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

!===============================================================================
! Copyright 2004-2016 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

!   Content : MKL PARDISO Fortran example
!
!********************************************************************************
!C----------------------------------------------------------------------
!C Example program to show the use of the "PARDISO" routine
!C for symmetric linear systems
!C---------------------------------------------------------------------
!C This program can be downloaded from the following site:
!C www.pardiso-project.org
!C
!C (C) Olaf Schenk, Department of Computer Science,
!C University of Basel, Switzerland.
!C Email: olaf.schenk@unibas.ch
!C
!C---------------------------------------------------------------------

      MODULE Solver
      use Base


      contains

      subroutine linsolve(n, error, nnz, a, ja, ia, b, x)
      implicit none
      include 'mkl_pardiso.fi'

      INTEGER :: n, nnz
      REAL (KIND=dp) :: x(n) ! Pardiso needs X(N) even if sol is returned in b

      REAL (KIND=dp) :: b(n)
!C.. Internal solver memory pointer for 64-bit architectures
!C.. INTEGER*8 pt(64)
!C.. Internal solver memory pointer for 32-bit architectures
!C.. INTEGER*4 pt(64)

!C.. This is OK in both cases
      TYPE(MKL_PARDISO_HANDLE) pt(64)
!C.. All other variables
      INTEGER maxfct, mnum, mtype, phase, nrhs, msglvl
      INTEGER iparm(64)
      INTEGER ia(nnz)
      INTEGER ja(nnz)
      INTEGER i, idum(1)
      REAL (KIND=dp) ::  ddum(1)
      integer error
      REAL (KIND=dp) a(nnz)

!C.. Fill all arrays containing matrix data.
      DATA nrhs /1/, maxfct /1/, mnum /1/
!C..
!C.. Setup PARDISO control parameter
!C..
      DO i = 1, 64
          iparm(i) = 0
      END DO
      iparm(1) = 1 ! no solver default
      iparm(2) = 0 ! 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 components 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 ! perturb the pivot elements with 1E-13
      iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
      iparm(12) = 0 ! not in use
      iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric)
      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 = 0 ! print statistical information
      mtype = 11 ! real unsymmetric


!C.. Initialize the internal solver memory pointer. This is only
!C necessary for the FIRST call of the PARDISO solver.
      DO i = 1, 64
          pt(i)%DUMMY = 0
      END DO


!C.. Reordering and Symbolic Factorization, This step also allocates
!C   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 ... '
      WRITE(sw,*) 'Reordering completed ... '


      IF (error .NE. 0) THEN
          WRITE(*,*) 'The following ERROR was detected: ', error
          WRITE(sw,*) 'The following ERROR was detected: ', error
          STOP 1
      END IF


     !WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
    !  WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
      WRITE(sw,*) 'Number of nonzeros in factors = ',iparm(18)
      WRITE(sw,*) 'Number of factorization MFLOPS = ',iparm(19)


!C.. 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 ... '
      WRITE(sw,*) 'Factorization completed ... '
      IF (error .NE. 0) THEN
          WRITE(*,*) 'The following ERROR was detected: ', error
          WRITE(sw,*) 'The following ERROR was detected: ', error
          STOP 1
      END IF

!C.. Back substitution and iterative refinement
      iparm(8) = 2 ! max numbers of iterative refinement steps
      phase = 33 ! only factorization

      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)
 !     END DO
!C.. 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)
      Return
      End subroutine linsolve

      !---------------------------------------------------------------------------------------------------------------------------
      !
      !
      !
      !---------------------------------------------------------------------------------------------------------------------------

      subroutine PardisoSolver(GK,nl,nk,ndf,X,Load)


      implicit none
      include 'mkl_pardiso.fi'

      integer nl,nk,ndf

      REAL (KIND=dp) GK(nl,nk),X(ndf),Load(ndf),GKTest(nl,nk)

      REAL (KIND=dp), ALLOCATABLE :: a(:), b(:), c(:,:), atemp(:)
      integer, allocatable :: ja(:), ia(:), iatemp(:), jatemp(:)
      INTEGER :: nnodes = 5
      integer error, iaN
      integer :: nnz = 1000, istat, i, j ! count
      integer :: size = 0
      integer :: size1 = 0
      integer :: size2 = 0
      integer :: flag = 0
      REAL (KIND=dp) :: delta = 0.000000000000001
      integer no_zero_count(nl), numA
    
      GKTest = 0.0d0
      no_zero_count = count( gk .NE. GKTEST, DIM=1 )  !  count values in a col 
      numA = sum(no_zero_count)
      !write(*,2040)numA
      write(sa,2040)numA
2040  Format('--------------------------------------------------------------------------------------------------------------',///&
            '           Count of the non-zero elements in the stiffness matrix :: ',i6,//)
       
      nnz =  numA+10     !90*ndf
      iaN = ndf+1

      
    ALLOCATE (a(nnz), ja(nnz), ia(iaN), jatemp(nnz), iatemp(iaN), b(ndf), c(ndf,ndf), atemp(nnz), STAT=istat)
    IF (istat.NE.0) THEN
        WRITE (*, *) '*** Could not allocate some arrays in LINSOLVE'
        STOP
    END IF

      b = 0.0D0
      a = 0.0d0
      c = 0.0d0
      iatemp = 0
      jatemp = 0
      atemp = 0
      ia = 0
      ja = 0
      size = 0
      size1 = 0
      size2 = 0
      do 200 i = 1, ndf
          flag = 0
          do 300 j = 1, ndf
              c(i,j) = gk(i,j)
              if(abs(c(i,j)) .gt. delta) then
                  size = size + 1
                  if(size .gt. nnz) then
                    Write(*,2030)nnz, ndf
                    Write(sa,2030)nnz, ndf
2030                Format( "        NNZ is set at :: ",I6, "    NDF is :: ",I6)                    
                    stop    '        Insufficient memory allocation in PARDISO to solve the inversion problem'
                endif
                  jatemp(size) = j
                  iatemp(size1 + 1) = size+1
                  atemp(size) = c(i,j)
                  if(flag .eq. 0) then
                      size1 = size1 + 1
                      iatemp(size1) = size
                      flag = 1
                  endif
              endif
300   end do
400   format(5(1x,F10.3))
200   end do

      do 500 i = 1, ndf
          b(i) = load(i)
500   end  do
100   Format(' Matrix Count :: ',i5, ' Size :: ',i5)


      call LinSolve(ndf, error,nnz, atemp, jatemp, iatemp, b, X)

      write(sa,2020)
2020  Format(    //'--------------------------------------------------------------',//,'      Solution Vector for the Structures Problem'//,'               i        X(i)',//)
      do 2000 i = 1,NDF
             write(sa,2010)i, x(i)
2000  End do
2010  Format( 12x,i4,5x, e10.3)

      return
      end subroutine PardisoSolver

      end module

Brian_Murphy · ‎03-03-2018

I've been using pardiso for a little over 2 years, and have been most pleased with it. Before that I was using UMFPACK. I concur with what you say about matlab. In the US, fortran is rarely required in engineering curriculums (maybe nowhere). It's interesting that pardiso has its roots in Basel. In June I will be visiting a turbomachinery company located in Basel.