- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello, I am trying to solve a sparse linear system in fortran 90 using the mkl pardiso solver as in the following code:
do i = 1, 64
iparm(i) = 0
end do
error = 0 ! initialize error flag
msglvl = 1 ! print statistical information
mtype = 13 ! complex unsymmetric
!C.. Initiliaze 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
maxfct = 1
mnum = 1
nrhs = 1
!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
print*, ' calling sparse solver'
CALL PARDISO (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns, &
perm, nrhs, iparm, msglvl, b, Ex, error)
WRITE(*,*) 'Reordering completed ... '
My code compiles (intel compiler ifort (IFORT) 12.1.3 20120130) without problem but once reached the solver call the program stops and I have to kill it.
The output on my screen is, before the program behave strangely, just
calling sparse solver
I've also tried with different parameters in iparm, as per the examples folder, but without different results.
I also print out the matrix, to be sure is in the correct form.
double complex values:(-1.000000,0.0000000E+00) (1.000000,0.0000000E+00) (0.5000000,0.0000000E+00)
(-1.500000,0.0000000E+00) (1.000000,0.0000000E+00) (-0.2222222,0.0000000E+00)
(-0.1111111,0.0000000E+00) (0.3333333,0.0000000E+00) (1.000000,0.0000000E+00)
integer rowIndex: 1 3 6 9 10
integer columns: 1 2 1 2 3 2 3 4 4
b is imaginary and the system has a solution.
Any suggestion?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I think that you have errors related to a misunderstanding of the nature of the intrinsic function CMPLX. Because kind(cmplx(1.0d0,1.0d0)) evaluates to 4, whereas kind(cmplx(1.0d0,1.0d0,kind=kind(1d0))) evaluates to 8, the latter of which you probably want, many errors arise. Note, as well, that Pardiso will not work correctly with default precision real and complex variables unless you set iparm(28) = 1. You need REAL(8) and COMPLEX(8) variables if you set the default value of zero in iparm(28).
Try the following program to see the implications:
[fortran]
program chkcmplx
implicit none
complex(8) :: i1
i1=cmplx(0d0,1d-1)
write(*,*)i1
i1=cmplx(0d0,1d-1,kind(0d0))
write(*,*)i1
write(*,*)kind(cmplx(1d0,1d0)),kind(cmplx(1.0d0,1.0d0,kind=kind(1d0)))
end [/fortran]
I have attached a zip file with corrected sources. The program built with these sources worked correctly with IFort 11.1.070 and 14.0.2, on Windows
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
into the code you show - all iparm are zero. Would you please give the example which we can compile and check the problem on our side?
please use iparm(27) = 1 to check the input matrix format.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
It is a bit long, I've put everything here
include 'mkl_pardiso.f90'
module d1isoModel
use kinddef
USE mkl_pardiso
implicit none
type :: d1isoGrid
private
real(rk), dimension(:), allocatable :: z
real(rk), dimension(:), allocatable :: sigma
integer(ik):: n
contains
procedure, public :: init => initd1isoModel
procedure, public :: printOnScreen => printd1isoModelOnScreen
procedure, public :: computeEx
procedure, public :: deallocate => clean_d1isoGrid
end type d1isoGrid
contains
subroutine initd1isoModel(this,n,z,sigma)
class(d1isoGrid) :: this
integer(ik) :: n
real(rk), dimension(n) :: z
real(rk), dimension(n-1) :: sigma
allocate(this%z(n))
allocate(this%sigma(n-1))
this%z = z
this%sigma = sigma
this%n = n
end subroutine initd1isoModel
subroutine clean_d1isoGrid(this)
class(d1isoGrid) :: this
if (allocated(this%sigma))then
deallocate(this%sigma)
end if
if (allocated(this%z))then
deallocate(this%z)
end if
print*, 'memory used for d1isoGrid relased ...'
end subroutine clean_d1isoGrid
function computeEx(this,freq,Ex) result(info)
use physicsparameters
class(d1isoGrid), intent(in) :: this
real(rk), intent(in) :: freq
complex(ck), dimension(:), allocatable :: Ex
integer :: info
! pardiso header file, parameters and variables, based on the intel example
TYPE(MKL_PARDISO_HANDLE), dimension(64) :: pt
INTEGER :: iparm(64)
INTEGER :: i, idum
complex(rk) :: ddum!, waltime1, waltime2
INTEGER :: maxfct, mnum, mtype, phase, nPardiso, nrhs, error, msglvl
! local arguments!
real(rk) :: omega
real(rk), dimension(this%n -3) :: znodes, C
integer :: j
! integer (ik) :: icol, irow
integer :: n
integer :: elementIndex
real(rk) :: zGhostNode
complex (ck), dimension ((this%n -3)*3 + 2 + 1) :: values
integer , dimension ((this%n -3)*3 + 2 + 1) :: columns
integer , dimension (this%n) :: rowIndex
complex (ck), dimension(this%n-1) :: b
integer, dimension(this%n -1) :: perm
n = this%n
nPardiso = n-1
! compute znodes coordiantes
do i = 2,n
znodes(i-1) = (this%z(i) + this%z(i-1))*0.5_rk
end do
omega = 2*pi*freq
allocate (Ex(this%n -1)) ! columns store respectively Ex, Ey and Ez.
b = cmplx(0._rk,0._rk)
b(1) = omega*mu0*i1
C = computeC(n,this%z,znodes)
!print*, znodes
!print*, C
zGhostNode = znodes(1) - (znodes(2) - znodes(1))
! print*, zGhostNode, znodes
print*, 'forming the finite-differences sparse matrix and r.h.s. ...'
values(1) = -((znodes(1) - zGhostNode)**(-1))
columns(1) = 1
rowIndex(1) = 1
values(2) = ((znodes(1) - zGhostNode)**(-1))
columns(2) = 2
elementIndex = 2
!print*, 'n = ', n
do i = 2,n-2
do j = -1,1
elementIndex = elementIndex+1
select case (j)
case(-1)
values(elementIndex) = C(i-1)*(znodes(i+1)-znodes(i))
!print*, c(i), values(elementIndex)
columns(elementIndex) = i+j
rowIndex(i) = elementIndex
case(0)
values(elementIndex) = -C(i-1)*(znodes(i+1)-znodes(i-1))
columns(elementIndex) = i+j
case(1)
values(elementIndex) = C(i-1)*(znodes(i)-znodes(i-1))
columns(elementIndex) = i+j
end select
end do
end do
elementIndex = elementIndex+1
values(elementIndex) = 1._rk
columns(elementIndex) = n-1
rowIndex(n-1) = elementIndex
rowIndex(n) = size(values)+1
print*, values
print*, columns
print*, rowIndex
! calling the sparse solver
!C..
!C.. Set up PARDISO control parameter
!C..
do i = 1, 64
iparm(i) = 0
end do
!!$ iparm(1) = 0 ! no solver default
!!$ iparm(2) = 2 ! uses openMP
!!$ iparm(3) = 0 ! numbers of processors
!!$ iparm(4) = 62 ! 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) = 0 ! 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) = 1 ! not in use
!!$ 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
iparm(27) = 1
error = 0 ! initialize error flag
msglvl = 1 ! print statistical information
mtype = 13 ! complex unsymmetric
!C.. Initiliaze 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
maxfct = 1
mnum = 1
nrhs = 1
!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
print*, ' calling sparse solver'
CALL PARDISO (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns, &
perm, nrhs, iparm, msglvl, b, Ex, 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)
!C.. Factorization.
phase = 22 ! only factorization
CALL pardiso (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns &
, perm, nrhs, iparm, msglvl, b, Ex, error)
WRITE(*,*) 'Factorization completed ... '
IF (error .NE. 0) THEN
WRITE(*,*) 'The following ERROR was detected: ', error
STOP
ENDIF
!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, nPardiso, values, rowIndex, columns &
, perm, nrhs, iparm, msglvl, b, Ex, error)
WRITE(*,*) 'Solve completed ... '
WRITE(*,*) 'The solution of the system is '
DO i = 1, n
WRITE(*,*) ' Ex(',i,') = ', Ex(i)
END DO
!C.. Termination and release of memory
phase = -1 ! release internal memory
CALL pardiso (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns &
, perm, nrhs, iparm, msglvl, b, Ex, error)
info = error
return
end function computeEx
function computeC(n,z,znodes) result(C)
integer(ik), intent(in) :: n
real(rk), dimension(n) :: z
real(rk), dimension(n-1) :: znodes
real(rk), dimension(2:n-2) :: C
integer(ik) :: i
!C = 0._rk
!print*, 'c ='
!print*, C
do i = 2,n-2
C(i) = (z(i+1)-z(i))**(-1) * ((znodes(i+1)-znodes(i))*(znodes(i)-znodes(i-1)))**(-1)
end do
!print*, 'c ='
!print*, C
return
end function computeC
end module d1isoModel
and the supporting modules
module kinddef
implicit none
integer, parameter :: ik = kind(1)
integer, parameter :: rk = kind(1.0d0)
integer, parameter :: ck = kind(cmplx(1.0d0,1.0d0))
save
end module
and
module physicsparameters
use kinddef
implicit none
complex(ck), parameter :: i1 = cmplx(0._rk,1._rk)
real(rk), parameter :: pi = acos(-1._rk)
real(rk), parameter :: mu0 = 4._rk*pi*1.E-7
save
end module physicsparameters
I call everything from main
program hello
use kinddef
use d1isoModel
implicit none
type(d1isoGrid) :: myGrid
real(rk), dimension(4) :: s
real(rk), dimension(5) :: z
complex(ck), dimension(:), allocatable :: E
integer(ik) :: info
z = (/-1., 0., 1., 2., 4./)
s = (/0., 0.001, .001, .1/)
call myGrid%init(5,z,s)
call myGrid%printOnScreen()
info = computeEx(myGrid,3._rk,E)
print*, 'info = ', info
call myGrid%deallocate()
deallocate(E)
end program
I use the flags
ifort -o test *.o -i8 -openmp -I$MKLROOT/include $MKLROOT/lib/libmkl_intel_ilp64.a $MKLROOT/lib/libmkl_core.a $MKLROOT/lib/libmkl_intel_thread.a -lpthread -l
to link
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Incomplete, can't compile because source code of printd1isoModelOnScreen is not provided.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
yes, you are right. you can comment the subroutine call and definition out as it is inessential to the mkl_pardiso problem. It only print sigma and z on screen for testing
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I found a couple of array bounds errors in your program. In Function ComputeEx, you need dimension znodes(this%n-1) rather than (this%n-3).
Secondly, I am not familiar with the F95 interface to Pardiso, so I replaced the F95 style calls with F77 calls. The changes can be found in the attached file. With these changes, your program ran to completion, with the solution
[bash]
Ex( 1 ) = (0.000000000000000E+000,-8.290467793798273E-005)
Ex( 2 ) = (0.000000000000000E+000,-5.921762709855910E-005)
Ex( 3 ) = (0.000000000000000E+000,-3.553057625913546E-005)
Ex( 4 ) = (0.000000000000000E+000,0.000000000000000E+000)
[/bash]
I think that if my assessment is correct, you can simplify the program by avoiding complex numbers. The assessment is based on observing that the matrix is real, and the R.H.S. vector b is purely imaginary. Therefore, instead of solving A x = b = i br, you can solve A y = br, and then take the solution as x = i y. If you can take advantage of this observation, you may pick up considerable improvement in execution speed.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I had some improvements on my system using the provided file.
At the moment I can compile and run the code without interruption if I link the libraries dynamically.
Nevertheless I obtain a sad set of (0., NaN) result.
ericmacbook:d1iso eric$ ./test
forming the finite-differences sparse matrix and r.h.s. ...
(-1.000000,0.0000000E+00) (1.000000,0.0000000E+00) (1.000000,0.0000000E+00)
(-2.000000,0.0000000E+00) (1.000000,0.0000000E+00) (1.000000,0.0000000E+00)
(-1.666667,0.0000000E+00) (0.6666667,0.0000000E+00) (1.000000,0.0000000E+00)
1 2 1 2 3 2
3 4 4
1 3 6 9 10
calling sparse solver
=== PARDISO: solving a complex nonsymetric system ===
The local (internal) PARDISO version is : 103911000
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON
Scaling is turned ON
Matching is turned ON
Summary: ( reordering phase )
================
Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000010 s
Time spent in reordering of the initial matrix (reorder) : 0.000123 s
Time spent in symbolic factorization (symbfct) : 0.000026 s
Time spent in allocation of internal data structures (malloc) : 0.000267 s
Time spent in additional calculations : 0.000017 s
Total time spent : 0.000443 s
Statistics:
===========
< Parallel Direct Factorization with number of processors: > 1
< Numerical Factorization with Level-3 BLAS performance >
< Linear system Ax = b >
number of equations: 4
number of non-zeros in A: 9
number of non-zeros in A (): 56.250000
number of right-hand sides: 1
< Factors L and U >
number of columns for each panel: 128
number of independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
number of supernodes: 2
size of largest supernode: 2
number of non-zeros in L: 10
number of non-zeros in U: 2
number of non-zeros in L+U: 12
Reordering completed ...
Number of nonzeros in factors = 12
Number of factorization MFLOPS = 0
=== PARDISO is running in In-Core mode, because iparam(60)=0 ===
Percentage of computed non-zeros for LL^T factorization
50 % 90 % 100 %
=== PARDISO: solving a complex nonsymetric system ===
Single-level factorization algorithm is turned ON
Summary: ( factorization phase )
================
Times:
======
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct) : 0.001583 s
Time spent in allocation of internal data structures (malloc) : 0.000015 s
Time spent in additional calculations : 0.000002 s
Total time spent : 0.001600 s
Statistics:
===========
< Parallel Direct Factorization with number of processors: > 1
< Numerical Factorization with Level-3 BLAS performance >
< Linear system Ax = b >
number of equations: 4
number of non-zeros in A: 9
number of non-zeros in A (): 56.250000
number of right-hand sides: 1
< Factors L and U >
number of columns for each panel: 128
number of independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
number of supernodes: 2
size of largest supernode: 2
number of non-zeros in L: 10
number of non-zeros in U: 2
number of non-zeros in L+U: 12
gflop for the numerical factorization: 0.000000
gflop/s for the numerical factorization: 0.000051
Factorization completed ...
=== PARDISO: solving a complex nonsymetric system ===
Summary: ( solution phase )
================
Times:
======
Time spent in direct solver at solve step (solve) : 0.000097 s
Time spent in additional calculations : 0.000038 s
Total time spent : 0.000135 s
Statistics:
===========
< Parallel Direct Factorization with number of processors: > 1
< Numerical Factorization with Level-3 BLAS performance >
< Linear system Ax = b >
number of equations: 4
number of non-zeros in A: 9
number of non-zeros in A (): 56.250000
number of right-hand sides: 1
< Factors L and U >
number of columns for each panel: 128
number of independent subgraphs: 0
< Preprocessing with state of the art partitioning metis>
number of supernodes: 2
size of largest supernode: 2
number of non-zeros in L: 10
number of non-zeros in U: 2
number of non-zeros in L+U: 12
gflop for the numerical factorization: 0.000000
gflop/s for the numerical factorization: 0.000051
Solve completed ...
The solution of the system is
Ex( 1 ) = (0.0000000E+00,NaN)
Ex( 2 ) = (0.0000000E+00,NaN)
Ex( 3 ) = (0.0000000E+00,NaN)
Ex( 4 ) = (0.0000000E+00,NaN)
info = 0
memory used for d1isoGrid relased ...
ericmacbook:d1iso eric$
maybe you know any smart debugging flag I can use to individuate/solve the problem?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I think that you have errors related to a misunderstanding of the nature of the intrinsic function CMPLX. Because kind(cmplx(1.0d0,1.0d0)) evaluates to 4, whereas kind(cmplx(1.0d0,1.0d0,kind=kind(1d0))) evaluates to 8, the latter of which you probably want, many errors arise. Note, as well, that Pardiso will not work correctly with default precision real and complex variables unless you set iparm(28) = 1. You need REAL(8) and COMPLEX(8) variables if you set the default value of zero in iparm(28).
Try the following program to see the implications:
[fortran]
program chkcmplx
implicit none
complex(8) :: i1
i1=cmplx(0d0,1d-1)
write(*,*)i1
i1=cmplx(0d0,1d-1,kind(0d0))
write(*,*)i1
write(*,*)kind(cmplx(1d0,1d0)),kind(cmplx(1.0d0,1.0d0,kind=kind(1d0)))
end [/fortran]
I have attached a zip file with corrected sources. The program built with these sources worked correctly with IFort 11.1.070 and 14.0.2, on Windows
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
oh, that was enlightening!
Thanks a million for your kind help.
As per the real vs complex you are right again, but both the matrix and the rhs doom is to become complex, so I will not be able to separate them again. If i'll be possible to solve independently for real and imaginary part, I'll keep your tip in mind.
Thanks again.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page