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(2,2), pw(2)
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(4*n), al(4*npl), au(4*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(4*n), cal(4*npcl), cau(4*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(4*ii-3) = alu(4*ii-3) + hecmat%cmat%pair(k)%val(1, 1)
167 alu(4*ii-2) = alu(4*ii-2) + hecmat%cmat%pair(k)%val(1, 2)
168 alu(4*ii-1) = alu(4*ii-1) + hecmat%cmat%pair(k)%val(2, 1)
169 alu(4*ii-0) = alu(4*ii-0) + hecmat%cmat%pair(k)%val(2, 2)
176 alutmp(1,1)= alu(4*ii-3) * sigma_diag
177 alutmp(1,2)= alu(4*ii-2)
178 alutmp(2,1)= alu(4*ii-1)
179 alutmp(2,2)= alu(4*ii-0) * sigma_diag
181 alutmp(k,k)= 1.d0/alutmp(k,k)
183 alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
185 pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
192 alu(4*ii-3)= alutmp(1,1)
193 alu(4*ii-2)= alutmp(1,2)
194 alu(4*ii-1)= alutmp(2,1)
195 alu(4*ii-0)= alutmp(2,2)
203 hecmat%Iarray(98) = 0
204 hecmat%Iarray(97) = 0
213 real(kind=
kreal),
intent(inout) :: zp(:)
214 integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k
215 real(kind=
kreal) :: sw(2), x(2)
218 integer(kind=kint),
parameter :: numofblockperthread = 100
219 integer(kind=kint),
save :: numofthread = 1, numofblock
220 integer(kind=kint),
save,
allocatable :: ictoblockindex(:)
221 integer(kind=kint),
save,
allocatable :: blockindextocolorindex(:)
222 integer(kind=kint),
save :: sectorcachesize0, sectorcachesize1
223 integer(kind=kint) :: blockindex, elementcount, numofelement, ii
224 real(kind=
kreal) :: numofelementperblock
225 integer(kind=kint) :: my_rank
229 numofblock = numofthread * numofblockperthread
230 if (
allocated(ictoblockindex))
deallocate(ictoblockindex)
231 if (
allocated(blockindextocolorindex))
deallocate(blockindextocolorindex)
232 allocate (ictoblockindex(0:ncolor), &
233 blockindextocolorindex(0:numofblock + ncolor))
234 numofelement = n + indexl(n) + indexu(n)
235 numofelementperblock = dble(numofelement) / numofblock
238 ictoblockindex(0) = 0
239 blockindextocolorindex = -1
240 blockindextocolorindex(0) = 0
249 do i = colorindex(ic-1)+1, colorindex(ic)
250 elementcount = elementcount + 1
251 elementcount = elementcount + (indexl(i) - indexl(i-1))
252 elementcount = elementcount + (indexu(i) - indexu(i-1))
253 if (elementcount > ii * numofelementperblock &
254 .or. i == colorindex(ic))
then
256 blockindex = blockindex + 1
257 blockindextocolorindex(blockindex) = i
262 ictoblockindex(ic) = blockindex
264 numofblock = blockindex
267 sectorcachesize0, sectorcachesize1 )
287 do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
288 do i = blockindextocolorindex(blockindex-1)+1, &
289 blockindextocolorindex(blockindex)
299 sw(1)= sw(1) - al(4*j-3)*x(1) - al(4*j-2)*x(2)
300 sw(2)= sw(2) - al(4*j-1)*x(1) - al(4*j-0)*x(2)
303 if (ncontact.ne.0)
then
310 sw(1)= sw(1) - cal(4*j-3)*x(1) - cal(4*j-2)*x(2)
311 sw(2)= sw(2) - cal(4*j-1)*x(1) - cal(4*j-0)*x(2)
316 x(2)= x(2) - alu(4*i-1)*x(1)
317 x(2)= alu(4*i )* x(2)
318 x(1)= alu(4*i-3)*( x(1) - alu(4*i-2)*x(2))
329 do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
330 do i = blockindextocolorindex(blockindex), &
331 blockindextocolorindex(blockindex-1)+1, -1
343 sw(1)= sw(1) + au(4*j-3)*x(1) + au(4*j-2)*x(2)
344 sw(2)= sw(2) + au(4*j-1)*x(1) + au(4*j-0)*x(2)
347 if (ncontact.gt.0)
then
348 isu= indexcu(i-1) + 1
354 sw(1)= sw(1) + cau(4*j-3)*x(1) + cau(4*j-2)*x(2)
355 sw(2)= sw(2) + cau(4*j-1)*x(1) + cau(4*j-0)*x(2)
360 x(2)= x(2) - alu(4*i-1)*x(1)
363 x(2)= alu(4*i )* x(2)
364 x(1)= alu(4*i-3)*( x(1) - alu(4*i-2)*x(2) )
367 zp(2*iold-1)= zp(2*iold-1) - x(1)
368 zp(2*iold )= zp(2*iold ) - x(2)
385 integer(kind=kint ) :: nthreads = 1
387 if (
associated(colorindex))
deallocate(colorindex)
388 if (
associated(perm))
deallocate(perm)
389 if (
associated(iperm))
deallocate(iperm)
390 if (
associated(alu))
deallocate(alu)
391 if (nthreads >= 1)
then
392 if (
associated(d))
deallocate(d)
393 if (
associated(al))
deallocate(al)
394 if (
associated(au))
deallocate(au)
395 if (
associated(indexl))
deallocate(indexl)
396 if (
associated(indexu))
deallocate(indexu)
397 if (
associated(iteml))
deallocate(iteml)
398 if (
associated(itemu))
deallocate(itemu)
399 if (ncontact.ne.0)
then
400 if (
associated(cal))
deallocate(cal)
401 if (
associated(cau))
deallocate(cau)
402 if (
associated(indexcl))
deallocate(indexcl)
403 if (
associated(indexcu))
deallocate(indexcu)
404 if (
associated(itemcl))
deallocate(itemcl)
405 if (
associated(itemcu))
deallocate(itemcu)
419 if (ncontact.ne.0)
then
429 initialized = .false.
432 subroutine write_debug_info
434 integer(kind=kint) :: my_rank, ic, in
437 if (my_rank.eq.0)
then
438 write(*,*)
'DEBUG: Output fort.19000+myrank and fort.29000+myrank for coloring information'
440 write(19000+my_rank,
'(a)')
'#NCOLORTot'
441 write(19000+my_rank,*) ncolor
442 write(19000+my_rank,
'(a)')
'#ic COLORindex(ic-1)+1 COLORindex(ic)'
444 write(19000+my_rank,*) ic, colorindex(ic-1)+1,colorindex(ic)
446 write(29000+my_rank,
'(a)')
'#n_node'
447 write(29000+my_rank,*) n
448 write(29000+my_rank,
'(a)')
'#in OLDtoNEW(in) NEWtoOLD(in)'
450 write(29000+my_rank,*) in, iperm(in), perm(in)
451 if (perm(iperm(in)) .ne. in)
then
452 write(29000+my_rank,*)
'** WARNING **: NEWtoOLD and OLDtoNEW: ',in
455 end subroutine write_debug_info
457 subroutine check_ordering
459 integer(kind=kint) :: ic, i, j, k
460 integer(kind=kint),
allocatable :: iicolor(:)
462 if (ncolor.gt.1)
then
465 do i= colorindex(ic-1)+1, colorindex(ic)
471 do i= colorindex(ic-1)+1, colorindex(ic)
472 do j= indexl(i-1)+1, indexl(i)
474 if (iicolor(i).eq.iicolor(k))
then
475 write(*,*) .eq.
'** ERROR **: L-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
482 do i= colorindex(ic), colorindex(ic-1)+1, -1
483 do j= indexu(i-1)+1, indexu(i)
485 if (iicolor(i).eq.iicolor(k))
then
486 write(*,*) .eq.
'** ERROR **: U-part: iicolor(i)iicolor(k)',i,k,iicolor(i)
494 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_22_clear(hecmat)
subroutine, public hecmw_precond_ssor_22_setup(hecmat)
subroutine, public hecmw_precond_ssor_22_apply(zp)
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)