27 integer(kind=kint) :: N
28 real(kind=
kreal),
pointer :: d(:) => null()
29 real(kind=
kreal),
pointer :: al(:) => null()
30 real(kind=
kreal),
pointer :: au(:) => null()
31 integer(kind=kint),
pointer :: indexL(:) => null()
32 integer(kind=kint),
pointer :: indexU(:) => null()
33 integer(kind=kint),
pointer :: itemL(:) => null()
34 integer(kind=kint),
pointer :: itemU(:) => null()
35 real(kind=
kreal),
pointer :: alu(:) => null()
37 integer(kind=kint) :: NColor
38 integer(kind=kint),
pointer :: COLORindex(:) => null()
39 integer(kind=kint),
pointer :: perm(:) => null()
40 integer(kind=kint),
pointer :: iperm(:) => null()
42 logical,
save :: isFirst = .true.
44 logical,
save :: INITIALIZED = .false.
50 real(kind=
kreal),
save :: omega = 1.d0
53 integer(kind=kint),
parameter :: numOfBlockPerThread = 100
54 integer(kind=kint),
save :: numOfThread = 1, numofblock
55 integer(kind=kint),
save,
allocatable :: icToBlockIndex(:)
56 integer(kind=kint),
save,
allocatable :: blockIndexToColorIndex(:)
57 integer(kind=kint),
save :: sectorCacheSize0, sectorCacheSize1
59 integer(kind=kint),
parameter :: SSOR_DEBUG = 0
66 integer(kind=kint ) :: npl, npu
67 integer(kind=kint ) :: ncolor_in
68 real (kind=
kreal) :: omega_inv
69 real (kind=
kreal) :: alutmp(4,4), pw(4)
70 integer(kind=kint ) :: ii, i, j, k
71 integer(kind=kint ) :: nthreads = 1
72 integer(kind=kint ),
allocatable :: perm_tmp(:)
73 real (kind=
kreal) :: t0
75 if (ssor_debug >= 1)
then
77 write(*,*)
'DEBUG: SSOR setup start',
hecmw_wtime()-t0
81 if (hecmat%Iarray(98) == 1)
then
83 else if (hecmat%Iarray(97) == 1)
then
98 omega_inv = 1.d0 / omega
101 allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
103 hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
104 if (ssor_debug >= 1)
write(*,*)
'DEBUG: RCM ordering done',
hecmw_wtime()-t0
106 hecmat%indexU, hecmat%itemU, perm_tmp, &
107 ncolor_in, ncolor, colorindex, perm, iperm)
108 if (ssor_debug >= 1)
write(*,*)
'DEBUG: MC ordering done',
hecmw_wtime()-t0
112 if (nthreads == 1)
then
114 allocate(colorindex(0:1), perm(n), iperm(n))
122 allocate(colorindex(0:n), perm_tmp(n), perm(n), iperm(n))
124 hecmat%indexU, hecmat%itemU, perm_tmp, iperm)
125 if (ssor_debug >= 1)
write(*,*)
'DEBUG: RCM ordering done',
hecmw_wtime()-t0
127 hecmat%indexU, hecmat%itemU, perm_tmp, &
128 ncolor_in, ncolor, colorindex, perm, iperm)
129 if (ssor_debug >= 1)
write(*,*)
'DEBUG: MC ordering done',
hecmw_wtime()-t0
137 do j=hecmat%indexU(i-1)+1,hecmat%indexU(i)
138 if( hecmat%itemU(j) > n )
exit
142 npl = max(hecmat%indexL(n),npl)
143 npu = hecmat%indexU(n)
144 allocate(indexl(0:n), indexu(0:n), iteml(npl), itemu(npu))
146 hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
147 indexl, indexu, iteml, itemu)
148 if (ssor_debug >= 1)
write(*,*)
'DEBUG: reordering profile done',
hecmw_wtime()-t0
151 allocate(d(16*n), al(16*npl), au(16*npu))
153 hecmat%indexL, hecmat%indexU, hecmat%itemL, hecmat%itemU, &
154 hecmat%AL, hecmat%AU, hecmat%D, &
155 indexl, indexu, iteml, itemu, al, au, d)
156 if (ssor_debug >= 1)
write(*,*)
'DEBUG: reordering values done',
hecmw_wtime()-t0
176 alutmp(1,1)= alu(16*ii-15) * omega_inv
177 alutmp(1,2)= alu(16*ii-14) * omega_inv
178 alutmp(1,3)= alu(16*ii-13) * omega_inv
179 alutmp(1,4)= alu(16*ii-12) * omega_inv
180 alutmp(2,1)= alu(16*ii-11) * omega_inv
181 alutmp(2,2)= alu(16*ii-10) * omega_inv
182 alutmp(2,3)= alu(16*ii- 9) * omega_inv
183 alutmp(2,4)= alu(16*ii- 8) * omega_inv
184 alutmp(3,1)= alu(16*ii- 7) * omega_inv
185 alutmp(3,2)= alu(16*ii- 6) * omega_inv
186 alutmp(3,3)= alu(16*ii- 5) * omega_inv
187 alutmp(3,4)= alu(16*ii- 4) * omega_inv
188 alutmp(4,1)= alu(16*ii- 3) * omega_inv
189 alutmp(4,2)= alu(16*ii- 2) * omega_inv
190 alutmp(4,3)= alu(16*ii- 1) * omega_inv
191 alutmp(4,4)= alu(16*ii ) * omega_inv
194 alutmp(k,k)= 1.d0/alutmp(k,k)
196 alutmp(i,k)= alutmp(i,k) * alutmp(k,k)
198 pw(j)= alutmp(i,j) - alutmp(i,k)*alutmp(k,j)
205 alu(16*ii-15)= alutmp(1,1)
206 alu(16*ii-14)= alutmp(1,2)
207 alu(16*ii-13)= alutmp(1,3)
208 alu(16*ii-12)= alutmp(1,4)
209 alu(16*ii-11)= alutmp(2,1)
210 alu(16*ii-10)= alutmp(2,2)
211 alu(16*ii- 9)= alutmp(2,3)
212 alu(16*ii- 8)= alutmp(2,4)
213 alu(16*ii- 7)= alutmp(3,1)
214 alu(16*ii- 6)= alutmp(3,2)
215 alu(16*ii- 5)= alutmp(3,3)
216 alu(16*ii- 4)= alutmp(3,4)
217 alu(16*ii- 3)= alutmp(4,1)
218 alu(16*ii- 2)= alutmp(4,2)
219 alu(16*ii- 1)= alutmp(4,3)
220 alu(16*ii )= alutmp(4,4)
233 hecmat%Iarray(98) = 0
234 hecmat%Iarray(97) = 0
236 if (ssor_debug >= 1)
write(*,*)
'DEBUG: SSOR setup done',
hecmw_wtime()-t0
240 subroutine setup_tuning_parameters
243 integer(kind=kint) :: blockindex, elementcount, numofelement, ii
244 real(kind=
kreal) :: numofelementperblock
245 integer(kind=kint) :: ic, i
246 if (ssor_debug >= 1)
write(*,*)
'DEBUG: setting up tuning parameters for SSOR'
251 numofblock = numofthread * numofblockperthread
252 if (
allocated(ictoblockindex))
deallocate(ictoblockindex)
253 if (
allocated(blockindextocolorindex))
deallocate(blockindextocolorindex)
254 allocate (ictoblockindex(0:ncolor), &
255 blockindextocolorindex(0:numofblock + ncolor))
256 numofelement = n + indexl(n) + indexu(n)
257 numofelementperblock = dble(numofelement) / numofblock
260 ictoblockindex(0) = 0
261 blockindextocolorindex = -1
262 blockindextocolorindex(0) = 0
266 do i = colorindex(ic-1)+1, colorindex(ic)
267 elementcount = elementcount + 1
268 elementcount = elementcount + (indexl(i) - indexl(i-1))
269 elementcount = elementcount + (indexu(i) - indexu(i-1))
270 if (elementcount > ii * numofelementperblock &
271 .or. i == colorindex(ic))
then
273 blockindex = blockindex + 1
274 blockindextocolorindex(blockindex) = i
277 ictoblockindex(ic) = blockindex
279 numofblock = blockindex
282 sectorcachesize0, sectorcachesize1 )
283 end subroutine setup_tuning_parameters
287 real(kind=
kreal),
intent(inout) :: zp(:)
288 integer(kind=kint) :: ic, i, iold, j, isl, iel, isu, ieu, k
289 real(kind=
kreal) :: sw1, sw2, sw3, sw4, x1, x2, x3, x4
291 real(kind=
kreal) :: omega_fac
294 integer(kind=kint) :: blockindex
298 call setup_tuning_parameters
304 omega_fac = 2.d0 - omega
323 do i = colorindex(ic-1)+1, colorindex(ic)
326 do blockindex = ictoblockindex(ic-1)+1, ictoblockindex(ic)
327 do i = blockindextocolorindex(blockindex-1)+1, &
328 blockindextocolorindex(blockindex)
332 sw1= omega_fac * zp(4*iold-3)
333 sw2= omega_fac * zp(4*iold-2)
334 sw3= omega_fac * zp(4*iold-1)
335 sw4= omega_fac * zp(4*iold )
345 sw1= sw1 - al(16*j-15)*x1 - al(16*j-14)*x2 - al(16*j-13)*x3 - al(16*j-12)*x4
346 sw2= sw2 - al(16*j-11)*x1 - al(16*j-10)*x2 - al(16*j- 9)*x3 - al(16*j- 8)*x4
347 sw3= sw3 - al(16*j- 7)*x1 - al(16*j- 6)*x2 - al(16*j- 5)*x3 - al(16*j- 4)*x4
348 sw4= sw4 - al(16*j- 3)*x1 - al(16*j- 2)*x2 - al(16*j- 1)*x3 - al(16*j- 0)*x4
355 x2= x2 - alu(16*i-11)*x1
356 x3= x3 - alu(16*i- 7)*x1 - alu(16*i- 6)*x2
357 x4= x4 - alu(16*i- 3)*x1 - alu(16*i- 2)*x2 - alu(16*i- 1)*x3
360 x3= alu(16*i- 5)*( x3 - alu(16*i- 4)*x4 )
361 x2= alu(16*i-10)*( x2 - alu(16*i- 8)*x4 - alu(16*i- 9)*x3 )
362 x1= alu(16*i-15)*( x1 - alu(16*i-12)*x4 - alu(16*i-13)*x3 - alu(16*i-14)*x2 )
383 do i = colorindex(ic-1)+1, colorindex(ic)
386 do blockindex = ictoblockindex(ic), ictoblockindex(ic-1)+1, -1
387 do i = blockindextocolorindex(blockindex), &
388 blockindextocolorindex(blockindex-1)+1, -1
407 sw1= sw1 + au(16*j-15)*x1 + au(16*j-14)*x2 + au(16*j-13)*x3 + au(16*j-12)*x4
408 sw2= sw2 + au(16*j-11)*x1 + au(16*j-10)*x2 + au(16*j- 9)*x3 + au(16*j- 8)*x4
409 sw3= sw3 + au(16*j- 7)*x1 + au(16*j- 6)*x2 + au(16*j- 5)*x3 + au(16*j- 4)*x4
410 sw4= sw4 + au(16*j- 3)*x1 + au(16*j- 2)*x2 + au(16*j- 1)*x3 + au(16*j- 0)*x4
418 x2= x2 - alu(16*i-11)*x1
419 x3= x3 - alu(16*i- 7)*x1 - alu(16*i- 6)*x2
420 x4= x4 - alu(16*i- 3)*x1 - alu(16*i- 2)*x2 - alu(16*i- 1)*x3
422 x3= alu(16*i- 5)*( x3 - alu(16*i- 4)*x4 )
423 x2= alu(16*i-10)*( x2 - alu(16*i- 8)*x4 - alu(16*i- 9)*x3 )
424 x1= alu(16*i-15)*( x1 - alu(16*i-12)*x4 - alu(16*i-13)*x3 - alu(16*i-14)*x2 )
426 zp(4*iold-3)= zp(4*iold-3) - x1
427 zp(4*iold-2)= zp(4*iold-2) - x2
428 zp(4*iold-1)= zp(4*iold-1) - x3
429 zp(4*iold )= zp(4*iold ) - x4
453 if (
associated(colorindex))
deallocate(colorindex)
454 if (
associated(perm))
deallocate(perm)
455 if (
associated(iperm))
deallocate(iperm)
456 if (
associated(alu))
deallocate(alu)
457 if (
associated(d))
deallocate(d)
458 if (
associated(al))
deallocate(al)
459 if (
associated(au))
deallocate(au)
460 if (
associated(indexl))
deallocate(indexl)
461 if (
associated(indexu))
deallocate(indexu)
462 if (
associated(iteml))
deallocate(iteml)
463 if (
associated(itemu))
deallocate(itemu)
475 initialized = .false.
real(kind=kreal) function, public hecmw_mat_get_omega(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_ncolor_in(hecMAT)
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_renum_item(N, perm, indexXp, itemXp)
subroutine, public hecmw_matrix_reorder_profile(N, perm, iperm, indexL, indexU, itemL, itemU, indexLp, indexUp, itemLp, itemUp)
subroutine, public hecmw_precond_ssor_44_clear(hecMAT)
subroutine, public hecmw_precond_ssor_44_apply(ZP)
subroutine, public hecmw_precond_ssor_44_setup(hecMAT)
subroutine, public hecmw_tuning_fx_calc_sector_cache(N, NDOF, sectorCacheSize0, sectorCacheSize1)
integer(kind=4), parameter kreal
real(kind=kreal) function hecmw_wtime()
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)