[Fortran] 纯文本查看 复制代码
subroutine RealPardiso(matrix,b,x,N,error)
use libpardiso600
implicit none
integer :: n !// depend on your equation
integer :: i, j, mm, tmp, nn, fileid, first, num, ncolumn, count
real(kind=8) :: matrix(n,n), b(n), x(n),dparm(64)
real(kind=8), allocatable :: aa(:)
integer, allocatable :: ia(:), ja(:)
!// =====================the variables of pardiso====================
integer(kind=8) :: pt(64) !// kind=8:x64; kind=4:win32
integer :: maxfct, mtype, phase, nrhs, error, msglvl, mnum, solver
integer, allocatable :: perm(:)
integer :: iparm(64)
integer, external :: mkl_get_max_threads
!// =====================the variables of pardiso====================
!// input data
allocate( perm( n ) )
!// =========================store data by CSR format=======================
num=0
do I=1,N
do J=1,N
if ( abs(matrix(i,j)) /= 0.0d0 ) then
num=num+1
endif
enddo
enddo
allocate( aa(num), ja(num), ia(n+1) ) !// store the no-zero element
ia(1)=1
count=1
do I=1,N
ncolumn=0
do J=1,N
if ( abs(matrix(I,J)) /= 0.0d0 ) then
ncolumn=ncolumn+1
ja(count)=J
aa(count)=abs(matrix(I,J))
count=count+1
endif
enddo
ia(I+1)=count
enddo
!// ===========================solve equations=============================
pt = 0 !// pointer initialization
maxfct = 1; mnum = 1; mtype = 11 !// mtype = -2: symmetric nonpositive definite matrix, mtype = 2: symmetric positive definite matrix
perm = 0; nrhs = 1
x = 0.d0
iparm(1) = 0 !// iparm use default values
iparm(3) = mkl_get_max_threads() !// iparm(3): parallel threads are obtained by MKL. In general, the number of cpus is the same as the number of threads
error = 0
msglvl = 1
solver=0
! iparm(32)=1
! iparm(33)=1
phase = 12 !// LU decompose
call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error,dparm )
! write(*,*)dparm(32),dparm(33),iparm(33)
iparm(8)=100
iparm(32)=1
! iparm(33)=1
solver=0
phase = 33 !// solve equations
call pardiso( pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error,dparm )
write(*,*)iparm(7),iparm(20),dparm(35)
do I=1,64
write(*,*)dparm(I)
enddo
phase = -1
call pardiso (pt, maxfct, mnum, mtype, phase, n, aa, ia, ja, perm, nrhs, iparm, msglvl, b, x, error,dparm)
deallocate( aa, ja, ia, perm )
!// ===========================stop solving equations=============================
end