- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Brian Murphy wrote:

Something has changed, but I don't know what. I can now putuse mkl_pardisoin a subprogram unit, and that will compile without errors. I was getting errors before, but not now. I do not haveinclude '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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page