26 integer(kind=kint) :: N
27 real(kind=
kreal),
pointer :: d(:) => null()
28 real(kind=
kreal),
pointer :: al(:) => null()
29 real(kind=
kreal),
pointer :: au(:) => null()
30 integer(kind=kint),
pointer :: indexL(:) => null()
31 integer(kind=kint),
pointer :: indexU(:) => null()
32 integer(kind=kint),
pointer :: itemL(:) => null()
33 integer(kind=kint),
pointer :: itemU(:) => null()
34 real(kind=
kreal),
pointer :: alu(:) => null()
36 integer(kind=kint) :: NContact = 0
37 real(kind=
kreal),
pointer :: cal(:) => null()
38 real(kind=
kreal),
pointer :: cau(:) => null()
39 integer(kind=kint),
pointer :: indexCL(:) => null()
40 integer(kind=kint),
pointer :: indexCU(:) => null()
41 integer(kind=kint),
pointer :: itemCL(:) => null()
42 integer(kind=kint),
pointer :: itemCU(:) => null()
44 integer(kind=kint) :: NColor
45 integer(kind=kint),
pointer :: COLORindex(:) => null()
46 integer(kind=kint),
pointer :: perm(:) => null()
47 integer(kind=kint),
pointer :: iperm(:) => null()
49 logical,
save :: isFirst = .true.
51 logical,
save :: INITIALIZED = .false.
58 integer(kind=kint ) :: npl, npu, npcl, npcu
59 real (kind=
kreal),
allocatable :: cd(:)
60 integer(kind=kint ) :: ncolor_in
61 real (kind=
kreal) :: sigma_diag
62 real (kind=
kreal) :: alutmp(1,1), pw(1)
63 integer(kind=kint ) :: ii, i, j, k
64 integer(kind=kint ) :: nthreads = 1
65 integer(kind=kint ),
allocatable :: perm_tmp(:)
72 if (hecmat%Iarray(98) == 1)
then
74 else if (hecmat%Iarray(97) == 1)
then
87 ncontact = hecmat%cmat%n_val
89 if (ncontact.gt.0)
then
93 if (nthreads == 1)
then
95 allocate(colorindex(0:1), perm(n), iperm(n))
103 allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
105 hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
108 hecmat%indexU, hecmat%itemU, perm_tmp, &
109 ncolor_in, ncolor, colorindex, perm, iperm)
116 npl = hecmat%indexL(n)
117 npu = hecmat%indexU(n)
118 allocate(indexl(0:n), indexu(0:n), iteml(npl), itemu(npu))
120 hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
121 indexl, indexu, iteml, itemu)
126 allocate(d(n), al(npl), au(npu))
128 hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
129 hecmat%AL, hecmat%AU, hecmat%D, &
130 indexl, indexu, iteml, itemu, al, au, d)
136 if (ncontact.gt.0)
then
137 npcl = hecmat%indexCL(n)
138 npcu = hecmat%indexCU(n)
139 allocate(indexcl(0:n), indexcu(0:n), itemcl(npcl), itemcu(npcu))
141 hecmat%indexCL, hecmat%indexCU, hecmat%itemCL, hecmat%itemCU, &
142 indexcl, indexcu, itemcl, itemcu)
144 allocate(cd(n), cal(npcl), cau(npcu))
146 hecmat%indexCL, hecmat%indexCU, hecmat%itemCL, hecmat%itemCU, &
147 hecmat%CAL, hecmat%CAU, hecmat%D, &
148 indexcl, indexcu, itemcl, itemcu, cal, cau, cd)
162 if (ncontact.gt.0)
then
163 do k= 1, hecmat%cmat%n_val
164 if (hecmat%cmat%pair(k)%i.ne.hecmat%cmat%pair(k)%j) cycle
165 ii = iperm( hecmat%cmat%pair(k)%i )
166 alu(ii) = alu(ii) + hecmat%cmat%pair(k)%val(1, 1)
173 alutmp(1,1)= alu(ii) * sigma_diag
174 alutmp(1,1)= 1.d0/alutmp(1,1)
183 hecmat%Iarray(98) = 0
184 hecmat%Iarray(97) = 0
193 real(kind=
kreal),
intent(inout) :: zp(:)
194 integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k
195 real(kind=
kreal) :: sw(1), x(1)
198 integer(kind=kint),
parameter :: numofblockperthread = 100
199 integer(kind=kint),
save :: numofthread = 1, numofblock
200 integer(kind=kint),
save,
allocatable :: ictoblockindex(:)
201 integer(kind=kint),
save,
allocatable :: blockindextocolorindex(:)
202 integer(kind=kint),
save :: sectorcachesize0, sectorcachesize1
203 integer(kind=kint) :: blockindex, elementcount, numofelement, ii
204 real(kind=
kreal) :: numofelementperblock
205 integer(kind=kint) :: my_rank
209 numofblock = numofthread * numofblockperthread
210 if (
allocated(ictoblockindex))
deallocate(ictoblockindex)
211 if (
allocated(blockindextocolorindex))
deallocate(blockindextocolorindex)
212 allocate (ictoblockindex(0:ncolor), &
213 blockindextocolorindex(0:numofblock + ncolor))
214 numofelement = n + indexl(n) + indexu(n)
215 numofelementperblock = dble(numofelement) / numofblock
218 ictoblockindex(0) = 0
219 blockindextocolorindex = -1
220 blockindextocolorindex(0) = 0
229 do i = colorindex(ic-1)+1, colorindex(ic)
230 elementcount = elementcount + 1
231 elementcount = elementcount + (indexl(i) - indexl(i-1))
232 elementcount = elementcount + (indexu(i) - indexu(i-1))
233 if (elementcount > ii * numofelementperblock &
234 .or. i == colorindex(ic))
then
236 blockindex = blockindex + 1
237 blockindextocolorindex(blockindex) = i
242 ictoblockindex(ic) = blockindex
244 numofblock = blockindex
247 sectorcachesize0, sectorcachesize1 )
267 do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
268 do i = blockindextocolorindex(blockindex-1)+1, &
269 blockindextocolorindex(blockindex)
277 sw(1)= sw(1) - al(j)*x(1)
280 if (ncontact.ne.0)
then
286 sw(1)= sw(1) - cal(j)*x(1)
301 do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
302 do i = blockindextocolorindex(blockindex), &
303 blockindextocolorindex(blockindex-1)+1, -1
314 sw(1)= sw(1) + au(j)*x(1)
317 if (ncontact.gt.0)
then
318 isu= indexcu(i-1) + 1
323 sw(1)= sw(1) + cau(j)*x(1)
331 zp(iold)= zp(iold) - x(1)
348 integer(kind=kint ) :: nthreads = 1
350 if (
associated(colorindex))
deallocate(colorindex)
351 if (
associated(perm))
deallocate(perm)
352 if (
associated(iperm))
deallocate(iperm)
353 if (
associated(alu))
deallocate(alu)
354 if (nthreads >= 1)
then
355 if (
associated(d))
deallocate(d)
356 if (
associated(al))
deallocate(al)
357 if (
associated(au))
deallocate(au)
358 if (
associated(indexl))
deallocate(indexl)
359 if (
associated(indexu))
deallocate(indexu)
360 if (
associated(iteml))
deallocate(iteml)
361 if (
associated(itemu))
deallocate(itemu)
362 if (ncontact.ne.0)
then
363 if (
associated(cal))
deallocate(cal)
364 if (
associated(cau))
deallocate(cau)
365 if (
associated(indexcl))
deallocate(indexcl)
366 if (
associated(indexcu))
deallocate(indexcu)
367 if (
associated(itemcl))
deallocate(itemcl)
368 if (
associated(itemcu))
deallocate(itemcu)
382 if (ncontact.ne.0)
then
392 initialized = .false.
395 subroutine write_debug_info
397 integer(kind=kint) :: my_rank, ic, in
400 if (my_rank.eq.0)
then
401 write(*,*)
'DEBUG: Output fort.19000+myrank and fort.29000+myrank for coloring information'
403 write(19000+my_rank,
'(a)')
'#NCOLORTot'
404 write(19000+my_rank,*) ncolor
405 write(19000+my_rank,
'(a)')
'#ic COLORindex(ic-1)+1 COLORindex(ic)'
407 write(19000+my_rank,*) ic, colorindex(ic-1)+1,colorindex(ic)
409 write(29000+my_rank,
'(a)')
'#n_node'
410 write(29000+my_rank,*) n
411 write(29000+my_rank,
'(a)')
'#in OLDtoNEW(in) NEWtoOLD(in)'
413 write(29000+my_rank,*) in, iperm(in), perm(in)
414 if (perm(iperm(in)) .ne. in)
then
415 write(29000+my_rank,*)
'** WARNING **: NEWtoOLD and OLDtoNEW: ',in
418 end subroutine write_debug_info
420 subroutine check_ordering
422 integer(kind=kint) :: ic, i, j, k
423 integer(kind=kint),
allocatable :: iicolor(:)
425 if (ncolor.gt.1)
then
428 do i= colorindex(ic-1)+1, colorindex(ic)
434 do i= colorindex(ic-1)+1, colorindex(ic)
435 do j= indexl(i-1)+1, indexl(i)
437 if (iicolor(i).eq.iicolor(k))
then
438 write(*,*) .eq.
'** ERROR **: L-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
445 do i= colorindex(ic), colorindex(ic-1)+1, -1
446 do j= indexu(i-1)+1, indexu(i)
448 if (iicolor(i).eq.iicolor(k))
then
449 write(*,*) .eq.
'** ERROR **: U-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
457 end subroutine check_ordering
integer(kind=kint) function, public hecmw_mat_get_ncolor_in(hecmat)
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecmat)
subroutine, public hecmw_matrix_reorder_renum_item(n, perm, indexxp, itemxp)
subroutine, public hecmw_matrix_reorder_values(n, ndof, perm, iperm, indexl, indexu, iteml, itemu, al, au, d, indexlp, indexup, itemlp, itemup, alp, aup, dp)
subroutine, public hecmw_matrix_reorder_profile(n, perm, iperm, indexl, indexu, iteml, itemu, indexlp, indexup, itemlp, itemup)
subroutine, public hecmw_precond_ssor_11_apply(zp)
subroutine, public hecmw_precond_ssor_11_setup(hecmat)
subroutine, public hecmw_precond_ssor_11_clear(hecmat)
subroutine, public hecmw_tuning_fx_calc_sector_cache(n, ndof, sectorcachesize0, sectorcachesize1)
integer(kind=4), parameter kreal
integer(kind=kint) function hecmw_comm_get_rank()
subroutine, public hecmw_matrix_ordering_rcm(n, indexl, iteml, indexu, itemu, perm, iperm)
subroutine, public hecmw_matrix_ordering_mc(n, indexl, iteml, indexu, itemu, perm_cur, ncolor_in, ncolor_out, colorindex, perm, iperm)