7 include
'mkl_cluster_sparse_solver.f90'
16 use mkl_cluster_sparse_solver
24 logical,
save :: INITIALIZED = .false.
26 type(MKL_CLUSTER_SPARSE_SOLVER_HANDLE) :: pt(64)
28 integer maxfct, mnum, mtype, nrhs, msglvl
30 data nrhs /1/, maxfct /1/, mnum /1/
32 integer,
save :: iparm(64)
33 integer(kind=kint),
save :: nn
34 integer(kind=kint),
pointer,
save :: irow(:), jcol(:)
35 real(kind=
kreal),
pointer,
save :: aval(:), rhs(:), solx(:)
37 integer(kind=kint),
parameter :: debug=0
43 type (sparse_matrix),
intent(inout) :: spmat
44 integer(kind=kint),
intent(in) :: phase_start
45 integer(kind=kint),
intent(out) :: istat
46 real(kind=
kreal),
pointer,
intent(inout) :: solution(:)
48 integer(kind=kint) :: myrank, phase
49 real(kind=
kreal) :: t1,t2,t3,t4,t5
55 if (.not. initialized)
then
76 if( phase_start == 1 )
then
78 call cluster_sparse_solver ( &
79 pt, maxfct, mnum, mtype, phase, nn, aval, irow, jcol, &
82 write(*,*)
'ERROR: MKL returned with error in phase -1', istat
91 if ( phase_start == 1 )
then
103 iparm(41) = spmat%DISPLS(myrank+1)+1
104 iparm(42) = spmat%DISPLS(myrank+1)+spmat%N_COUNTS(myrank+1)
105 call sort_column_ascending_order(spmat,myrank)
115 call export_spmat_to_centralizedcrs(spmat,myrank,nn,irow,jcol,aval,rhs,solx)
119 if (myrank==0 .and. spmat%timelog > 0) &
120 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: Additional Setup completed. time(sec)=',t2-t1
122 if ( phase_start == 1 )
then
127 call cluster_sparse_solver ( &
128 pt, maxfct, mnum, mtype, phase, nn, aval, irow, jcol, &
131 if (myrank==0 .and. spmat%timelog > 0)
call print_iparm_paramters()
132 write(*,*)
'ERROR: MKL returned with error in phase 1', istat
137 if (myrank==0 .and. spmat%timelog > 0) &
138 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: Analysis completed. time(sec)=',t3-t2
142 if ( phase_start .le. 2 )
then
146 call cluster_sparse_solver ( &
147 pt, maxfct, mnum, mtype, phase, nn, aval, irow, jcol, &
150 if (myrank==0 .and. spmat%timelog > 0)
call print_iparm_paramters()
151 write(*,*)
'ERROR: MKL returned with error in phase 2', istat
155 if (myrank==0 .and. spmat%timelog > 0) &
156 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: Factorization completed. time(sec)=',t4-t3
163 call cluster_sparse_solver ( &
164 pt, maxfct, mnum, mtype, phase, nn, aval, irow, jcol, &
167 if (myrank==0 .and. spmat%timelog > 0)
call print_iparm_paramters()
168 write(*,*)
'ERROR: MKL returned with error in phase 3', istat
175 solution(i) = spmat%rhs(i)
183 if (myrank==0 .and. spmat%timelog > 0)
then
184 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: Solution completed. time(sec)=',t5-t4
186 if( debug>0 .and. myrank==0 )
call print_iparm_paramters()
189 stop
"MKL Pardiso not available"
194 subroutine print_iparm_paramters()
195 write(*,
'(A60,I8)')
'Number of iterative refinement steps performed: ',iparm(7)
196 write(*,
'(A60,I8)')
'Number of perturbed pivots: ',iparm(14)
197 write(*,
'(A60,I8)')
'Peak memory on symbolic factorization: ',iparm(15)
198 write(*,
'(A60,I8)')
'Permanent memory on symbolic factorization: ',iparm(16)
199 write(*,
'(A60,I8)')
'Size of factors/Peak memory on num. fact. and sol: ',iparm(17)
200 write(*,
'(A60,I8)')
'The number of non-zero elements in the factors: ',iparm(18)
201 if( mtype < 11 )
then
202 write(*,
'(A60,I8)')
'Number of positive eigenvalues: ',iparm(22)
203 write(*,
'(A60,I8)')
'Number of negative eigenvalues: ',iparm(23)
207 subroutine export_spmat_to_centralizedcrs(spMAT,myrank,n,ia,ja,a,b,x)
208 type (sparse_matrix),
intent(in) :: spmat
209 integer(kind=kint),
intent(in) :: myrank
210 integer(kind=kint),
intent(out) :: n
211 integer(kind=kint),
pointer,
intent(inout) :: ia(:), ja(:)
212 real(kind=
kreal),
pointer,
intent(inout) :: a(:), b(:), x(:)
214 integer(kind=kint) :: i, nprocs, ierr, nnz
215 integer(kind=kint),
allocatable :: dispmat(:), ncounts(:)
216 real(kind=
kreal) :: t1,t2,t3,t4
221 allocate(dispmat(nprocs),ncounts(nprocs))
229 dispmat(i+1)=dispmat(i)+ncounts(i)
230 nnz = nnz + ncounts(i)
232 nnz = nnz + ncounts(nprocs)
234 if( myrank == 0 )
then
236 if( .not.
associated(ia) )
allocate(ia(nnz+n), stat=ierr)
238 write(*,*)
" Allocation error in Setup Cluster MKL, ia",ierr,nnz
241 if( .not.
associated(ja) )
allocate(ja(nnz+n), stat=ierr)
243 write(*,*)
" Allocation error in Setup Cluster MKL, ja",ierr,nnz
246 if( .not.
associated(a) )
allocate(a(nnz+n), stat=ierr)
248 write(*,*)
" Allocation error in Setup Cluster MKL, a",ierr,nnz
251 if( .not.
associated(b) )
allocate(b(n), stat=ierr)
253 write(*,*)
" Allocation error in Setup Cluster MKL, b",ierr,n
256 if( .not.
associated(x) )
allocate(x(n), stat=ierr)
258 write(*,*)
" Allocation error in Setup Cluster MKL, x",ierr,n
262 allocate(ia(1),ja(1),a(1),b(1),x(1))
266 if (myrank==0 .and. spmat%timelog > 0 .and. debug > 0 ) &
267 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: - Allocate Matrix time(sec)=',t2-t1
274 if( myrank == 0 )
then
285 if (myrank==0 .and. spmat%timelog > 0 .and. debug > 0 ) &
286 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: - Gather Matrix time(sec)=',t3-t2
290 call coo2csr(n, nnz, ia, ja, a)
291 if( debug>0 )
call check_csr(n, nnz, ia, ja, a)
295 if (myrank==0 .and. spmat%timelog > 0 .and. debug > 0 ) &
296 write(*,
'(A,f10.3)')
' [Cluster Pardiso]: - Convert Matrix Format time(sec)=',t4-t3
298 deallocate(dispmat,ncounts)
301 subroutine coo2csr(n, nnz, ia, ja, a)
302 integer(kind=kint),
intent(in) :: n
303 integer(kind=kint),
intent(inout) :: nnz
304 integer(kind=kint),
pointer,
intent(inout) :: ia(:)
305 integer(kind=kint),
pointer,
intent(inout) :: ja(:)
306 real(kind=
kreal),
pointer,
intent(inout) :: a(:)
308 integer(kind=kint) :: i, k, idx, is, ie, nnz_new, prev
309 integer(kind=kint),
allocatable :: counter(:), counter_curr(:), ia0(:), ja0(:)
310 real(kind=
kreal),
allocatable :: a0(:)
312 allocate(counter(0:n),counter_curr(n),ia0(n+1),ja0(nnz),a0(nnz))
317 counter(idx) = counter(idx) + 1
321 ia0(i+1) = ia0(i)+counter(i)
327 counter_curr(idx) = counter_curr(idx) + 1
328 idx = ia0(idx)+counter_curr(idx)
338 call sort_column_ascending_order0(ie-is+1,ja0(is:ie),a0(is:ie))
345 ia(i+1) = ia(i)+counter(i)
352 do k=ia0(i)+1,ia0(i+1)
353 if( ja0(k) == prev )
then
354 a(nnz_new) = a(nnz_new) + a0(k)
356 nnz_new = nnz_new + 1
362 ia(i+1) = nnz_new + 1
368 deallocate(counter,counter_curr,ia0,ja0,a0)
372 subroutine check_csr(n, nnz, ia, ja, a)
373 integer(kind=kint),
intent(in) :: n
374 integer(kind=kint),
intent(in) :: nnz
375 integer(kind=kint),
pointer,
intent(in) :: ia(:)
376 integer(kind=kint),
pointer,
intent(in) :: ja(:)
377 real(kind=
kreal),
pointer,
intent(in) :: a(:)
381 if( ia(n+1)-1 /= nnz )
then
382 write(*,*)
"Error in check_csr(1): ia(n+1)-1 /= nnz"
388 if( ja(k) <= 0 )
then
389 write(*,*)
"Error in check_csr(2): ja(k) <= 0"
393 write(*,*)
"Error in check_csr(3): ja(k) > n"
397 if( ja(k) <= ja(k-1) )
then
398 write(*,*)
"Error in check_csr(4): ja(k) <= ja(k-1)"
407 subroutine sort_column_ascending_order(spMAT,myrank)
408 type (sparse_matrix),
intent(inout) :: spmat
409 integer(kind=kint),
intent(in) :: myrank
411 integer(kind=kint) :: i, is, ie
414 if(spmat%IRN(i)-spmat%IRN(i-1)>10000) &
415 &
write(*,*)
'warning: n may be too large for set_mkl_set_prof0 for row ',i-1
422 ie = spmat%IRN(i+1)-1
423 call sort_column_ascending_order0(ie-is+1,spmat%JCN(is:ie),spmat%A(is:ie))
427 end subroutine sort_column_ascending_order
429 subroutine sort_column_ascending_order0(n,ja,a)
430 integer(kind=kint),
intent(in) :: n
431 integer(kind=kint),
intent(inout) :: ja(:)
432 real(kind=
kreal),
intent(inout) :: a(:)
434 integer(kind=kint) :: i, work(2,n)
435 real(kind=
kreal) :: oa(n)
450 end subroutine sort_column_ascending_order0
452 recursive subroutine myqsort(n,work)
453 integer(kind=kint),
intent(in) :: n
454 integer(kind=kint),
intent(inout) :: work(:,:)
456 integer(kind=kint) :: i, key, work1(2,n), n_low, n_high, next, minidx, maxidx
466 if( work(1,i) > next ) sorted = .false.
467 if( next < minidx ) minidx = next
468 if( next > maxidx ) maxidx = next
472 if( n<5 .or. maxidx-minidx < 2)
then
473 call myinssort(n,work)
477 key = (minidx+maxidx)/2
481 if(work(1,i)<key)
then
483 work1(1:2,n_low) = work(1:2,i)
489 if(work(1,i)>=key)
then
491 work1(1:2,n_low+n_high) = work(1:2,i)
495 if(n_low>0)
call myqsort(n_low,work1(1:2,1:n_low))
496 if(n_high>0)
call myqsort(n_high,work1(1:2,n_low+1:n_low+n_high))
499 work(1:2,i) = work1(1:2,i)
502 end subroutine myqsort
504 subroutine myinssort(n,work)
505 integer(kind=kint),
intent(in) :: n
506 integer(kind=kint),
intent(inout) :: work(:,:)
508 integer(kind=kint) :: i,j,tmp(2)
512 if(work(1,j)<work(1,j-1))
then
513 tmp(1:2) = work(1:2,j)
514 work(1:2,j) = work(1:2,j-1)
515 work(1:2,j-1) = tmp(1:2)
520 end subroutine myinssort
integer(kind=kint) function hecmw_comm_get_size()
integer(kind=kint) function hecmw_comm_get_comm()
integer(kind=4), parameter kreal
integer(kind=kint) function hecmw_comm_get_rank()
subroutine hecmw_abort(comm)
real(kind=kreal) function hecmw_wtime()
This module provides linear equation solver interface for Cluster Pardiso.
subroutine, public hecmw_clustermkl_wrapper(spMAT, phase_start, solution, istat)
subroutine hecmw_gatherv_real(sbuf, sc, rbuf, rcs, disp, root, comm)
subroutine hecmw_allgather_int_1(sval, rbuf, comm)
subroutine hecmw_gatherv_int(sbuf, sc, rbuf, rcs, disp, root, comm)
This module provides DOF based sparse matrix data structure (CSR and COO)
integer(kind=kint), parameter, public sparse_matrix_type_coo
integer(kind=kint), parameter, public sparse_matrix_type_csr
integer(kind=kint), parameter, public sparse_matrix_symtype_sym
subroutine, public sparse_matrix_gather_rhs(spMAT, rhs_all)
subroutine, public sparse_matrix_scatter_rhs(spMAT, rhs_all)
integer(kind=kint), parameter, public sparse_matrix_symtype_asym
integer(kind=kint), parameter, public sparse_matrix_symtype_spd