I am trying to solve a large sparse matrix system using the PARDISO subroutine in Intel MKL. I completed writing a code in gfortran and makefile. But when I execute the makefile, the error message pops up saying "undefined reference to 'mkl_pardiso_' ". I am running the code on Ubuntu 22.04 LTS, and the gfortran code and makefile are below:
program large_sparse_solver
use omp_lib
! Include Intel MKL headers
include 'mkl.fi'
! Define the matrix dimensions
integer, parameter :: eps = 5.0d-15
integer, parameter :: p = 12
integer, parameter :: m = 2*(p+1) + 4*p
integer, parameter :: nb = p*(m-p-1) ! Block size
integer, parameter :: n = (m-p-1)**2 ! Size of the system
integer, parameter :: num_blocks = n/nb ! Number of blocks in the sparse matrix
! Define arrays
real(8), dimension(n,n) :: A
real(8), dimension(n) :: b, x
! Declarations and variables for timing and residual
integer :: i, j, iTimes1, iTimes2, rate
real(8) :: elapsed_time
real(8) :: rsdl(n) ! Residual
! Initialize MKL parameters
integer, parameter :: max_threads = 28 ! Number of OpenMP threads
! Load A, b
OPEN(1, FILE='./data/stifk1n5', STATUS='old', ACCESS='SEQUENTIAL', FORM='FORMATTED', ACTION='READ')
OPEN(2, FILE='./data/rhsk1n5', STATUS='old', ACCESS='SEQUENTIAL', FORM='FORMATTED', ACTION='READ')
do i = 1, n
read(1, *) (A(i,j), j=1,n)
read(2, *) b(i)
enddo
close(1); close(2)
! Initialize x
x(:) = 0.0d0
call system_clock(count_rate=rate)
! Start timing
call system_clock(iTimes1)
! Set the number of threads for OpenMP
call omp_set_num_threads(max_threads)
! Initialize MKL
call mkl_set_num_threads(max_threads)
! Use the PARDISO solver from MKL for solving
call solve_sparse_system(A, b, x, n)
! Stop timing
call system_clock(iTimes2)
! Calculate and print the execution time
elapsed_time = real(iTimes2-iTimes1)/real(rate)
print *, "Execution time (seconds): ", elapsed_time
! Clean up
call mkl_free_buffers()
contains
! PARDISO in Intel MKL for Fortran
! Calculates the solution of a set of sparse linear equations with single or multiple right-hand sides.
! https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2023-0/pardiso.html#GUID-431916D5-B76D-48A1-ABB5-1A0613FDC0FA
subroutine solve_sparse_system(A, b, x, n)
integer, intent(in) :: n
real(8), intent(inout) :: A(n,n), b(n)
real(8), intent(out) :: x(n)
! Input parameters of PARDISO
integer :: mtype, nrhs, maxfct, mnum, phase
! Output parameters of PARDISO
integer :: perm(n), iparm(64), error, msglvl
integer :: i, j, k, nz, ia(n+1)
integer, allocatable :: ja(:)
real(8), allocatable :: nza(:)
TYPE(mkl_pardiso_handle) :: pt(64)
! Convert the matrix system into three array variation of CSR format
call get_CSR(A, b, n, nza, ia, ja, nz)
! Initialize PARDISO solver parameters
! iparm(:) = 0
! pt(:) = 0
mtype = 11
maxfct = 1
mnum = 1
nrhs = 1
phase = 23
msglvl = 0
call pardisoinit (pt, mtype, iparm)
! iparm(2) = 2 ! Parallel factorization
! iparm(3) = 3 ! Parallel solve
! iparm(4) = 0 ! No iterative refinement
! iparm(8) = 2 ! Use Intel MKL PARDISO
! Solve the matrix system using PARDISO
call mkl_pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, perm, nrhs, iparm, msglvl, b, x, error)
if (error /= 0) then
print *, "PARDISO solver returned error code: ", error
stop
end if
! Deallocate memory
! deallocate(ia, ja, nza)
end subroutine solve_sparse_system
! Convert the matrix system into three array variation of CSR format
! Output :
! nz = the number of non-zero elements in A(n,n)
! nza = the arrya of non-zero elements of A, (nz by 1) array
! ia(i) = the index of element in 'nza' corresponding to the index of column of the first non-zero element in the i-th row of A (rowIndex), (n+1 by 1) array, ia(n+1) = nz + 1
! ja(j) = the index of column in A corresponding to the index of element in 'nza', (nz by 1) array
subroutine get_CSR(A, b, n, nza, ia, ja, nz)
integer, intent(in) :: n
real*8, intent(in) :: A(n,n), b(n)
integer, intent(inout) :: nz, ia(n+1)
integer, allocatable, intent(inout) :: ja(:)
real*8, allocatable, intent(inout) :: nza(:)
logical :: non_mask(n,n)
integer :: i, j, k, ii, jj, kk
! Count the number of non-zero elements in matrix A
non_mask = A(:,:) < -eps .or. A(:,:) > eps
nz = count(non_mask)
! Allocate the size of ja, nza
allocate(ja(nz), nza(nz))
! k = 0
! do i = 1, n
! do j = 1, n
! if (dabs(A(i,j))>eps) then
! k = k + 1
! endif
! enddo
! enddo
k = 0
do i = 1, n
kk = 0
do j = 1, n
if (dabs(A(i,j))>eps) then
k = k + 1
nza(k) = A(i,j)
ja(k) = j
if (kk==0) then
ia(i) = k
kk = 1
endif
endif
enddo
enddo
ia(n+1) = nz + 1
end subroutine get_CSR
end program large_sparse_solver
Here is the makefile.
# Makefile for compiling and running the large_sparse_solver Fortran program
# Compiler settings
FC = gfortran
FLAGS = -O3 -fopenmp -ffixed-line-length-none -ffree-line-length-none
# Libraries (Intel MKL)
MKLROOT = /opt/intel/oneapi/mkl/2023.2.0
MKL_LIB = -L$(MKLROOT)/lib/intel64 -lmkl_intel_lp64 -lmkl_gnu_thread -lmkl_core -lgomp -lpthread -lm -ldl
# MKL_LIB = -L$(MKLROOT)/lib/intel64 -libmkl_intel_lp64 -libmkl_gnu_thread -libmkl_core -libgomp -lpthread -lm -ldl
# MKL_PARDISO = -L$(MKLROOT)/include/mkl_pardiso.f90
# Include directory for MKL headers
MKL_INCLUDE = -I$(MKLROOT)/include
# Source files
SRC = large_sparse_solver.f90
# Executable name
EXEC = test_run
# Targets
all: $(EXEC)
$(EXEC): $(SRC)
$(FC) $(FLAGS) $(MKL_INCLUDE) $^ -o $@ $(MKL_LIB)
clean:
rm -f *.mod $(EXEC)
I need a makefile to execute gfortran files calling Intel MKL subroutines.